From 3daa745dfd612da1034e5cc8e83cb0da43b0ca1a Mon Sep 17 00:00:00 2001
From: jlchan <jesse.lee.chan@gmail.com>
Date: Mon, 9 Mar 2026 13:46:26 -0500
Subject: [PATCH 1/6] Refactor DGMulti `create_cache` to use shared
 array-allocation helper

Introduce `DGMultiCommonArrays` struct and `allocate_common_dgmulti_arrays`
to eliminate duplicated allocation of `u_values`, `u_face_values`,
`flux_face_values`, and `local_values_threaded` across the `DGMultiWeakForm`
and `DGMultiFluxDiffSBP` `create_cache` methods.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/solvers/dgmulti/dg.jl                | 39 ++++++++++++++++++------
 src/solvers/dgmulti/flux_differencing.jl | 14 +++------
 2 files changed, 34 insertions(+), 19 deletions(-)

diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index 4f302a196a9..c2172efe5ac 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -161,6 +161,30 @@ function set_zero!(du, dg::DGMulti, other_args...)
     return nothing
 end
 
+# Holds arrays shared across most DGMulti cache types:
+# solution values at volume/face quadrature points and thread-local scratch storage.
+struct DGMultiCommonArrays{uType, fType, lType}
+    u_values::uType
+    u_face_values::fType
+    flux_face_values::fType
+    local_values_threaded::lType
+end
+
+# Allocates arrays shared across most DGMulti cache types.
+function allocate_common_dgmulti_arrays(mesh::DGMultiMesh, equations, dg::DGMulti,
+                                        uEltype)
+    rd = dg.basis
+    md = mesh.md
+    nvars = nvariables(equations)
+    u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
+    u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
+    flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
+    local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
+                             for _ in 1:Threads.maxthreadid()]
+    return DGMultiCommonArrays(u_values, u_face_values, flux_face_values,
+                               local_values_threaded)
+end
+
 # Constructs cache variables for both affine and non-affine (curved) DGMultiMeshes
 function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm, RealT,
                       uEltype) where {NDIMS}
@@ -173,12 +197,6 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
     # ∫f(u) * dv/dx_i = ∑_j (Vq*Drst[i])'*diagm(wq)*(rstxyzJ[i,j].*f(Vq*u))
     weak_differentiation_matrices = map(D -> -M \ ((Vq * D)' * Diagonal(wq)), Drst)
 
-    nvars = nvariables(equations)
-
-    # storage for volume quadrature values, face quadrature values, flux values
-    u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
-    u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
-    flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
     if typeof(rd.approximation_type) <:
        Union{SBP, AbstractNonperiodicDerivativeOperator}
         lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators
@@ -186,10 +204,6 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
         lift_scalings = nothing
     end
 
-    # local storage for volume integral and source computations
-    local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.maxthreadid()]
-
     # 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)
@@ -198,11 +212,16 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
     invJ = inv.(rd.Vq * md.J)
 
     # for scaling by curved geometric terms (not used by affine DGMultiMesh)
+    nvars = nvariables(equations)
     flux_threaded = [[allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
                       for _ in 1:NDIMS] for _ in 1:Threads.maxthreadid()]
     rotated_flux_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
                              for _ in 1:Threads.maxthreadid()]
 
+    (; u_values, u_face_values, flux_face_values,
+    local_values_threaded) = allocate_common_dgmulti_arrays(mesh, equations, dg,
+                                                            uEltype)
+
     return (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj,
             u_values, u_face_values, flux_face_values,
             local_values_threaded, flux_threaded, rotated_flux_threaded)
diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl
index 56d4e1f3752..db5f4df8db4 100644
--- a/src/solvers/dgmulti/flux_differencing.jl
+++ b/src/solvers/dgmulti/flux_differencing.jl
@@ -302,21 +302,17 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP,
     # for use with flux differencing schemes
     Qrst_skew = compute_flux_differencing_SBP_matrices(dg)
 
-    # Todo: DGMulti. Factor common storage into a struct (MeshDataCache?) for reuse across solvers?
-    # storage for volume quadrature values, face quadrature values, flux values
-    nvars = nvariables(equations)
-    u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
-    u_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
-    flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
     lift_scalings = rd.wf ./ rd.wq[rd.Fmask] # lift scalings for diag-norm SBP operators
 
-    local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
-                             for _ in 1:Threads.maxthreadid()]
-
+    nvars = nvariables(equations)
     # Use an array of SVectors (chunks of `nvars` are contiguous in memory) to speed up flux differencing
     fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq)
                                for _ in 1:Threads.maxthreadid()]
 
