Use VectorOfArray in wrap_array for DGMulti solvers

docs/src/visualization.md

@@ -174,7 +174,7 @@ julia> compute_vorticity(velocity, semi) =
          compute_vorticity(velocity, Trixi.mesh_equations_solver_cache(semi)...);
 
 julia> function get_velocity(sol)
-         rho, rhou, rhov, E = StructArrays.components(sol.u[end])
+         rho, rhou, rhov, E = StructArrays.components(Base.parent(sol.u[end]))
          v1 = rhou ./ rho
          v2 = rhov ./ rho
          return v1, v2

examples/dgmulti_1d/elixir_euler_shu_osher_gauss_shock_capturing.jl

@@ -5,8 +5,8 @@ gamma_gas = 1.4
 equations = CompressibleEulerEquations1D(gamma_gas)
 
 ###############################################################################
-# setup the GSBP DG discretization that uses the Gauss operators from 
-# Chan, Del Rey Fernandez, Carpenter (2019). 
+# setup the GSBP DG discretization that uses the Gauss operators from
+# Chan, Del Rey Fernandez, Carpenter (2019).
 # [https://doi.org/10.1137/18M1209234](https://doi.org/10.1137/18M1209234)
 
 # Shu-Osher initial condition for 1D compressible Euler equations
@@ -19,8 +19,7 @@ function initial_condition_shu_osher(x, t, equations::CompressibleEulerEquations
     v_left = 4 * sqrt(35) / 9
     p_left = 31 / 3
 
-    # Replaced v_right = 0 to v_right = 0.1 to avoid positivity issues.
-    v_right = 0.1
+    v_right = 0.0
     p_right = 1.0
 
     rho = ifelse(x[1] > x0, 1 + 1 / 5 * sin(5 * x[1]), rho_left)
@@ -82,12 +81,14 @@ summary_callback = SummaryCallback()
 analysis_callback = AnalysisCallback(semi, interval = 100, uEltype = real(dg))
 
 # handles the re-calculation of the maximum Δt after each time step
-stepsize_callback = StepsizeCallback(cfl = 0.1)
+stepsize_callback = StepsizeCallback(cfl = 0.2)
 
 # collect all callbacks such that they can be passed to the ODE solver
 callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)
 
 # ###############################################################################
 # # run the simulation
 
-sol = solve(ode, SSPRK43(), adaptive = true, callback = callbacks, save_everystep = false)
+sol = solve(ode, SSPRK43(), adaptive = false,
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            callback = callbacks, save_everystep = false)

src/Trixi.jl

@@ -60,6 +60,7 @@ using OffsetArrays: OffsetArray, OffsetVector
 using P4est
 using T8code
 using RecipesBase: RecipesBase
+using RecursiveArrayTools: VectorOfArray
 using Requires: @require
 using Static: Static, One, True, False
 @reexport using StaticArrays: SVector

src/semidiscretization/semidiscretization.jl

@@ -277,6 +277,11 @@ function _jacobian_ad_forward(semi, t0, u0_ode, du_ode, config)
     return J
 end
 
+# unpack u if it is wrapped in VectorOfArray (mainly for DGMulti solvers)
+jacobian_ad_forward(semi::AbstractSemidiscretization, t0, u0_ode::VectorOfArray) = jacobian_ad_forward(semi,
+                                                                                                       t0,
+                                                                                                       parent(u0_ode))
+
 # This version is specialized to `StructArray`s used by some `DGMulti` solvers.
 # We need to convert the numerical solution vectors since ForwardDiff cannot
 # handle arrays of `SVector`s.

src/solvers/dgmulti/dg.jl

