Add solution_container and auxiliary_container for covariant DGMulti solver to match upstream refactor

src/callbacks_step/analysis_covariant.jl

@@ -26,8 +26,8 @@ function Trixi.integrate(func::Func, u,
                          equations::AbstractCovariantEquations,
                          dg::DGMulti, cache; normalize = true) where {Func}
     rd = dg.basis
-    md = mesh.md
-    @unpack u_values, aux_quad_values = cache
+    (; u_values) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
 
     # interpolate u to quadrature points
     Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u)
@@ -108,7 +108,8 @@ function Trixi.analyze(::typeof(Trixi.entropy_timederivative), du, u, t,
                        dg::DGMulti, cache)
     rd = dg.basis
     md = mesh.md
-    @unpack u_values, aux_quad_values = cache
+    (; u_values) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
 
     # interpolate u, du to quadrature points
     du_values = similar(u_values)
@@ -186,7 +187,8 @@ function Trixi.calc_error_norms(func, u, t, analyzer,
                                 cache_analysis) where {NDIMS, NDIMS_AMBIENT}
     rd = dg.basis
     md = mesh.md
-    @unpack u_values, aux_quad_values = cache
+    (; u_values) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
 
     # interpolate u to quadrature points
     Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u)

src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl

@@ -214,7 +214,7 @@ end
 @inline function calc_vorticity_node(u, equations::AbstractEquations,
                                      dg::DGMulti, cache, i, element)
     rd = dg.basis
-    (; aux_values) = cache
+    (; aux_values) = cache.auxiliary_container
     Dr, Ds = rd.Drst
 
     dv2dxi1 = dv1dxi2 = zero(eltype(u))

src/callbacks_step/save_solution_covariant.jl

@@ -36,7 +36,7 @@ end
 # Convert to another set of variables using a solution_variables function
 function convert_variables(u, solution_variables, mesh::DGMultiMesh,
                            equations::AbstractCovariantEquations, dg, cache)
-    (; aux_values) = cache
+    (; aux_values) = cache.auxiliary_container
     # Extract the number of solution variables to be output 
     # (may be different than the number of conservative variables) 
     n_vars = length(Trixi.varnames(solution_variables, equations))
@@ -76,7 +76,7 @@ end
 @inline function cons2prim_and_vorticity(u, mesh::DGMultiMesh,
                                          equations::AbstractCovariantEquations{2},
                                          dg::DGMulti, cache, i)
-    (; aux_values) = cache
+    (; aux_values) = cache.auxiliary_container
     u_node = u[:, i]
     aux_node = aux_values[i]
     element = (i - 1) ÷ nnodes(dg) + 1

src/callbacks_step/stepsize_dg2d.jl

@@ -74,7 +74,7 @@ function Trixi.max_dt(u, t, mesh::DGMultiMesh,
                       constant_speed::False,
                       equations::AbstractCovariantEquations{NDIMS},
                       dg::DGMulti, cache) where {NDIMS}
-    (; aux_values) = cache
+    (; aux_values) = cache.auxiliary_container
 
     dt_min = Inf
     for e in eachelement(mesh, dg, cache)

src/solvers/dgmulti/dg.jl

@@ -15,13 +15,21 @@ function Trixi.create_cache(mesh::DGMultiMesh{NDIMS}, equations::AbstractCovaria
     nvars = nvariables(equations)
     naux = n_aux_node_vars(equations)
 
-    # storage for volume quadrature values, face quadrature values, flux values
+    # We are duplicating the contents of solution_container in the top-level cache, but 
+    # note that no actual data is being copied here, just references to the same arrays.
     u_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xq), dg)
     u_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg)
     flux_face_values = Trixi.allocate_nested_array(uEltype, nvars, size(md.xf), dg)
+    local_values_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
+                             for _ in 1:Threads.maxthreadid()]
+    solution_container = (; u_values, u_face_values, flux_face_values,
+                          local_values_threaded)
+
+    # To parallel the solution container, we create an auxiliary container.
     aux_values = Trixi.allocate_nested_array(uEltype, naux, size(md.x), dg)
     aux_quad_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xq), dg)
     aux_face_values = Trixi.allocate_nested_array(uEltype, naux, size(md.xf), dg)