+    (; u_values, u_face_values, flux_face_values,
+    local_values_threaded) = allocate_common_dgmulti_arrays(mesh, equations, dg,
+                                                            uEltype)
+
     return (; md, Qrst_skew, dxidxhatj = md.rstxyzJ,
             invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
             u_values, u_face_values, flux_face_values,

From 70106799e1b25074c29ea0fb3c8426943ed15236 Mon Sep 17 00:00:00 2001
From: jlchan <jesse.lee.chan@gmail.com>
Date: Mon, 9 Mar 2026 16:53:02 -0500
Subject: [PATCH 2/6] Store DGMultiCommonArrays as named field in cache, unpack
 where needed

Replace splatting of `DGMultiCommonArrays` fields into the cache named
tuple with storing the struct directly as `cache.common_arrays`. Update
all access sites across `dg.jl`, `flux_differencing.jl`,
`flux_differencing_gauss_sbp.jl`, `sbp.jl`, and
`analysis_dgmulti.jl` to unpack from `cache.common_arrays`.

Also applies `DGMultiCommonArrays` consistently to all four `create_cache`
variants (`DGMultiWeakForm`, `DGMultiFluxDiffSBP`, `DGMultiFluxDiff`,
`DGMultiFluxDiffPeriodicFDSBP`).

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/callbacks_step/analysis_dgmulti.jl        |  6 +--
 src/solvers/dgmulti/dg.jl                     | 51 ++++++++++---------
 src/solvers/dgmulti/flux_differencing.jl      | 16 +++---
 .../dgmulti/flux_differencing_gauss_sbp.jl    |  5 +-
 src/solvers/dgmulti/sbp.jl                    |  6 +--
 5 files changed, 42 insertions(+), 42 deletions(-)

diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl
index e720d283ffa..58ebccb39ae 100644
--- a/src/callbacks_step/analysis_dgmulti.jl
+++ b/src/callbacks_step/analysis_dgmulti.jl
@@ -10,7 +10,7 @@ function calc_error_norms(func, u, t, analyzer,
                           dg::DGMulti{NDIMS}, cache, cache_analysis) where {NDIMS}
     rd = dg.basis
     md = mesh.md
-    @unpack u_values = cache
+    (; u_values) = cache.common_arrays
 
     # interpolate u to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
@@ -31,7 +31,7 @@ function integrate(func::Func, u, mesh::DGMultiMesh,
                    equations, dg::DGMulti, cache; normalize = true) where {Func}
     rd = dg.basis
     md = mesh.md
-    @unpack u_values = cache
+    (; u_values) = cache.common_arrays
 
     # interpolate u to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
@@ -47,7 +47,7 @@ function analyze(::typeof(entropy_timederivative), du, u, t,
                  mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     rd = dg.basis
     md = mesh.md
-    @unpack u_values = cache
+    (; u_values) = cache.common_arrays
 
     # interpolate u, du to quadrature points
     du_values = similar(u_values) # Todo: DGMulti. Can we move this to the analysis cache somehow?
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index c2172efe5ac..ff5da8775d4 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -42,7 +42,7 @@ end
     end
 end
 
-@inline nelements(dg::DGMulti, cache) = size(cache.u_values)[end]
+@inline nelements(dg::DGMulti, cache) = size(cache.common_arrays.u_values)[end]
 
 """
     eachdim(mesh)
@@ -218,20 +218,17 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
     rotated_flux_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
                              for _ in 1:Threads.maxthreadid()]
 
-    (; u_values, u_face_values, flux_face_values,
-    local_values_threaded) = allocate_common_dgmulti_arrays(mesh, equations, dg,
-                                                            uEltype)
+    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
 
     return (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj,
-            u_values, u_face_values, flux_face_values,
-            local_values_threaded, flux_threaded, rotated_flux_threaded)
+            common_arrays, flux_threaded, rotated_flux_threaded)
 end
 
 function compute_coefficients!(::Nothing, u, initial_condition, t,
                                mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     md = mesh.md
     rd = dg.basis
-    @unpack u_values = cache
+    (; u_values) = cache.common_arrays
 
     # evaluate the initial condition at quadrature points
     @threaded for i in each_quad_node_global(mesh, dg, cache)
@@ -364,7 +361,7 @@ end
 function prolong2interfaces!(cache, u,
                              mesh::DGMultiMesh, equations, dg::DGMulti)
     rd = dg.basis
-    @unpack u_face_values = cache
+    (; u_face_values) = cache.common_arrays
     apply_to_each_field(mul_by!(rd.Vf), u_face_values, u)
 
     return nothing
@@ -377,7 +374,8 @@ end
                                          have_nonconservative_terms::False, equations,
                                          volume_integral::VolumeIntegralWeakForm,
                                          dg::DGMulti, cache)
-    @unpack weak_differentiation_matrices, dxidxhatj, u_values, local_values_threaded = cache
+    @unpack weak_differentiation_matrices, dxidxhatj = cache
+    (; u_values, local_values_threaded) = cache.common_arrays
 
     flux_values = local_values_threaded[Threads.threadid()]
     for i in eachdim(mesh)
@@ -405,7 +403,8 @@ end
                                          have_nonconservative_terms::False, equations,
                                          volume_integral::VolumeIntegralWeakForm,
                                          dg::DGMulti, cache) where {NDIMS}
-    (; weak_differentiation_matrices, dxidxhatj, u_values) = cache
+    (; weak_differentiation_matrices, dxidxhatj) = cache
+    (; u_values) = cache.common_arrays
 
     flux_values = cache.flux_threaded[Threads.threadid()]
     for i in eachdim(mesh)
@@ -446,7 +445,7 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
                                volume_integral::VolumeIntegralWeakForm, dg::DGMulti,
                                cache)
     rd = dg.basis
-    (; u_values) = cache
+    (; u_values) = cache.common_arrays
     # interpolate to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
 
@@ -466,7 +465,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
     @unpack surface_flux = surface_integral
     md = mesh.md
     @unpack mapM, mapP, nxyzJ, Jf = md
-    @unpack u_face_values, flux_face_values = cache
+    (; u_face_values, flux_face_values) = cache.common_arrays
 
     @threaded for face_node_index in each_face_node_global(mesh, dg, cache)
 
@@ -488,7 +487,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
     flux_conservative, flux_nonconservative = surface_integral.surface_flux
     md = mesh.md
     @unpack mapM, mapP, nxyzJ, Jf = md
-    @unpack u_face_values, flux_face_values = cache
+    (; u_face_values, flux_face_values) = cache.common_arrays
 
     @threaded for face_node_index in each_face_node_global(mesh, dg, cache)
 
@@ -522,7 +521,7 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGMulti, cache)
     rd = dg.basis
-    apply_to_each_field(mul_by_accum!(rd.LIFT), du, cache.flux_face_values)
+    apply_to_each_field(mul_by_accum!(rd.LIFT), du, cache.common_arrays.flux_face_values)
 
     return nothing
 end
@@ -532,7 +531,7 @@ function prolong2interfaces!(cache, u,
                              mesh::DGMultiMesh, equations, dg::DGMultiSBP)
     rd = dg.basis
     @unpack Fmask = rd
-    @unpack u_face_values = cache
+    (; u_face_values) = cache.common_arrays
     @threaded for e in eachelement(mesh, dg, cache)
         for (i, fid) in enumerate(Fmask)
             u_face_values[i, e] = u[fid, e]
@@ -548,7 +547,8 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGMultiSBP, cache)
     rd = dg.basis
-    @unpack flux_face_values, lift_scalings = cache
+    (; flux_face_values) = cache.common_arrays
+    @unpack lift_scalings = cache
 
     @threaded for e in eachelement(mesh, dg, cache)
         for i in each_face_node(mesh, dg, cache)
@@ -582,7 +582,7 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                     dg::DGMulti{NDIMS}) where {NDIMS}
     rd = dg.basis
     md = mesh.md
-    @unpack u_face_values, flux_face_values = cache
+    (; u_face_values, flux_face_values) = cache.common_arrays
     @unpack xyzf, nxyzJ, Jf = md
     @unpack surface_flux = dg.surface_integral
 
@@ -617,8 +617,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
         end
     end
 
-    # Note: modifying the values of the reshaped array modifies the values of cache.flux_face_values.
-    # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape.
+    # Note: modifying the values of the reshaped array modifies the values of cache.common_arrays.flux_face_values.
+    # However, we don't have to re-reshape, since cache.common_arrays.flux_face_values still retains its original shape.
 
     return nothing
 end
@@ -641,8 +641,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
     # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300.
     reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ())
 
-    u_face_values = reshape_by_face(cache.u_face_values)
-    flux_face_values = reshape_by_face(cache.flux_face_values)
+    u_face_values = reshape_by_face(cache.common_arrays.u_face_values)
+    flux_face_values = reshape_by_face(cache.common_arrays.flux_face_values)
     Jf = reshape_by_face(md.Jf)
     nxyzJ, xyzf = reshape_by_face.(md.nxyzJ), reshape_by_face.(md.xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix}
 
@@ -667,8 +667,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
         end
     end
 
-    # Note: modifying the values of the reshaped array modifies the values of cache.flux_face_values.
-    # However, we don't have to re-reshape, since cache.flux_face_values still retains its original shape.
+    # Note: modifying the values of the reshaped array modifies the values of cache.common_arrays.flux_face_values.
+    # However, we don't have to re-reshape, since cache.common_arrays.flux_face_values still retains its original shape.
 
     return nothing
 end
@@ -694,7 +694,8 @@ function invert_jacobian!(du, mesh::DGMultiMesh{NDIMS, <:NonAffine}, equations,
                           dg::DGMulti, cache; scaling = -1) where {NDIMS}
     # Vq = interpolation matrix to quadrature points, Pq = quadrature-based L2 projection matrix
     (; Pq, Vq) = dg.basis
-    (; local_values_threaded, invJ) = cache
+    (; invJ) = cache
+    (; local_values_threaded) = cache.common_arrays
 
     @threaded for e in eachelement(mesh, dg, cache)
         du_at_quad_points = local_values_threaded[Threads.threadid()]
@@ -730,7 +731,7 @@ function calc_sources!(du, u, t, source_terms,
     rd = dg.basis
     md = mesh.md
     @unpack Pq = rd
-    @unpack u_values, local_values_threaded = cache
+    (; u_values, local_values_threaded) = cache.common_arrays
     @threaded for e in eachelement(mesh, dg, cache)
         source_values = local_values_threaded[Threads.threadid()]
 
diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl
index db5f4df8db4..787bc5b02f7 100644
--- a/src/solvers/dgmulti/flux_differencing.jl
+++ b/src/solvers/dgmulti/flux_differencing.jl
@@ -309,14 +309,11 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP,
     fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq)
                                for _ in 1:Threads.maxthreadid()]
 
-    (; u_values, u_face_values, flux_face_values,
-    local_values_threaded) = allocate_common_dgmulti_arrays(mesh, equations, dg,
-                                                            uEltype)
+    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
 
     return (; md, Qrst_skew, dxidxhatj = md.rstxyzJ,
             invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
-            u_values, u_face_values, flux_face_values,
-            local_values_threaded, fluxdiff_local_threaded)
+            common_arrays, fluxdiff_local_threaded)
 end
 
 # most general create_cache: works for `DGMultiFluxDiff{<:Polynomial}`
@@ -368,20 +365,23 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiff, RealT,
     interpolated_geometric_terms = map(x -> [Vq; Vf] * x, mesh.md.rstxyzJ)
     J = rd.Vq * md.J
 
+    common_arrays = DGMultiCommonArrays(u_values, u_face_values, flux_face_values,
+                                        local_values_threaded)
+
     return (; md, Qrst_skew, VhP, Ph,
             invJ = inv.(J), dxidxhatj = interpolated_geometric_terms,
             entropy_var_values, projected_entropy_var_values,
             entropy_projected_u_values,
-            u_values, u_face_values, flux_face_values,
-            local_values_threaded, fluxdiff_local_threaded, rhs_local_threaded)
+            common_arrays, fluxdiff_local_threaded, rhs_local_threaded)
 end
 
 # TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function.
 function entropy_projection!(cache, u, mesh::DGMultiMesh, equations, dg::DGMulti)
     rd = dg.basis
     @unpack Vq = rd
-    @unpack VhP, entropy_var_values, u_values = cache
+    @unpack VhP, entropy_var_values = cache
     @unpack projected_entropy_var_values, entropy_projected_u_values = cache
+    (; u_values) = cache.common_arrays
 
     apply_to_each_field(mul_by!(Vq), u_values, u)
 
diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
index 5185a7dbec5..61f2e74ed01 100644
--- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
+++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
@@ -454,9 +454,10 @@ function entropy_projection!(cache, u, mesh::DGMultiMesh, equations,
                              dg::DGMultiFluxDiff{<:GaussSBP})
     rd = dg.basis
     @unpack Vq = rd
-    @unpack VhP, entropy_var_values, u_values = cache
+    @unpack VhP, entropy_var_values = cache
     @unpack projected_entropy_var_values, entropy_projected_u_values = cache
     @unpack interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_face = cache
+    (; u_values) = cache.common_arrays
 
     @threaded for e in eachelement(mesh, dg, cache)
         apply_to_each_field(mul_by!(interp_matrix_lobatto_to_gauss),
@@ -501,7 +502,7 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
         # applies LIFT matrix, output is stored at Gauss nodes
         gauss_volume_local = gauss_volume_local_threaded[Threads.threadid()]
         apply_to_each_field(mul_by!(gauss_LIFT), gauss_volume_local,
-                            view(cache.flux_face_values, :, e))
+                            view(cache.common_arrays.flux_face_values, :, e))
 
         for i in eachindex(gauss_volume_local)
             du[i, e] = du[i, e] + gauss_volume_local[i]
diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl
index a6f7fbf844a..160c5d7e23e 100644
--- a/src/solvers/dgmulti/sbp.jl
+++ b/src/solvers/dgmulti/sbp.jl
@@ -179,10 +179,8 @@ function create_cache(mesh::DGMultiMesh, equations,
                       dg::DGMultiFluxDiffPeriodicFDSBP, RealT, uEltype)
     md = mesh.md
 
-    # storage for volume quadrature values, face quadrature values, flux values
-    nvars = nvariables(equations)
-    u_values = allocate_nested_array(uEltype, nvars, size(md.xq), dg)
-    return (; u_values, invJ = inv.(md.J))
+    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
+    return (; common_arrays, invJ = inv.(md.J))
 end
 
 # Specialize calc_volume_integral for periodic SBP operators (assumes the operator is sparse).

From e080a740bb48ae87b16f0e0ac21b883322593ec9 Mon Sep 17 00:00:00 2001
From: Jesse Chan <1156048+jlchan@users.noreply.github.com>
Date: Mon, 9 Mar 2026 16:54:49 -0500
Subject: [PATCH 3/6] Apply suggestions from code review

Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com>
---
 src/solvers/dgmulti/dg.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index ff5da8775d4..e7680b14010 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -521,7 +521,8 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGMulti, cache)
     rd = dg.basis
-    apply_to_each_field(mul_by_accum!(rd.LIFT), du, cache.common_arrays.flux_face_values)
+    apply_to_each_field(mul_by_accum!(rd.LIFT), du,
+                        cache.common_arrays.flux_face_values)
 
     return nothing
 end

From d00b51584489650d039fbdbbe7b07ef94b8f06be Mon Sep 17 00:00:00 2001
From: jlchan <jesse.lee.chan@gmail.com>
Date: Mon, 9 Mar 2026 17:19:42 -0500
Subject: [PATCH 4/6] Fix DGMultiCommonArrays type parameter for mismatched
 face array types
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

`u_face_values` in the `DGMultiFluxDiff` cache is a view into
`entropy_projected_u_values`, while `flux_face_values = similar(...)`
allocates a fresh array — these have different concrete types.
Split the single `fType` parameter into separate `ufType` and `ffType`
to accommodate this.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
---
 src/solvers/dgmulti/dg.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index e7680b14010..2149b7f3ede 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -163,10 +163,10 @@ end
 
 # Holds arrays shared across most DGMulti cache types:
 # solution values at volume/face quadrature points and thread-local scratch storage.
-struct DGMultiCommonArrays{uType, fType, lType}
+struct DGMultiCommonArrays{uType, ufType, ffType, lType}
     u_values::uType
-    u_face_values::fType
-    flux_face_values::fType
+    u_face_values::ufType
+    flux_face_values::ffType
     local_values_threaded::lType
 end
 

From ea3f5038de2d87fbd715caffa9555e17a23684ad Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Tue, 10 Mar 2026 13:06:02 -0500
Subject: [PATCH 5/6] renaming

---
 src/callbacks_step/analysis_dgmulti.jl        |  6 +-
 src/solvers/dgmulti/dg.jl                     | 56 ++++++++++---------
 src/solvers/dgmulti/flux_differencing.jl      | 14 +++--
 .../dgmulti/flux_differencing_gauss_sbp.jl    |  4 +-
 src/solvers/dgmulti/sbp.jl                    |  5 +-
 5 files changed, 45 insertions(+), 40 deletions(-)

diff --git a/src/callbacks_step/analysis_dgmulti.jl b/src/callbacks_step/analysis_dgmulti.jl
index 58ebccb39ae..281301499c6 100644
--- a/src/callbacks_step/analysis_dgmulti.jl
+++ b/src/callbacks_step/analysis_dgmulti.jl
@@ -10,7 +10,7 @@ function calc_error_norms(func, u, t, analyzer,
                           dg::DGMulti{NDIMS}, cache, cache_analysis) where {NDIMS}
     rd = dg.basis
     md = mesh.md
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     # interpolate u to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
@@ -31,7 +31,7 @@ function integrate(func::Func, u, mesh::DGMultiMesh,
                    equations, dg::DGMulti, cache; normalize = true) where {Func}
     rd = dg.basis
     md = mesh.md
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     # interpolate u to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
@@ -47,7 +47,7 @@ function analyze(::typeof(entropy_timederivative), du, u, t,
                  mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     rd = dg.basis
     md = mesh.md
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     # interpolate u, du to quadrature points
     du_values = similar(u_values) # Todo: DGMulti. Can we move this to the analysis cache somehow?
diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl
index 2149b7f3ede..3d7aba9536f 100644
--- a/src/solvers/dgmulti/dg.jl
+++ b/src/solvers/dgmulti/dg.jl
@@ -42,7 +42,7 @@ end
     end
 end
 
-@inline nelements(dg::DGMulti, cache) = size(cache.common_arrays.u_values)[end]
+@inline nelements(dg::DGMulti, cache) = size(cache.solution_container.u_values)[end]
 
 """
     eachdim(mesh)
@@ -163,7 +163,7 @@ end
 
 # Holds arrays shared across most DGMulti cache types:
 # solution values at volume/face quadrature points and thread-local scratch storage.
-struct DGMultiCommonArrays{uType, ufType, ffType, lType}
+struct DGMultiSolutionContainer{uType, ufType, ffType, lType}
     u_values::uType
     u_face_values::ufType
     flux_face_values::ffType
@@ -171,8 +171,9 @@ struct DGMultiCommonArrays{uType, ufType, ffType, lType}
 end
 
 # Allocates arrays shared across most DGMulti cache types.
-function allocate_common_dgmulti_arrays(mesh::DGMultiMesh, equations, dg::DGMulti,
-                                        uEltype)
+function initialize_dgmulti_solution_container(mesh::DGMultiMesh, equations,
+                                               dg::DGMulti,
+                                               uEltype)
     rd = dg.basis
     md = mesh.md
     nvars = nvariables(equations)
@@ -181,8 +182,8 @@ function allocate_common_dgmulti_arrays(mesh::DGMultiMesh, equations, dg::DGMult
     flux_face_values = allocate_nested_array(uEltype, nvars, size(md.xf), dg)
     local_values_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
                              for _ in 1:Threads.maxthreadid()]
-    return DGMultiCommonArrays(u_values, u_face_values, flux_face_values,
-                               local_values_threaded)
+    return DGMultiSolutionContainer(u_values, u_face_values, flux_face_values,
+                                    local_values_threaded)
 end
 
 # Constructs cache variables for both affine and non-affine (curved) DGMultiMeshes
@@ -218,17 +219,18 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
     rotated_flux_threaded = [allocate_nested_array(uEltype, nvars, (rd.Nq,), dg)
                              for _ in 1:Threads.maxthreadid()]
 
-    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
+    solution_container = initialize_dgmulti_solution_container(mesh, equations, dg,
+                                                               uEltype)
 
     return (; md, weak_differentiation_matrices, lift_scalings, invJ, dxidxhatj,
-            common_arrays, flux_threaded, rotated_flux_threaded)
+            solution_container, flux_threaded, rotated_flux_threaded)
 end
 
 function compute_coefficients!(::Nothing, u, initial_condition, t,
                                mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     md = mesh.md
     rd = dg.basis
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     # evaluate the initial condition at quadrature points
     @threaded for i in each_quad_node_global(mesh, dg, cache)
@@ -361,7 +363,7 @@ end
 function prolong2interfaces!(cache, u,
                              mesh::DGMultiMesh, equations, dg::DGMulti)
     rd = dg.basis
-    (; u_face_values) = cache.common_arrays
+    (; u_face_values) = cache.solution_container
     apply_to_each_field(mul_by!(rd.Vf), u_face_values, u)
 
     return nothing
@@ -375,7 +377,7 @@ end
                                          volume_integral::VolumeIntegralWeakForm,
                                          dg::DGMulti, cache)
     @unpack weak_differentiation_matrices, dxidxhatj = cache
-    (; u_values, local_values_threaded) = cache.common_arrays
+    (; u_values, local_values_threaded) = cache.solution_container
 
     flux_values = local_values_threaded[Threads.threadid()]
     for i in eachdim(mesh)
@@ -404,7 +406,7 @@ end
                                          volume_integral::VolumeIntegralWeakForm,
                                          dg::DGMulti, cache) where {NDIMS}
     (; weak_differentiation_matrices, dxidxhatj) = cache
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     flux_values = cache.flux_threaded[Threads.threadid()]
     for i in eachdim(mesh)