@@ -128,6 +128,14 @@ end
 # interface with semidiscretization_hyperbolic
 wrap_array(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
 wrap_array_native(u_ode, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = u_ode
+
+# used to initialize `u_ode` in `semidiscretize`
+function allocate_coefficients(mesh::DGMultiMesh, equations, dg::DGMulti, cache)
+    return VectorOfArray(allocate_nested_array(real(dg), nvariables(equations),
+                                               size(mesh.md.x), dg))
+end
+wrap_array(u_ode::VectorOfArray, mesh::DGMultiMesh, equations, dg::DGMulti, cache) = parent(u_ode)
+
 function digest_boundary_conditions(boundary_conditions::NamedTuple{Keys, ValueTypes},
                                     mesh::DGMultiMesh,
                                     dg::DGMulti,
@@ -199,10 +207,6 @@ function create_cache(mesh::DGMultiMesh{NDIMS}, equations, dg::DGMultiWeakForm,
             local_values_threaded, flux_threaded, rotated_flux_threaded)
 end
 
-function allocate_coefficients(mesh::DGMultiMesh, equations, dg::DGMulti, cache)
-    return allocate_nested_array(real(dg), nvariables(equations), size(mesh.md.x), dg)
-end
-
 function compute_coefficients!(u, initial_condition, t,
                                mesh::DGMultiMesh, equations, dg::DGMulti, cache)
     md = mesh.md

src/solvers/dgmulti/flux_differencing_gauss_sbp.jl

@@ -176,7 +176,7 @@ end
 @inline function tensor_product_gauss_face_operator!(out::AbstractVector,
                                                      A::TensorProductGaussFaceOperator{2, Interpolation},
                                                      x_in::AbstractVector)
-#! format: on                                                     
+#! format: on
     (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A
     (; nnodes_1d) = A
 
@@ -215,7 +215,7 @@ end
 @inline function tensor_product_gauss_face_operator!(out::AbstractVector,
                                                      A::TensorProductGaussFaceOperator{3, Interpolation},
                                                      x::AbstractVector)
-#! format: on                                                     
+#! format: on
     (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A
     (; nnodes_1d) = A
 
@@ -268,7 +268,7 @@ end
 @inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector,
                                                      A::TensorProductGaussFaceOperator{2, Projection{ApplyFaceWeights}},
                                                      x::AbstractVector) where {ApplyFaceWeights}
-#! format: on                                                     
+#! format: on
     (; interp_matrix_gauss_to_face_1d, face_indices_tensor_product) = A
     (; inv_volume_weights_1d, nnodes_1d) = A
 
@@ -322,7 +322,7 @@ end
 @inline function tensor_product_gauss_face_operator!(out_vec::AbstractVector,
                                                      A::TensorProductGaussFaceOperator{3, Projection{ApplyFaceWeights}},
                                                      x::AbstractVector) where {ApplyFaceWeights}
-#! format: on                                                                               
+#! format: on
     @unpack interp_matrix_gauss_to_face_1d, face_indices_tensor_product = A
     @unpack inv_volume_weights_1d, nnodes_1d, nfaces = A
 

src/visualization/types.jl

@@ -635,6 +635,20 @@ function PlotData1D(u, mesh::TreeMesh, equations, solver, cache;
                       orientation_x)
 end
 
+# unwrap u if it is VectorOfArray 
+PlotData1D(u::VectorOfArray, mesh, equations, dg::DGMulti{1}, cache; kwargs...) = PlotData1D(parent(u),
+                                                                                             mesh,
+                                                                                             equations,
+                                                                                             dg,
+                                                                                             cache;
+                                                                                             kwargs...)
+PlotData2D(u::VectorOfArray, mesh, equations, dg::DGMulti{2}, cache; kwargs...) = PlotData2D(parent(u),
+                                                                                             mesh,
+                                                                                             equations,
+                                                                                             dg,
+                                                                                             cache;
+                                                                                             kwargs...)
+
 function PlotData1D(u, mesh, equations, solver, cache;
                     solution_variables = nothing, nvisnodes = nothing,
                     reinterpolate = default_reinterpolate(solver),

test/test_dgmulti_1d.jl

@@ -74,14 +74,14 @@ end
                                  "elixir_euler_shu_osher_gauss_shock_capturing.jl"),
                         cells_per_dimension=(64,), tspan=(0.0, 1.0),
                         l2=[
-                            1.673813320412685,
-                            5.980737909458242,
-                            21.587822949251173
+                            1.6967151731067875,
+                            6.018445633981826,
+                            21.77425594743242
                         ],
                         linf=[
-                            3.1388039126918064,
-                            10.630952212105246,
-                            37.682826521024865
+                            3.2229876650556477,
+                            10.702690533393842,
+                            38.37424900889908
                         ])
     # Ensure that we do not have excessive memory allocations
     # (e.g., from type instabilities)
@@ -207,14 +207,14 @@ end
                         cells_per_dimension=(8,),
                         approximation_type=SBP(),
                         l2=[
-                            3.03001101100507e-6,
-                            1.692177335948727e-5,
-                            3.002634351734614e-16,
-                            1.1636653574178203e-15
+                            3.0300196635805022e-6,
+                            1.6921833812545857e-5,
+                            2.9844594164368975e-16,
+                            1.1012004949980629e-15
                         ],
                         linf=[
-                            1.2043401988570679e-5,
-                            5.346847010329059e-5,
+                            1.2043309307818717e-5,
+                            5.346754311919e-5,
                             9.43689570931383e-16,
                             2.220446049250313e-15
                         ])