Use VectorOfArray in wrap_array for DGMulti solvers (#2150)

* update bound on RecursiveArrayTools to have access to 3.27.3 d * comment out broken precompile statements * adding VoA for DGMulti * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * run `Trixi.rhs!` twice to try to avoid excessive allocations * bump lower compat of RecursiveArrayTools * run Trixi.rhs! twice to reduce allocation count * update l2, linf errors * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update test/test_dgmulti_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Update test/test_dgmulti_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Apply suggestions from code review Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * bump LinearAlgebra lower compat for Downgrade CI * fixing one precompile statement * Revert "fixing one precompile statement" This reverts commit 1b700bc. * change mu::Float64 to mu::Function in elixir * update CI * Formatting suggestions Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * fix threaded CI test values * bump LinearAlgebra compat to fix threaded_legacy tests * Update test/test_threaded.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * unwrap VoA for jacobian computations * unwrap VoA for PlotData1D/PlotData2D * Update src/semidiscretization/semidiscretization.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update src/visualization/types.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Unpack VoA for visualization test * bump Julia compat and ci.yml to v1.10 * Update test/test_visualization.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * fix VoA in test_visualization * Update test/test_dgmulti_3d.jl * Update test/test_dgmulti_3d.jl * uncommenting precompile statement * Auto stash before merge of "jc/wrap_VectorOfArray" and "origin/jc/wrap_VectorOfArray" * removing IfElse.jl from Project.toml (why did I add this?) * mark allocation tests as broken for now * mark allocation tests as broken for now mark all 1D allocated tests as broken for now * update dgmulti tests * remove @test_broken for passing tests * update DGMulti FDSBP 3D CI values * update parabolic dispatch * update parabolic dispatch update parabolic dispatch * Apply suggestions from code review Formatting Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Remove additional Trixi.rhs! call in tests * Reduce number of allowed allocations * remove additional Trixi.rhs! calls and reduce number of allowed allocations * remove more Trixi.rhs! calls * format * clean up parabolic rhs test * Update test/test_parabolic_2d.jl Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> * Update test/test_dgmulti_1d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * use CFL based time stepping to fix positivity issue * format * fix docs * fix test values (use the ones I get locally) * Keep mu as a constant instead of a function * try to fix precompiling of summary callback * fix precompile * try compat "1" for LinearAlgebra * revert precompile stuff * remove #= FIXME comments * Update test/test_parabolic_3d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * Update test/test_parabolic_3d.jl Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> * add news and bump dev version --------- Co-authored-by: github-actions[bot] <41898282+github-actions[bot]@users.noreply.github.com> Co-authored-by: Hendrik Ranocha <ranocha@users.noreply.github.com> Co-authored-by: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com>
trixi-framework · Feb 7, 2025 · 9f2a5be · 9f2a5be
1 parent d107764
commit 9f2a5be
Show file tree

Hide file tree

Showing 16 changed files with 348 additions and 284 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -4,6 +4,23 @@ Trixi.jl follows the interpretation of [semantic versioning (semver)](https://ju
 used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
+
+## Changes when updating to v0.10 from v0.9.x
+
+#### Added
+
+#### Changed
+
+- The numerical solution is wrapped in a `VectorOfArrays` from
+  [RecursiveArrayTools.jl](https://github.com/SciML/RecursiveArrayTools.jl)
+  for `DGMulti` solvers ([#2150]). You can use `Base.parent` to unwrap
+  the original data.
+
+#### Deprecated
+
+#### Removed
+
+
 ## Changes in the v0.9 lifecycle
 
 #### Added
@@ -24,6 +41,7 @@ for human readability.
 
 - The required Julia version is updated to v1.10.
 
+
 ## Changes when updating to v0.9 from v0.8.x
 
 #### Added

diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "Trixi"
 uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
 authors = ["Michael Schlottke-Lakemper <m.schlottke-lakemper@hlrs.de>", "Gregor Gassner <ggassner@uni-koeln.de>", "Hendrik Ranocha <mail@ranocha.de>", "Andrew R. Winters <andrew.ross.winters@liu.se>", "Jesse Chan <jesse.chan@rice.edu>"]
-version = "0.9.18-DEV"
+version = "0.10.0-DEV"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

diff --git a/docs/src/visualization.md b/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

diff --git a/examples/dgmulti_1d/elixir_euler_shu_osher_gauss_shock_capturing.jl b/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)
diff --git a/src/Trixi.jl b/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

diff --git a/src/semidiscretization/semidiscretization.jl b/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.

diff --git a/src/solvers/dgmulti/dg.jl b/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

diff --git a/src/solvers/dgmulti/flux_differencing_gauss_sbp.jl b/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
 

diff --git a/src/visualization/types.jl b/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),

diff --git a/test/test_dgmulti_1d.jl b/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
                         ])