@@ -445,7 +447,7 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
                                volume_integral::VolumeIntegralWeakForm, dg::DGMulti,
                                cache)
     rd = dg.basis
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
     # interpolate to quadrature points
     apply_to_each_field(mul_by!(rd.Vq), u_values, u)
 
@@ -465,7 +467,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
     @unpack surface_flux = surface_integral
     md = mesh.md
     @unpack mapM, mapP, nxyzJ, Jf = md
-    (; u_face_values, flux_face_values) = cache.common_arrays
+    (; u_face_values, flux_face_values) = cache.solution_container
 
     @threaded for face_node_index in each_face_node_global(mesh, dg, cache)
 
@@ -487,7 +489,7 @@ function calc_interface_flux!(cache, surface_integral::SurfaceIntegralWeakForm,
     flux_conservative, flux_nonconservative = surface_integral.surface_flux
     md = mesh.md
     @unpack mapM, mapP, nxyzJ, Jf = md
-    (; u_face_values, flux_face_values) = cache.common_arrays
+    (; u_face_values, flux_face_values) = cache.solution_container
 
     @threaded for face_node_index in each_face_node_global(mesh, dg, cache)
 
@@ -522,7 +524,7 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
                                 dg::DGMulti, cache)
     rd = dg.basis
     apply_to_each_field(mul_by_accum!(rd.LIFT), du,
-                        cache.common_arrays.flux_face_values)
+                        cache.solution_container.flux_face_values)
 
     return nothing
 end
