diff --git a/src/callbacks_step/analysis_covariant.jl b/src/callbacks_step/analysis_covariant.jl index a3b311c0..9cc247c2 100644 --- a/src/callbacks_step/analysis_covariant.jl +++ b/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) diff --git a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl b/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl index baabb9d7..694549ae 100644 --- a/src/callbacks_step/save_solution_2d_manifold_in_3d_cartesian.jl +++ b/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)) diff --git a/src/callbacks_step/save_solution_covariant.jl b/src/callbacks_step/save_solution_covariant.jl index 677ff8a4..5b160214 100644 --- a/src/callbacks_step/save_solution_covariant.jl +++ b/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 diff --git a/src/callbacks_step/stepsize_dg2d.jl b/src/callbacks_step/stepsize_dg2d.jl index b5e5069c..6a36a07a 100644 --- a/src/callbacks_step/stepsize_dg2d.jl +++ b/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) diff --git a/src/solvers/dgmulti/dg.jl b/src/solvers/dgmulti/dg.jl index 65fa3e4f..96a71582 100644 --- a/src/solvers/dgmulti/dg.jl +++ b/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 diff --git a/src/solvers/dgmulti/dg_manifold_covariant.jl b/src/solvers/dgmulti/dg_manifold_covariant.jl index 8d597ee0..c923a2a3 100644 --- a/src/solvers/dgmulti/dg_manifold_covariant.jl +++ b/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