diff --git a/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl new file mode 100644 index 00000000000..269f40924c3 --- /dev/null +++ b/examples/p4est_3d_dgsem/elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl @@ -0,0 +1,76 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +boundary_condition = BoundaryConditionDirichlet(initial_condition) +boundary_conditions = (; Bottom = boundary_condition, + Top = boundary_condition, + Circle = boundary_condition, + Cut = boundary_condition) + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 4 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + local_twosided_variables_cons = [], + local_onesided_variables_nonlinear = [], + max_iterations_newton = 10) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +# Unstructured 3D half circle mesh from HOHQMesh +mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/11461efbfb02c42e06aca338b3d0b645/raw/81deeb1ebc4945952c30af5bb75fe222a18d975c/abaqus_half_circle_3d.inp", + joinpath(@__DIR__, "abaqus_half_circle_3d.inp")) + +mesh = P4estMesh{3}(mesh_file, initial_refinement_level = 0) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_conditions) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_integrals = (entropy,)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 1, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim, + extra_node_variables = (:limiting_coefficient,)) + +stepsize_callback = StepsizeCallback(cfl = 0.5) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); diff --git a/examples/tree_3d_dgsem/elixir_euler_convergence_sc_subcell.jl b/examples/tree_3d_dgsem/elixir_euler_convergence_sc_subcell.jl new file mode 100644 index 00000000000..fea1afea654 --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_euler_convergence_sc_subcell.jl @@ -0,0 +1,70 @@ +using OrdinaryDiffEqLowStorageRK +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +initial_condition = initial_condition_convergence_test + +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +polydeg = 3 +basis = LobattoLegendreBasis(polydeg) +limiter_idp = SubcellLimiterIDP(equations, basis; + positivity_variables_cons = ["rho"], + positivity_variables_nonlinear = [pressure], + local_twosided_variables_cons = [], + local_onesided_variables_nonlinear = [], + max_iterations_newton = 40, # Default parameters are not sufficient to fulfill bounds properly. + newton_tolerances = (1.0e-14, 1.0e-15)) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (0.0, 0.0, 0.0) +coordinates_max = (2.0, 2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 2, + n_cells_max = 10_000, + periodicity = true) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_convergence_test, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 1.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); diff --git a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl index ead0b003e81..61452ce3174 100644 --- a/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl +++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -59,11 +59,11 @@ coordinates_max = (1.0, 1.0, 1.0) mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 3, n_cells_max = 100_000, - periodicity = true) + periodicity = false) # create the semi discretization object semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - boundary_conditions = boundary_condition_periodic) + boundary_conditions = BoundaryConditionDirichlet(initial_condition)) ############################################################################### # ODE solvers, callbacks etc. diff --git a/src/equations/ideal_glm_mhd_3d.jl b/src/equations/ideal_glm_mhd_3d.jl index f6cff2c9851..4b8348f087e 100644 --- a/src/equations/ideal_glm_mhd_3d.jl +++ b/src/equations/ideal_glm_mhd_3d.jl @@ -364,6 +364,9 @@ end # For `VolumeIntegralSubcellLimiting` the nonconservative flux is created as a callable struct to # enable dispatch on the type of the nonconservative term (symmetric / jump). """ + flux_nonconservative_powell_local_symmetric(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations3D) flux_nonconservative_powell_local_symmetric(u_ll, u_rr, normal_direction::AbstractVector, equations::IdealGlmMhdEquations3D) @@ -389,6 +392,58 @@ compute the subcell fluxes in dg_3d_subcell_limiters.jl. [DOI: 10.1016/j.jcp.2023.112607](https://doi.org/10.1016/j.jcp.2023.112607). [arXiv: 2211.14009](https://doi.org/10.48550/arXiv.2211.14009). """ +@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations3D) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + if orientation == 1 + B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_ll * B1_avg, + B2_ll * B1_avg, + B3_ll * B1_avg, + v_dot_B_ll * B1_avg + v1_ll * psi_ll * psi_avg, + v1_ll * B1_avg, + v2_ll * B1_avg, + v3_ll * B1_avg, + v1_ll * psi_avg) + elseif orientation == 2 + B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_ll * B2_avg, + B2_ll * B2_avg, + B3_ll * B2_avg, + v_dot_B_ll * B2_avg + v2_ll * psi_ll * psi_avg, + v1_ll * B2_avg, + v2_ll * B2_avg, + v3_ll * B2_avg, + v2_ll * psi_avg) + else # orientation == 3 + B3_avg = (B3_ll + B3_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_ll * B3_avg, + B2_ll * B3_avg, + B3_ll * B3_avg, + v_dot_B_ll * B3_avg + v3_ll * psi_ll * psi_avg, + v1_ll * B3_avg, + v2_ll * B3_avg, + v3_ll * B3_avg, + v3_ll * psi_avg) + end + + return f +end + @inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr, normal_direction::AbstractVector, equations::IdealGlmMhdEquations3D) @@ -428,6 +483,10 @@ compute the subcell fluxes in dg_3d_subcell_limiters.jl. end """ + flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations3D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) flux_nonconservative_powell_local_symmetric(u_ll, normal_direction_ll::AbstractVector, equations::IdealGlmMhdEquations3D, nonconservative_type::NonConservativeLocal, @@ -442,6 +501,68 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl. """ +@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, + orientation::Integer, + equations::IdealGlmMhdEquations3D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + v1_ll = rho_v1_ll / rho_ll + v2_ll = rho_v2_ll / rho_ll + v3_ll = rho_v3_ll / rho_ll + v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll + f = SVector(0, + B1_ll, + B2_ll, + B3_ll, + v_dot_B_ll, + v1_ll, + v2_ll, + v3_ll, + 0) + else #nonconservative_term ==2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + if orientation == 1 + v1_ll = rho_v1_ll / rho_ll + f = SVector(0, + 0, + 0, + 0, + v1_ll * psi_ll, + 0, + 0, + 0, + v1_ll) + elseif orientation == 2 + v2_ll = rho_v2_ll / rho_ll + f = SVector(0, + 0, + 0, + 0, + v2_ll * psi_ll, + 0, + 0, + 0, + v2_ll) + else #orientation == 3 + v3_ll = rho_v3_ll / rho_ll + f = SVector(0, + 0, + 0, + 0, + v3_ll * psi_ll, + 0, + 0, + 0, + v3_ll) + end + end + return f +end + @inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, normal_direction_ll::AbstractVector, equations::IdealGlmMhdEquations3D, @@ -487,6 +608,10 @@ This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl end """ + flux_nonconservative_powell_local_symmetric(u_ll, orientation::Integer, + equations::IdealGlmMhdEquations3D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) flux_nonconservative_powell_local_symmetric(u_ll, normal_direction_avg::AbstractVector, equations::IdealGlmMhdEquations3D, nonconservative_type::NonConservativeSymmetric, @@ -501,6 +626,67 @@ the non-conservative staggered "fluxes" for subcell limiting. See, e.g., This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl. """ +@inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdEquations3D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr + + if nonconservative_term == 1 + # Powell nonconservative term: (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0) + if orientation == 1 + B1_avg = (B1_ll + B1_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + B1_avg, + 0) + elseif orientation == 2 + B2_avg = (B2_ll + B2_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + B2_avg, + 0) + else # orientation == 3 + B3_avg = (B3_ll + B3_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + B3_avg, + B3_avg, + B3_avg, + B3_avg, + B3_avg, + B3_avg, + B3_avg, + 0) + end + else #nonconservative_term == 2 + # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2}) + psi_avg = (psi_ll + psi_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code + f = SVector(0, + 0, + 0, + 0, + psi_avg, + 0, + 0, + 0, + psi_avg) + end + + return f +end + @inline function (noncons_flux::FluxNonConservativePowellLocalSymmetric)(u_ll, u_rr, normal_direction_avg::AbstractVector, equations::IdealGlmMhdEquations3D, diff --git a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl index f07a852606e..50a27514f2c 100644 --- a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl @@ -433,4 +433,412 @@ end return nothing end + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**with non-conservative terms in "local * jump" form**). +# +# See also `flux_differencing_kernel!`. +# +# The calculation of the non-conservative staggered "fluxes" requires non-conservative +# terms that can be written as a product of local and jump contributions. +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, + u, mesh::P4estMesh{3}, + nonconservative_terms::True, equations, + volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM, + element, + cache) where { + F_CONS <: Function, + F_NONCONS <: + FluxNonConservative{NonConservativeJump()} + } + (; contravariant_vectors) = cache.elements + (; weights, derivative_split) = dg.basis + (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache + (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()] + + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()] + phi = phi_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + + # pull the contravariant vectors in each coordinate direction + Ja1_node = get_contravariant_vector(1, contravariant_vectors, + i, j, k, element) # x direction + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux + # computations. + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element) + # pull the contravariant vectors and compute the average + Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, + ii, j, k, element) + Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii) + + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde1 = volume_flux_cons(u_node, u_node_ii, Ja1_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1, + equations, dg, ii, j, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_avg, + equations, + NonConservativeJump(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[i, ii], + flux1_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5f0 * derivative_split[ii, i], + flux1_noncons, + equations, dg, noncons, ii, j, k) + end + end + end + + # FV-form flux `fhat` in x direction + fhat_temp[:, 1, :, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + # pull the local contravariant vector + Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, Ja1_node, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[i] * flux_temp[v, i, j, k] + fhat_temp[v, i + 1, j, k] = value + fhat1_L[v, i + 1, j, k] = value + fhat1_R[v, i + 1, j, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[i] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i + 1, j, k] = value + + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + + phi[v, noncons, i, j, k] * value + fhat1_R[v, i + 1, j, k] = fhat1_R[v, i + 1, j, k] + + phi[v, noncons, i + 1, j, k] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for k in eachnode(dg), j in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, 1, j, k, element) + Ja1_node_0 = get_contravariant_vector(1, contravariant_vectors, + 1, j, k, element) + + for i in 2:(nnodes(dg) - 1) + u_i = get_node_vars(u, equations, dg, i, j, k, element) + Ja1_node_i = get_contravariant_vector(1, contravariant_vectors, + i, j, k, element) + Ja1_avg = 0.5f0 * (Ja1_node_0 + Ja1_node_i) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_i, Ja1_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat1_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat1_L[v, i + 1, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, nnodes(dg), j, k, element) + Ja1_node_N = get_contravariant_vector(1, contravariant_vectors, + nnodes(dg), j, k, element) + Ja1_avg = 0.5f0 * (Ja1_node_0 + Ja1_node_N) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, Ja1_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing because Trixi multiplies all the non-cons terms with 0.5 + fhat1_R[v, nnodes(dg), j, k] -= phi[v, noncons, nnodes(dg), j, k] * + phi_jump[v] + end + end + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + + # pull the contravariant vectors in each coordinate direction + Ja2_node = get_contravariant_vector(2, contravariant_vectors, + i, j, k, element) + + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element) + # pull the contravariant vectors and compute the average + Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, + i, jj, k, element) + Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj) + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde2 = volume_flux_cons(u_node, u_node_jj, Ja2_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2, + equations, dg, i, jj, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_avg, + equations, + NonConservativeJump(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[j, jj], + flux2_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5f0 * derivative_split[jj, j], + flux2_noncons, + equations, dg, noncons, i, jj, k) + end + end + end + + # FV-form flux `fhat` in y direction + fhat_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + # pull the local contravariant vector + Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, Ja2_node, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[j] * flux_temp[v, i, j, k] + fhat_temp[v, i, j + 1, k] = value + fhat2_L[v, i, j + 1, k] = value + fhat2_R[v, i, j + 1, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[j] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j + 1, k] = value + + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + + phi[v, noncons, i, j, k] * value + fhat2_R[v, i, j + 1, k] = fhat2_R[v, i, j + 1, k] + + phi[v, noncons, i, j + 1, k] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for k in eachnode(dg), i in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, i, 1, k, element) + Ja2_node_0 = get_contravariant_vector(2, contravariant_vectors, + i, 1, k, element) + + for j in 2:(nnodes(dg) - 1) + u_j = get_node_vars(u, equations, dg, i, j, k, element) + Ja2_node_j = get_contravariant_vector(2, contravariant_vectors, + i, j, k, element) + Ja2_avg = 0.5f0 * (Ja2_node_0 + Ja2_node_j) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_j, Ja2_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat2_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat2_L[v, i, j + 1, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, i, nnodes(dg), k, element) + Ja2_node_N = get_contravariant_vector(2, contravariant_vectors, + i, nnodes(dg), k, element) + Ja2_avg = 0.5f0 * (Ja2_node_0 + Ja2_node_N) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, Ja2_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing cause Trixi multiplies all the non-cons terms with 0.5 + fhat2_R[v, i, nnodes(dg), k] -= phi[v, noncons, i, nnodes(dg), k] * + phi_jump[v] + end + end + end + + # Split form volume flux in orientation 3: z direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + + # pull the contravariant vectors in each coordinate direction + Ja3_node = get_contravariant_vector(3, contravariant_vectors, + i, j, k, element) + + for kk in (k + 1):nnodes(dg) + u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element) + # pull the contravariant vectors and compute the average + Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, + i, j, kk, element) + Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk) + # compute the contravariant sharp flux in the direction of the averaged contravariant vector + fluxtilde3 = volume_flux_cons(u_node, u_node_kk, Ja3_avg, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[k, kk], fluxtilde3, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[kk, k], fluxtilde3, + equations, dg, i, j, kk) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux3_noncons = volume_flux_noncons(u_node, u_node_kk, Ja3_avg, + equations, + NonConservativeJump(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[k, kk], + flux3_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5f0 * derivative_split[kk, k], + flux3_noncons, + equations, dg, noncons, i, j, kk) + end + end + end + + # FV-form flux `fhat` in z direction + fhat_temp[:, :, :, 1] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, :, 1] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + # pull the local contravariant vector + Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, Ja3_node, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[k] * flux_temp[v, i, j, k] + fhat_temp[v, i, j, k + 1] = value + fhat3_L[v, i, j, k + 1] = value + fhat3_R[v, i, j, k + 1] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[k] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j, k + 1] = value + + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + + phi[v, noncons, i, j, k] * value + fhat3_R[v, i, j, k + 1] = fhat3_R[v, i, j, k + 1] + + phi[v, noncons, i, j, k + 1] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for j in eachnode(dg), i in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, i, j, 1, element) + Ja3_node_0 = get_contravariant_vector(3, contravariant_vectors, + i, j, 1, element) + + for k in 2:(nnodes(dg) - 1) + u_k = get_node_vars(u, equations, dg, i, j, k, element) + Ja3_node_k = get_contravariant_vector(3, contravariant_vectors, + i, j, k, element) + Ja3_avg = 0.5f0 * (Ja3_node_0 + Ja3_node_k) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_k, Ja3_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat3_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat3_L[v, i, j, k + 1] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, i, j, nnodes(dg), element) + Ja3_node_N = get_contravariant_vector(3, contravariant_vectors, + i, j, nnodes(dg), element) + Ja3_avg = 0.5f0 * (Ja3_node_0 + Ja3_node_N) + + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, Ja3_avg, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing cause Trixi multiplies all the non-cons terms with 0.5 + fhat3_R[v, i, j, nnodes(dg)] -= phi[v, noncons, i, j, nnodes(dg)] * + phi_jump[v] + end + end + end + + return nothing +end end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl index 5d3d9d86069..eb4f09b4004 100644 --- a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl @@ -217,6 +217,602 @@ end return nothing end +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**with non-conservative terms in "local * symmetric" form**). +# +# See also `flux_differencing_kernel!`. +# +# The calculation of the non-conservative staggered "fluxes" requires non-conservative +# terms that can be written as a product of local and a symmetric contributions. See, e.g., +# +# - Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts +# Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +# +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, + u, mesh::TreeMesh{2}, + have_nonconservative_terms::True, equations, + volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM, + element, + cache) where { + F_CONS <: Function, + F_NONCONS <: + FluxNonConservative{NonConservativeSymmetric()} + } + @unpack weights, derivative_split = dg.basis + @unpack flux_temp_threaded, flux_nonconservative_temp_threaded = cache + @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()] + + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()] + phi = phi_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of `volume_flux_cons` and `volume_flux_noncons` to save half of the possible two-point flux + # computations. + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element) + flux1 = volume_flux_cons(u_node, u_node_ii, 1, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], flux1, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, + equations, dg, ii, j, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[i, ii], + flux1_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[ii, i], + flux1_noncons, + equations, dg, noncons, ii, j, k) + end + end + end + + # FV-form flux `fhat` in x direction + fhat_temp[:, 1, :, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 1, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[i] * flux_temp[v, i, j, k] + fhat_temp[v, i + 1, j, k] = value + fhat1_L[v, i + 1, j, k] = value + fhat1_R[v, i + 1, j, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[i] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i + 1, j, k] = value + + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + + phi[v, noncons, i, j, k] * value + fhat1_R[v, i + 1, j, k] = fhat1_R[v, i + 1, j, k] + + phi[v, noncons, i + 1, j, k] * value + end + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element) + flux2 = volume_flux_cons(u_node, u_node_jj, 2, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], flux2, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, + equations, dg, i, jj, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[j, jj], + flux2_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[jj, j], + flux2_noncons, + equations, dg, noncons, i, jj, k) + end + end + end + + # FV-form flux `fhat` in y direction + fhat_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 2, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[j] * flux_temp[v, i, j, k] + fhat_temp[v, i, j + 1, k] = value + fhat2_L[v, i, j + 1, k] = value + fhat2_R[v, i, j + 1, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[j] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j + 1, k] = value + + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + + phi[v, noncons, i, j, k] * value + fhat2_R[v, i, j + 1, k] = fhat2_R[v, i, j + 1, k] + + phi[v, noncons, i, j + 1, k] * value + end + end + + # Split form volume flux in orientation 3: z direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + for kk in (k + 1):nnodes(dg) + u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element) + flux3 = volume_flux_cons(u_node, u_node_kk, 3, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[k, kk], flux3, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[kk, k], flux3, + equations, dg, i, j, kk) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux3_noncons = volume_flux_noncons(u_node, u_node_kk, 3, equations, + NonConservativeSymmetric(), noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[k, kk], + flux3_noncons, + equations, dg, noncons, i, j, k) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[kk, k], + flux3_noncons, + equations, dg, noncons, i, j, kk) + end + end + end + + # FV-form flux `fhat` in z direction + fhat_temp[:, :, :, 1] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, :, 1] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 3, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[k] * flux_temp[v, i, j, k] + fhat_temp[v, i, j, k + 1] = value + fhat3_L[v, i, j, k + 1] = value + fhat3_R[v, i, j, k + 1] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[k] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j, k + 1] = value + + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + + phi[v, noncons, i, j, k] * value + fhat3_R[v, i, j, k + 1] = fhat3_R[v, i, j, k + 1] + + phi[v, noncons, i, j, k + 1] * value + end + end + + return nothing +end + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**with non-conservative terms in "local * jump" form**). +# +# See also `flux_differencing_kernel!`. +# +# The calculation of the non-conservative staggered "fluxes" requires non-conservative +# terms that can be written as a product of local and jump contributions. +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, + u, mesh::TreeMesh{2}, + nonconservative_terms::True, equations, + volume_flux::Tuple{F_CONS, F_NONCONS}, dg::DGSEM, + element, + cache) where { + F_CONS <: Function, + F_NONCONS <: + FluxNonConservative{NonConservativeJump()} + } + @unpack weights, derivative_split = dg.basis + @unpack flux_temp_threaded, flux_nonconservative_temp_threaded = cache + @unpack fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded = cache + + volume_flux_cons, volume_flux_noncons = volume_flux + + flux_temp = flux_temp_threaded[Threads.threadid()] + flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()] + + fhat_temp = fhat_temp_threaded[Threads.threadid()] + fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()] + phi = phi_threaded[Threads.threadid()] + + # The FV-form fluxes are calculated in a recursive manner, i.e.: + # fhat_(0,1) = w_0 * FVol_0, + # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j, for j=1,...,N-1, + # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i). + + # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated + # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing` + # and saved in in `flux_temp`. + + # Split form volume flux in orientation 1: x direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + + # All diagonal entries of `derivative_split` are zero. Thus, we can skip + # the computation of the diagonal terms. In addition, we use the symmetry + # of `volume_flux_cons` and skew-symmetry of `volume_flux_noncons` to save half of the possible two-point flux + # computations. + for ii in (i + 1):nnodes(dg) + u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element) + flux1 = volume_flux_cons(u_node, u_node_ii, 1, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], flux1, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1, + equations, dg, ii, j, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations, + NonConservativeJump(), + noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5f0 * derivative_split[i, ii], + flux1_noncons, + equations, dg, noncons, i, j, k) + # Note the sign flip due to skew-symmetry when argument order is swapped + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5f0 * derivative_split[ii, i], + flux1_noncons, + equations, dg, noncons, ii, j, k) + end + end + end + + # FV-form flux `fhat` in x direction + fhat_temp[:, 1, :, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, 1, :, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 1, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[i] * flux_temp[v, i, j, k] + fhat_temp[v, i + 1, j, k] = value + fhat1_L[v, i + 1, j, k] = value + fhat1_R[v, i + 1, j, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[i] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i + 1, j, k] = value + + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + + phi[v, noncons, i, j, k] * value + fhat1_R[v, i + 1, j, k] = fhat1_R[v, i + 1, j, k] + + phi[v, noncons, i + 1, j, k] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for k in eachnode(dg), j in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, 1, j, k, element) + for i in 2:(nnodes(dg) - 1) + u_i = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_i, 1, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat1_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat1_L[v, i + 1, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, nnodes(dg), j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, 1, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing because Trixi multiplies all the non-cons terms with 0.5 + fhat1_R[v, nnodes(dg), j, k] -= phi[v, noncons, nnodes(dg), j, k] * + phi_jump[v] + end + end + end + + ######## + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + for jj in (j + 1):nnodes(dg) + u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element) + flux2 = volume_flux_cons(u_node, u_node_jj, 2, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], flux2, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2, + equations, dg, i, jj, k) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations, + NonConservativeJump(), + noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[j, jj], + flux2_noncons, + equations, dg, noncons, i, j, k) + # Note the sign flip due to skew-symmetry when argument order is swapped + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5 * derivative_split[jj, j], + flux2_noncons, + equations, dg, noncons, i, jj, k) + end + end + end + + # FV-form flux `fhat` in y direction + fhat_temp[:, :, 1, :] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, 1, :] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 2, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[j] * flux_temp[v, i, j, k] + fhat_temp[v, i, j + 1, k] = value + fhat2_L[v, i, j + 1, k] = value + fhat2_R[v, i, j + 1, k] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[j] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j + 1, k] = value + + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + + phi[v, noncons, i, j, k] * value + fhat2_R[v, i, j + 1, k] = fhat2_R[v, i, j + 1, k] + + phi[v, noncons, i, j + 1, k] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for k in eachnode(dg), i in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, i, 1, k, element) + for j in 2:(nnodes(dg) - 1) + u_j = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_j, 2, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat2_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat2_L[v, i, j + 1, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, i, nnodes(dg), k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, 2, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing cause Trixi multiplies all the non-cons terms with 0.5 + fhat2_R[v, i, nnodes(dg), k] -= phi[v, noncons, i, nnodes(dg), k] * + phi_jump[v] + end + end + end + + ######## + + # Split form volume flux in orientation 3: z direction + flux_temp .= zero(eltype(flux_temp)) + flux_noncons_temp .= zero(eltype(flux_noncons_temp)) + + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_node = get_node_vars(u, equations, dg, i, j, k, element) + for kk in (k + 1):nnodes(dg) + u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element) + flux3 = volume_flux_cons(u_node, u_node_kk, 3, equations) + multiply_add_to_node_vars!(flux_temp, derivative_split[k, kk], flux3, + equations, dg, i, j, k) + multiply_add_to_node_vars!(flux_temp, derivative_split[kk, k], flux3, + equations, dg, i, j, kk) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + # We multiply by 0.5 because that is done in other parts of Trixi + flux3_noncons = volume_flux_noncons(u_node, u_node_kk, 3, equations, + NonConservativeJump(), + noncons) + multiply_add_to_node_vars!(flux_noncons_temp, + 0.5 * derivative_split[k, kk], + flux3_noncons, + equations, dg, noncons, i, j, k) + # Note the sign flip due to skew-symmetry when argument order is swapped + multiply_add_to_node_vars!(flux_noncons_temp, + -0.5 * derivative_split[kk, k], + flux3_noncons, + equations, dg, noncons, i, j, kk) + end + end + end + + # FV-form flux `fhat` in y direction + fhat_temp[:, :, :, 1] .= zero(eltype(fhat1_L)) + fhat_noncons_temp[:, :, :, :, 1] .= zero(eltype(fhat1_L)) + + # Compute local contribution to non-conservative flux + for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg) + u_local = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + set_node_vars!(phi, + volume_flux_noncons(u_local, 3, equations, + NonConservativeLocal(), noncons), + equations, dg, noncons, i, j, k) + end + end + + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg) + # Conservative part + for v in eachvariable(equations) + value = fhat_temp[v, i, j, k] + weights[k] * flux_temp[v, i, j, k] + fhat_temp[v, i, j, k + 1] = value + fhat3_L[v, i, j, k + 1] = value + fhat3_R[v, i, j, k + 1] = value + end + # Nonconservative part + for noncons in 1:n_nonconservative_terms(volume_flux_noncons), + v in eachvariable(equations) + + value = fhat_noncons_temp[v, noncons, i, j, k] + + weights[k] * flux_noncons_temp[v, noncons, i, j, k] + fhat_noncons_temp[v, noncons, i, j, k + 1] = value + + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + + phi[v, noncons, i, j, k] * value + fhat3_R[v, i, j, k + 1] = fhat3_R[v, i, j, k + 1] + + phi[v, noncons, i, j, k + 1] * value + end + end + + # Apply correction term to the flux-differencing formula for nonconservative local * jump fluxes. + for j in eachnode(dg), i in eachnode(dg) + u_0 = get_node_vars(u, equations, dg, i, j, 1, element) + for k in 2:(nnodes(dg) - 1) + u_k = get_node_vars(u, equations, dg, i, j, k, element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_k, 2, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing on each term because Trixi multiplies all the non-cons terms with 0.5 + fhat3_R[v, i, j, k] -= phi[v, noncons, i, j, k] * phi_jump[v] + fhat3_L[v, i, j, k + 1] -= phi[v, noncons, i, j, k] * phi_jump[v] + end + end + end + u_N = get_node_vars(u, equations, dg, i, j, nnodes(dg), element) + for noncons in 1:n_nonconservative_terms(volume_flux_noncons) + phi_jump = volume_flux_noncons(u_0, u_N, 3, equations, + NonConservativeJump(), noncons) + + for v in eachvariable(equations) + # The factor of 2 is missing cause Trixi multiplies all the non-cons terms with 0.5 + fhat3_R[v, i, j, nnodes(dg)] -= phi[v, noncons, i, j, nnodes(dg)] * + phi_jump[v] + end + end + end + return nothing +end + # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems. @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, diff --git a/src/solvers/dgsem_tree/subcell_limiters_3d.jl b/src/solvers/dgsem_tree/subcell_limiters_3d.jl index 5c8d6c06277..7a5081afde5 100644 --- a/src/solvers/dgsem_tree/subcell_limiters_3d.jl +++ b/src/solvers/dgsem_tree/subcell_limiters_3d.jl @@ -68,7 +68,7 @@ end u, t, semi, mesh::TreeMesh3D, equations) _, _, dg, cache = mesh_equations_solver_cache(semi) - + (; boundary_conditions) = semi # Calc bounds at interfaces and periodic boundaries for interface in eachinterface(dg, cache) # Get neighboring element ids @@ -109,6 +109,50 @@ end end end + # Calc bounds at physical boundaries + for boundary in eachboundary(dg, cache) + element = cache.boundaries.neighbor_ids[boundary] + orientation = cache.boundaries.orientations[boundary] + neighbor_side = cache.boundaries.neighbor_sides[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Define node indices and boundary index based on the orientation and neighbor_side + if neighbor_side == 2 # Element is on the right, boundary on the left + if orientation == 1 # boundary in x-direction + index = (1, i, j) + boundary_index = 1 + elseif orientation == 2 # boundary in y-direction + index = (i, 1, j) + boundary_index = 3 + else # orientation == 3 # boundary in z-direction + index = (i, j, 1) + boundary_index = 5 + end + else # Element is on the left, boundary on the right + if orientation == 1 # boundary in x-direction + index = (nnodes(dg), i, j) + boundary_index = 2 + elseif orientation == 2 # boundary in y-direction + index = (i, nnodes(dg), j) + boundary_index = 4 + else # orientation == 3 # boundary in z-direction + index = (i, j, nnodes(dg)) + boundary_index = 6 + end + end + u_inner = get_node_vars(u, equations, dg, index..., element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[boundary_index], + orientation, boundary_index, + mesh, equations, dg, cache, + index..., element) + var_outer = u_outer[variable] + + var_min[index..., element] = min(var_min[index..., element], var_outer) + var_max[index..., element] = max(var_max[index..., element], var_outer) + end + end + return nothing end @@ -171,6 +215,7 @@ end @inline function calc_bounds_onesided_interface!(var_minmax, min_or_max, variable, u, t, semi, mesh::TreeMesh{3}) _, equations, dg, cache = mesh_equations_solver_cache(semi) + (; boundary_conditions) = semi # Calc bounds at interfaces and periodic boundaries for interface in eachinterface(dg, cache) @@ -211,6 +256,50 @@ end end end + # Calc bounds at physical boundaries + for boundary in eachboundary(dg, cache) + element = cache.boundaries.neighbor_ids[boundary] + orientation = cache.boundaries.orientations[boundary] + neighbor_side = cache.boundaries.neighbor_sides[boundary] + + for j in eachnode(dg), i in eachnode(dg) + # Define node indices and boundary index based on the orientation and neighbor_side + if neighbor_side == 2 # Element is on the right, boundary on the left + if orientation == 1 # boundary in x-direction + index = (1, i, j) + boundary_index = 1 + elseif orientation == 2 # boundary in y-direction + index = (i, 1, j) + boundary_index = 3 + else # orientation == 3 # boundary in z-direction + index = (i, j, 1) + boundary_index = 5 + end + else # Element is on the left, boundary on the right + if orientation == 1 # boundary in x-direction + index = (nnodes(dg), i, j) + boundary_index = 2 + elseif orientation == 2 # boundary in y-direction + index = (i, nnodes(dg), j) + boundary_index = 4 + else # orientation == 3 # boundary in z-direction + index = (i, j, nnodes(dg)) + boundary_index = 6 + end + end + u_inner = get_node_vars(u, equations, dg, index..., element) + u_outer = get_boundary_outer_state(u_inner, t, + boundary_conditions[boundary_index], + orientation, boundary_index, + mesh, equations, dg, cache, + index..., element) + var_outer = variable(u_outer, equations) + + var_minmax[index..., element] = min_or_max(var_minmax[index..., element], + var_outer) + end + end + return nothing end diff --git a/test/test_p4est_3d.jl b/test/test_p4est_3d.jl index 12b2fa92e91..bf05e823f59 100644 --- a/test/test_p4est_3d.jl +++ b/test/test_p4est_3d.jl @@ -622,6 +622,64 @@ end @test_allocations(Trixi.rhs!, semi, sol, 1000) end +@trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl (positivity bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl"), + l2=[ + 0.0010625880666563514, + 0.0011220178360278588, + 0.001175570902488149, + 0.0013866847769005606, + 0.002942603300074813 + ], + linf=[ + 0.019666290518305374, + 0.019355047676718584, + 0.02588920954597862, + 0.022936612741258244, + 0.04983574326863449 + ], + tspan=(0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test_allocations(Trixi.rhs!, semi, sol, 15_000) +end + +@trixi_testset "elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl (local bounds)" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_source_terms_nonperiodic_hohqmesh_sc_subcell.jl"), + local_twosided_variables_cons=["rho"], + local_onesided_variables_nonlinear=[(entropy_guermond_etal, + min)], + max_iterations_newton=30, + l2=[ + 0.01573068932681815, + 0.01591679096422271, + 0.015795998889146637, + 0.01746834678809081, + 0.06766354772185222 + ], + linf=[ + 0.06511438601036024, + 0.09961666854491202, + 0.10474759248351417, + 0.07545726460168245, + 0.26949807234431944 + ], + tspan=(0.0, 0.3)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + # Larger values for allowed allocations due to usage of custom + # integrator which are not *recorded* for the methods from + # OrdinaryDiffEq.jl + # Corresponding issue: https://github.com/trixi-framework/Trixi.jl/issues/1877 + @test_allocations(Trixi.rhs!, semi, sol, 15_000) +end + @trixi_testset "elixir_mhd_alfven_wave_er.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_er.jl"), diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl index 17866e74f47..ab0d0925dcd 100644 --- a/test/test_tree_3d_euler.jl +++ b/test/test_tree_3d_euler.jl @@ -30,7 +30,7 @@ EXAMPLES_DIR = joinpath(examples_dir(), "tree_3d_dgsem") # (e.g., from type instabilities) @test_allocations(Trixi.rhs!, semi, sol, 1000) # Extra test to make sure the "TimeSeriesCallback" made correct data. - # Extracts data at all points from the first step of the time series and compares it to the + # Extracts data at all points from the first step of the time series and compares it to the # exact solution and an interpolated reference solution point_data = [getindex(time_series.affect!.point_data[i], 1:5) for i in 1:3] exact_data = [initial_condition_convergence_test(time_series.affect!.point_coordinates[:,