@@ -532,7 +534,7 @@ function prolong2interfaces!(cache, u,
                              mesh::DGMultiMesh, equations, dg::DGMultiSBP)
     rd = dg.basis
     @unpack Fmask = rd
-    (; u_face_values) = cache.common_arrays
+    (; u_face_values) = cache.solution_container
     @threaded for e in eachelement(mesh, dg, cache)
         for (i, fid) in enumerate(Fmask)
             u_face_values[i, e] = u[fid, e]
@@ -548,7 +550,7 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
                                 surface_integral::SurfaceIntegralWeakForm,
                                 dg::DGMultiSBP, cache)
     rd = dg.basis
-    (; flux_face_values) = cache.common_arrays
+    (; flux_face_values) = cache.solution_container
     @unpack lift_scalings = cache
 
     @threaded for e in eachelement(mesh, dg, cache)
@@ -583,7 +585,7 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
                                     dg::DGMulti{NDIMS}) where {NDIMS}
     rd = dg.basis
     md = mesh.md
-    (; u_face_values, flux_face_values) = cache.common_arrays
+    (; u_face_values, flux_face_values) = cache.solution_container
     @unpack xyzf, nxyzJ, Jf = md
     @unpack surface_flux = dg.surface_integral
 
@@ -618,8 +620,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
         end
     end
 