+    auxiliary_container = (; aux_values, aux_quad_values, aux_face_values)
 
     if typeof(rd.approximation_type) <:
        Union{SBP, Trixi.AbstractNonperiodicDerivativeOperator}
@@ -30,10 +38,6 @@ function Trixi.create_cache(mesh::DGMultiMesh{NDIMS}, equations::AbstractCovaria
         lift_scalings = nothing
     end
 
-    # local storage for volume integral and source computations
-    local_values_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
-
     # For curved meshes, we interpolate geometric terms from nodal points to quadrature points.
     # For affine meshes, we just access one element of this interpolated data.
     dxidxhatj = map(x -> rd.Vq * x, md.rstxyzJ)
@@ -49,14 +53,16 @@ function Trixi.create_cache(mesh::DGMultiMesh{NDIMS}, equations::AbstractCovaria
 
     # for scaling by curved geometric terms (not used by affine DGMultiMesh)
     flux_threaded = [[Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                      for _ in 1:NDIMS] for _ in 1:Threads.nthreads()]
+                      for _ in 1:NDIMS] for _ in 1:Threads.maxthreadid()]
     rotated_flux_threaded = [Trixi.allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.nthreads()]
+                             for _ in 1:Threads.maxthreadid()]
 
+    # For backwards compatibility with older DGMulti code, the solution is included in the 
+    # top-level cache. TODO: remove once DGMulti refactor is complete and stable.
     cache = (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj,
-             u_values, u_face_values, flux_face_values,
-             aux_values, aux_quad_values, aux_face_values,
-             local_values_threaded, flux_threaded, rotated_flux_threaded)
+             solution_container, u_values, u_face_values, flux_face_values,
+             auxiliary_container, local_values_threaded, flux_threaded,
+             rotated_flux_threaded)
     return cache
 end
 

src/solvers/dgmulti/dg_manifold_covariant.jl

@@ -5,7 +5,8 @@ function Trixi.compute_coefficients!(::Nothing, u, initial_condition, t,
                                      dg::DGMulti, cache)
     md = mesh.md
     rd = dg.basis
-    @unpack u_values, aux_quad_values = cache
+    (; u_values) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
     # evaluate the initial condition at quadrature points
     Trixi.@threaded for i in Trixi.each_quad_node_global(mesh, dg, cache)
         x_node = SVector(getindex.(md.xyzq, i))
@@ -24,7 +25,8 @@ function Trixi.calc_sources!(du, u, t, source_terms,
     rd = dg.basis
     md = mesh.md
     @unpack Pq = rd
-    @unpack u_values, aux_quad_values, local_values_threaded = cache
+    (; u_values, local_values_threaded) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
     Trixi.@threaded for e in Trixi.eachelement(mesh, dg, cache)
         source_values = local_values_threaded[Threads.threadid()]
 
@@ -56,7 +58,9 @@ function Trixi.calc_volume_integral!(du, u,
                                      cache) where {NDIMS_AMBIENT, NDIMS}
     rd = dg.basis
     md = mesh.md
-    (; weak_differentiation_matrices, u_values, aux_quad_values, local_values_threaded) = cache
+    (; weak_differentiation_matrices) = cache
+    (; u_values, local_values_threaded) = cache.solution_container
+    (; aux_quad_values) = cache.auxiliary_container
 
     # interpolate to quadrature points
     Trixi.apply_to_each_field(Trixi.mul_by!(rd.Vq), u_values, u)
@@ -88,7 +92,8 @@ function Trixi.calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeak
     rd = dg.basis
     @unpack mapM, mapP, Jf = md
     @unpack nrstJ, Nfq = rd
-    @unpack u_face_values, flux_face_values, aux_face_values = cache
+    (; u_face_values, flux_face_values) = cache.solution_container
+    (; aux_face_values) = cache.auxiliary_container
     Trixi.@threaded for face_node_index in Trixi.each_face_node_global(mesh, dg, cache)
 
         # inner (idM -> minus) and outer (idP -> plus) indices