From 1a686d21c4d563e56a323b7e1ab22cd7fe67a9be Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 4 Mar 2026 09:34:27 +0100 Subject: [PATCH 1/3] Add initial support for 3d TreeMesh --- ...lixir_euler_sedov_blast_wave_sc_subcell.jl | 98 ++++++++++++++++ .../subcell_limiter_idp_correction_3d.jl | 2 +- .../dgsem_p4est/dg_3d_subcell_limiters.jl | 8 +- src/solvers/dgsem_tree/dg.jl | 1 + .../dgsem_tree/dg_3d_subcell_limiters.jl | 109 ++++++++++++++++++ test/test_tree_3d_euler.jl | 27 +++++ 6 files changed, 240 insertions(+), 5 deletions(-) create mode 100644 examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl create mode 100644 src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl 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 new file mode 100644 index 00000000000..852fd5ffdbe --- /dev/null +++ b/examples/tree_3d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl @@ -0,0 +1,98 @@ +using Trixi + +############################################################################### +# semidiscretization of the compressible Euler equations + +equations = CompressibleEulerEquations3D(1.4) + +""" + initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D) + +The Sedov blast wave setup based on example 35.1.4 from Flash +- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf +with smaller strength of the initial discontinuity. +""" +function initial_condition_sedov_blast_wave(x, t, + equations::CompressibleEulerEquations3D) + # Set up polar coordinates + inicenter = SVector(0.0, 0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + z_norm = x[3] - inicenter[3] + r = sqrt(x_norm^2 + y_norm^2 + z_norm^2) + + # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000 + r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6) + E = 1.0 + p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2) + p0_outer = 1.0e-5 + + # Calculate primitive variables + rho = 1.0 + v1 = 0.0 + v2 = 0.0 + v3 = 0.0 + p = r > r0 ? p0_outer : p0_inner + + return prim2cons(SVector(rho, v1, v2, v3, p), equations) +end +initial_condition = initial_condition_sedov_blast_wave + +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], + max_iterations_newton = 20) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + +coordinates_min = (-1.0, -1.0, -1.0) +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) + +# create the semi discretization object +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 3.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 = 10, + save_initial_solution = true, + save_final_solution = true, + 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/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl b/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl index f96a491487c..8390c47821a 100644 --- a/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl +++ b/src/callbacks_stage/subcell_limiter_idp_correction_3d.jl @@ -6,7 +6,7 @@ #! format: noindent function perform_idp_correction!(u, dt, - mesh::P4estMesh{3}, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, dg, cache) @unpack inverse_weights = dg.basis # Plays role of inverse DG-subcell sizes @unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes diff --git a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl index 00462d05786..aee9efb52a6 100644 --- a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl @@ -5,7 +5,7 @@ @muladd begin #! format: noindent -function create_cache(mesh::P4estMesh{3}, +function create_cache(mesh::Union{TreeMesh{3}, P4estMesh{3}}, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DG, cache_containers, uEltype) cache = create_cache(mesh, equations, @@ -61,7 +61,7 @@ end # Subcell limiting currently only implemented for certain mesh types @inline function volume_integral_kernel!(du, u, element, - mesh::P4estMesh{3}, + mesh::Union{TreeMesh{3}, P4estMesh{3}}, nonconservative_terms, equations, volume_integral::VolumeIntegralSubcellLimiting, dg::DGSEM, cache) @@ -551,7 +551,7 @@ end fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, - u, mesh::P4estMesh{3}, + u, mesh::Union{TreeMesh{3}, P4estMesh{3}}, nonconservative_terms::False, equations, limiter::SubcellLimiterIDP, dg, element, cache) @unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes @@ -602,7 +602,7 @@ end fstar1_L, fstar1_R, fstar2_L, fstar2_R, fstar3_L, fstar3_R, - u, mesh::P4estMesh{3}, + u, mesh::Union{TreeMesh{3}, P4estMesh{3}}, nonconservative_terms::True, equations, limiter::SubcellLimiterIDP, dg, element, cache) @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes diff --git a/src/solvers/dgsem_tree/dg.jl b/src/solvers/dgsem_tree/dg.jl index 54d71af2b66..463f01ad734 100644 --- a/src/solvers/dgsem_tree/dg.jl +++ b/src/solvers/dgsem_tree/dg.jl @@ -63,4 +63,5 @@ include("dg_3d_compressible_euler.jl") include("subcell_limiters.jl") include("subcell_limiters_2d.jl") include("dg_2d_subcell_limiters.jl") +include("dg_3d_subcell_limiters.jl") end # @muladd diff --git a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl new file mode 100644 index 00000000000..6229054c569 --- /dev/null +++ b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl @@ -0,0 +1,109 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element +# (**without non-conservative terms**). +# +# See also `flux_differencing_kernel!`. +@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, + u, mesh::TreeMesh{3}, + have_nonconservative_terms::False, equations, + volume_flux, dg::DGSEM, element, cache) + @unpack weights, derivative_split = dg.basis + @unpack flux_temp_threaded = cache + + flux_temp = flux_temp_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)) + + 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 the `volume_flux` 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(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) + end + end + + # FV-form flux `fhat` in x direction + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1), + v in eachvariable(equations) + + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] + + weights[i] * flux_temp[v, i, j, k] + fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + end + + # Split form volume flux in orientation 2: y direction + flux_temp .= zero(eltype(flux_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(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) + end + end + + # FV-form flux `fhat` in y direction + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg), + v in eachvariable(equations) + + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] + + weights[j] * flux_temp[v, i, j, k] + fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + end + + # Split form volume flux in orientation 3: z direction + flux_temp .= zero(eltype(flux_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(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) + end + end + + # FV-form flux `fhat` in z direction + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg), + v in eachvariable(equations) + + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] + + weights[k] * flux_temp[v, i, j, k] + fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + end + + return nothing +end +end # @muladd diff --git a/test/test_tree_3d_euler.jl b/test/test_tree_3d_euler.jl index 1bbcc7ea558..b3ea6a05e4d 100644 --- a/test/test_tree_3d_euler.jl +++ b/test/test_tree_3d_euler.jl @@ -498,6 +498,33 @@ end # (e.g., from type instabilities) @test_allocations(Trixi.rhs!, semi, sol, 1000) end + +@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_sedov_blast_wave_sc_subcell.jl"), + l2=[ + 0.30386298122471, + 0.08632041574897817, + 0.086320415748978, + 0.08632041574897818, + 0.36209285851713263 + ], + linf=[ + 2.7876219089976964, + 1.3892133662695192, + 1.3892133662695196, + 1.3892133662695192, + 4.854504215641112 + ], + tspan=(0.0, 0.5)) + # 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 end end # module From 7ce64a2e877ca08937b78585669868f03005fc9c Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 4 Mar 2026 09:55:48 +0100 Subject: [PATCH 2/3] Update NEWS.md --- NEWS.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/NEWS.md b/NEWS.md index edddc607a42..f02da93dbf4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -27,9 +27,10 @@ This enables in particular adaptive mesh refinement for that solver-mesh combina Instead, it is possible to reconstruct in custom variables, if functions `cons2recon` and `recon2cons` are provided to `VolumeIntegralPureLGLFiniteVolumeO2` and `VolumeIntegralShockCapturingRRG`([#2817]). - Add Legendre-Gauss basis for DGSEM and implement solver (`VolumeIntegralWeakForm` and `SurfaceIntegralWeakForm` only) support for conforming 1D & 2D `TreeMesh`es ([#1965]). -- Extended 3D support for subcell limiting with `P4estMesh` was added ([#2763]). +- Extended 3D support for subcell limiting was added ([#2763], [#2846]). In the new version, local (minimum and maximum) limiting for conservative variables (using the - keyword `local_twosided_variables_cons` in `SubcellLimiterIDP()`) is supported. + keyword `local_twosided_variables_cons` in `SubcellLimiterIDP()`) with `P4estMesh` is supported. + Moreover, initial support for `TreeMesh{3}` including positivity limiting on periodic meshes was added. ## Changes when updating to v0.15 from v0.14.x From f59e8d128be92c60614c1299cc30d5018d51b3d3 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 5 Mar 2026 08:54:55 +0100 Subject: [PATCH 3/3] Split `for` loops over nodes and variables --- .../dgsem_p4est/dg_3d_subcell_limiters.jl | 36 +++++++++---------- .../dg_2d_subcell_limiters.jl | 16 +++++---- .../dgsem_tree/dg_2d_subcell_limiters.jl | 16 +++++---- .../dgsem_tree/dg_3d_subcell_limiters.jl | 36 +++++++++---------- 4 files changed, 56 insertions(+), 48 deletions(-) diff --git a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl index aee9efb52a6..189e1e85b78 100644 --- a/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_p4est/dg_3d_subcell_limiters.jl @@ -171,12 +171,12 @@ end end # FV-form flux `fhat` in x direction - for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1), - v in eachvariable(equations) - - fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] + - weights[i] * flux_temp[v, i, j, k] - fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1) + for v in eachvariable(equations) + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] + + weights[i] * flux_temp[v, i, j, k] + fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + end end # Split form volume flux in orientation 2: y direction @@ -206,12 +206,12 @@ end end # FV-form flux `fhat` in y direction - for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg), - v in eachvariable(equations) - - fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] + - weights[j] * flux_temp[v, i, j, k] - fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg) + for v in eachvariable(equations) + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] + + weights[j] * flux_temp[v, i, j, k] + fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + end end # Split form volume flux in orientation 3: z direction @@ -241,12 +241,12 @@ end end # FV-form flux `fhat` in z direction - for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg), - v in eachvariable(equations) - - fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] + - weights[k] * flux_temp[v, i, j, k] - fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] + + weights[k] * flux_temp[v, i, j, k] + fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + end end return nothing diff --git a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl index f29de2d912d..bddf6cbf618 100644 --- a/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl @@ -60,9 +60,11 @@ end # FV-form flux `fhat` in x direction - for j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations) - fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] - fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] + for j in eachnode(dg), i in 1:(nnodes(dg) - 1) + for v in eachvariable(equations) + fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] + end end # Split form volume flux in orientation 2: y direction @@ -91,9 +93,11 @@ end # FV-form flux `fhat` in y direction - for j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations) - fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] - fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + for j in 1:(nnodes(dg) - 1), i in eachnode(dg) + for v in eachvariable(equations) + fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + end end return nothing diff --git a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl index 292a7952a40..de87fadae85 100644 --- a/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl @@ -151,9 +151,11 @@ end end # FV-form flux `fhat` in x direction - for j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations) - fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] - fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] + for j in eachnode(dg), i in 1:(nnodes(dg) - 1) + for v in eachvariable(equations) + fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j] + fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j] + end end # Split form volume flux in orientation 2: y direction @@ -172,9 +174,11 @@ end end # FV-form flux `fhat` in y direction - for j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations) - fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] - fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + for j in 1:(nnodes(dg) - 1), i in eachnode(dg) + for v in eachvariable(equations) + fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j] + fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1] + end end return nothing diff --git a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl index 6229054c569..9154bed5c58 100644 --- a/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl +++ b/src/solvers/dgsem_tree/dg_3d_subcell_limiters.jl @@ -48,12 +48,12 @@ end # FV-form flux `fhat` in x direction - for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1), - v in eachvariable(equations) - - fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] + - weights[i] * flux_temp[v, i, j, k] - fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1) + for v in eachvariable(equations) + fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] + + weights[i] * flux_temp[v, i, j, k] + fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k] + end end # Split form volume flux in orientation 2: y direction @@ -72,12 +72,12 @@ end # FV-form flux `fhat` in y direction - for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg), - v in eachvariable(equations) - - fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] + - weights[j] * flux_temp[v, i, j, k] - fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg) + for v in eachvariable(equations) + fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] + + weights[j] * flux_temp[v, i, j, k] + fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k] + end end # Split form volume flux in orientation 3: z direction @@ -96,12 +96,12 @@ end # FV-form flux `fhat` in z direction - for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg), - v in eachvariable(equations) - - fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] + - weights[k] * flux_temp[v, i, j, k] - fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg) + for v in eachvariable(equations) + fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] + + weights[k] * flux_temp[v, i, j, k] + fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1] + end end return nothing