-    # Note: modifying the values of the reshaped array modifies the values of cache.common_arrays.flux_face_values.
-    # However, we don't have to re-reshape, since cache.common_arrays.flux_face_values still retains its original shape.
+    # Note: modifying the values of the reshaped array modifies the values of cache.solution_container.flux_face_values.
+    # However, we don't have to re-reshape, since cache.solution_container.flux_face_values still retains its original shape.
 
     return nothing
 end
@@ -642,8 +644,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
     # https://github.com/JuliaLang/julia/issues/36313#issuecomment-782336300.
     reshape_by_face(u) = Base.ReshapedArray(u, (num_pts_per_face, num_faces_total), ())
 
-    u_face_values = reshape_by_face(cache.common_arrays.u_face_values)
-    flux_face_values = reshape_by_face(cache.common_arrays.flux_face_values)
+    u_face_values = reshape_by_face(cache.solution_container.u_face_values)
+    flux_face_values = reshape_by_face(cache.solution_container.flux_face_values)
     Jf = reshape_by_face(md.Jf)
     nxyzJ, xyzf = reshape_by_face.(md.nxyzJ), reshape_by_face.(md.xyzf) # broadcast over nxyzJ::NTuple{NDIMS,Matrix}
 
@@ -668,8 +670,8 @@ function calc_single_boundary_flux!(cache, t, boundary_condition, boundary_key,
         end
     end
 
-    # Note: modifying the values of the reshaped array modifies the values of cache.common_arrays.flux_face_values.
-    # However, we don't have to re-reshape, since cache.common_arrays.flux_face_values still retains its original shape.
+    # Note: modifying the values of the reshaped array modifies the values of cache.solution_container.flux_face_values.
+    # However, we don't have to re-reshape, since cache.solution_container.flux_face_values still retains its original shape.
 
     return nothing
 end
@@ -696,7 +698,7 @@ function invert_jacobian!(du, mesh::DGMultiMesh{NDIMS, <:NonAffine}, equations,
     # Vq = interpolation matrix to quadrature points, Pq = quadrature-based L2 projection matrix
     (; Pq, Vq) = dg.basis
     (; invJ) = cache
-    (; local_values_threaded) = cache.common_arrays
+    (; local_values_threaded) = cache.solution_container
 
     @threaded for e in eachelement(mesh, dg, cache)
         du_at_quad_points = local_values_threaded[Threads.threadid()]
@@ -732,7 +734,7 @@ function calc_sources!(du, u, t, source_terms,
     rd = dg.basis
     md = mesh.md
     @unpack Pq = rd
-    (; u_values, local_values_threaded) = cache.common_arrays
+    (; u_values, local_values_threaded) = cache.solution_container
     @threaded for e in eachelement(mesh, dg, cache)
         source_values = local_values_threaded[Threads.threadid()]
 
diff --git a/src/solvers/dgmulti/flux_differencing.jl b/src/solvers/dgmulti/flux_differencing.jl
index 787bc5b02f7..9b3bc211d1b 100644
--- a/src/solvers/dgmulti/flux_differencing.jl
+++ b/src/solvers/dgmulti/flux_differencing.jl
@@ -309,11 +309,12 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiffSBP,
     fluxdiff_local_threaded = [zeros(SVector{nvars, uEltype}, rd.Nq)
                                for _ in 1:Threads.maxthreadid()]
 
-    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
+    solution_container = initialize_dgmulti_solution_container(mesh, equations, dg,
+                                                               uEltype)
 
     return (; md, Qrst_skew, dxidxhatj = md.rstxyzJ,
             invJ = inv.(md.J), lift_scalings, inv_wq = inv.(rd.wq),
-            common_arrays, fluxdiff_local_threaded)
+            solution_container, fluxdiff_local_threaded)
 end
 
 # most general create_cache: works for `DGMultiFluxDiff{<:Polynomial}`
@@ -365,14 +366,15 @@ function create_cache(mesh::DGMultiMesh, equations, dg::DGMultiFluxDiff, RealT,
     interpolated_geometric_terms = map(x -> [Vq; Vf] * x, mesh.md.rstxyzJ)
     J = rd.Vq * md.J
 
-    common_arrays = DGMultiCommonArrays(u_values, u_face_values, flux_face_values,
-                                        local_values_threaded)
+    solution_container = DGMultiSolutionContainer(u_values, u_face_values,
+                                                  flux_face_values,
+                                                  local_values_threaded)
 
     return (; md, Qrst_skew, VhP, Ph,
             invJ = inv.(J), dxidxhatj = interpolated_geometric_terms,
             entropy_var_values, projected_entropy_var_values,
             entropy_projected_u_values,
-            common_arrays, fluxdiff_local_threaded, rhs_local_threaded)
+            solution_container, fluxdiff_local_threaded, rhs_local_threaded)
 end
 
 # TODO: DGMulti. Address hard-coding of `entropy2cons!` and `cons2entropy!` for this function.
@@ -381,7 +383,7 @@ function entropy_projection!(cache, u, mesh::DGMultiMesh, equations, dg::DGMulti
     @unpack Vq = rd
     @unpack VhP, entropy_var_values = cache
     @unpack projected_entropy_var_values, entropy_projected_u_values = cache
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     apply_to_each_field(mul_by!(Vq), u_values, u)
 
diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
index 61f2e74ed01..48d75938d86 100644
--- a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
+++ b/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl
@@ -457,7 +457,7 @@ function entropy_projection!(cache, u, mesh::DGMultiMesh, equations,
     @unpack VhP, entropy_var_values = cache
     @unpack projected_entropy_var_values, entropy_projected_u_values = cache
     @unpack interp_matrix_lobatto_to_gauss, interp_matrix_gauss_to_face = cache
-    (; u_values) = cache.common_arrays
+    (; u_values) = cache.solution_container
 
     @threaded for e in eachelement(mesh, dg, cache)
         apply_to_each_field(mul_by!(interp_matrix_lobatto_to_gauss),
@@ -502,7 +502,7 @@ function calc_surface_integral!(du, u, mesh::DGMultiMesh, equations,
         # applies LIFT matrix, output is stored at Gauss nodes
         gauss_volume_local = gauss_volume_local_threaded[Threads.threadid()]
         apply_to_each_field(mul_by!(gauss_LIFT), gauss_volume_local,
-                            view(cache.common_arrays.flux_face_values, :, e))
+                            view(cache.solution_container.flux_face_values, :, e))
 
         for i in eachindex(gauss_volume_local)
             du[i, e] = du[i, e] + gauss_volume_local[i]
diff --git a/src/solvers/dgmulti/sbp.jl b/src/solvers/dgmulti/sbp.jl
index 160c5d7e23e..bb0aea56b28 100644
--- a/src/solvers/dgmulti/sbp.jl
+++ b/src/solvers/dgmulti/sbp.jl
@@ -179,8 +179,9 @@ function create_cache(mesh::DGMultiMesh, equations,
                       dg::DGMultiFluxDiffPeriodicFDSBP, RealT, uEltype)
     md = mesh.md
 
-    common_arrays = allocate_common_dgmulti_arrays(mesh, equations, dg, uEltype)
-    return (; common_arrays, invJ = inv.(md.J))
+    solution_container = initialize_dgmulti_solution_container(mesh, equations, dg,
+                                                               uEltype)
+    return (; solution_container, invJ = inv.(md.J))
 end
 
 # Specialize calc_volume_integral for periodic SBP operators (assumes the operator is sparse).

From 05caa50d1d80768ef1b506552e2fe779e147c766 Mon Sep 17 00:00:00 2001
From: Jesse Chan <jesse.lee.chan@gmail.com>
Date: Wed, 11 Mar 2026 10:32:43 -0500
Subject: [PATCH 6/6] fix unpack

---
 src/solvers/dgmulti/volume_integral_adaptive.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/solvers/dgmulti/volume_integral_adaptive.jl b/src/solvers/dgmulti/volume_integral_adaptive.jl
index 0b29fdc3174..2fc2eaa079f 100644
--- a/src/solvers/dgmulti/volume_integral_adaptive.jl
+++ b/src/solvers/dgmulti/volume_integral_adaptive.jl
@@ -28,7 +28,7 @@ function create_cache(mesh::DGMultiMesh, equations,
     nvars = nvariables(equations)
 
     # Required for entropy change computation (`entropy_change_reference_element`)
-    du_values = similar(cache_FD.u_values)
+    du_values = similar(cache_FD.solution_container.u_values)
 
     # Thread-local buffer for face interpolation, which is required
     # for computation of entropy potential at interpolated face nodes
@@ -54,7 +54,7 @@ function calc_volume_integral!(du, u, mesh::DGMultiMesh,
     @unpack maximum_entropy_increase = volume_integral.indicator
 
     # For weak form integral
-    @unpack u_values = cache
+    @unpack u_values = cache.solution_container
 
     # For entropy production computation
     rd = dg.basis