diff --git a/NEWS.md b/NEWS.md index dc498557381..dc811bf1e19 100644 --- 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 index 53fae1bd842..95d3c42d395 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Trixi" uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb" authors = ["Michael Schlottke-Lakemper ", "Gregor Gassner ", "Hendrik Ranocha ", "Andrew R. Winters ", "Jesse Chan "] -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 index 95f7ed6b684..694c6039876 100644 --- 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 index 79c92656176..fac9b24113e 100644 --- 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,7 +81,7 @@ 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) @@ -90,4 +89,6 @@ 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 index 6cfc9a46fe6..ffe759d0750 100644 --- 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 index dd5c3c4791d..f41c7ea4a7f 100644 --- 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 index 2d588c5c79d..a48a6288eba 100644 --- 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 index 5be36a62d7b..e2a534839d9 100644 --- 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 index e291b314bc2..4a3d384c209 100644 --- 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 index 1c3cd604df1..d35b88d88a8 100644 --- 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 ]) diff --git a/test/test_dgmulti_2d.jl b/test/test_dgmulti_2d.jl index 2e41591d52c..75584940013 100644 --- a/test/test_dgmulti_2d.jl +++ b/test/test_dgmulti_2d.jl @@ -31,8 +31,8 @@ isdir(outdir) && rm(outdir, recursive = true) 0.002062399194485476, 0.004897700392503701 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -59,8 +59,8 @@ end 0.013593978324845324, 0.03270995869587523 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -87,8 +87,8 @@ end 0.0013568139808290969, 0.0032249020004324613 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -115,8 +115,8 @@ end 0.01396529873508534, 0.04227683691807904 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -144,8 +144,8 @@ end 0.05321027922608157, 0.13392025411844033 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -174,8 +174,8 @@ end 0.012674028479251254, 0.02210545278615017 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -187,15 +187,19 @@ end @trixi_testset "elixir_euler_bilinear.jl (Bilinear quadrilateral elements, SBP, flux differencing)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_bilinear.jl"), l2=[ - 1.0259432774540821e-5, 9.014087689495575e-6, - 9.01408768888544e-6, 2.738953324859446e-5 + 1.0267413589968656e-5, + 9.03069720963081e-6, + 9.030697209721065e-6, + 2.7436672091049314e-5 ], linf=[ - 7.362605996297233e-5, 6.874189724781488e-5, - 6.874189703509614e-5, 0.00019124355334110277 + 7.36251369879426e-5, + 6.874041557969335e-5, + 6.874041552329402e-5, + 0.00019123932693609902 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -207,15 +211,19 @@ end @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, SBP, flux differencing)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"), l2=[ - 1.7204593127904542e-5, 1.5921547179522804e-5, - 1.5921547180107928e-5, 4.894071422525737e-5 + 1.7209164346836478e-5, + 1.5928649356474767e-5, + 1.5928649356802847e-5, + 4.8963394546089164e-5 ], linf=[ - 0.00010525416937667842, 0.00010003778102718464, - 0.00010003778071832059, 0.0003642628211952825 + 0.00010525404319428056, + 0.00010003768703326088, + 0.00010003768694910598, + 0.0003642622844113319 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -229,20 +237,20 @@ end approximation_type=GaussSBP(), surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)), l2=[ - 3.4666312079259457e-6, - 3.4392774480368986e-6, - 3.439277447953705e-6, - 1.0965598424665836e-5 + 3.4664508443541302e-6, + 3.4389354928807557e-6, + 3.438935492692069e-6, + 1.0965259031107001e-5 ], linf=[ - 1.1327280377004811e-5, - 1.1343911926253725e-5, - 1.1343911906935844e-5, - 3.679582619220412e-5 + 1.1326776948594741e-5, + 1.1343379410666543e-5, + 1.1343379308081936e-5, + 3.679395547040443e-5 ], rtol=2 * sqrt(eps())) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -257,19 +265,19 @@ end volume_integral=VolumeIntegralWeakForm(), surface_integral=SurfaceIntegralWeakForm(FluxHLL(min_max_speed_naive)), l2=[ - 7.905498158659466e-6, - 8.731690809663625e-6, - 8.731690811576996e-6, - 2.9113296018693953e-5 + 7.906577233358824e-6, + 8.733496764180975e-6, + 8.733496764698532e-6, + 2.911852322169076e-5 ], linf=[ - 3.298811230090237e-5, - 4.032272476939269e-5, - 4.032272526011127e-5, - 0.00012013725458537294 + 3.298755256198049e-5, + 4.0322966492922774e-5, + 4.032296598488472e-5, + 0.00012013778942154829 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -281,19 +289,19 @@ end @trixi_testset "elixir_euler_hohqmesh.jl (Quadrilateral elements, SBP, flux differencing)" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_hohqmesh.jl"), l2=[ - 0.0008153911341517156, - 0.0007768159701964676, - 0.00047902606811690694, - 0.0015551846076348535 + 0.0008153911341539523, + 0.0007768159702011952, + 0.0004790260681142826, + 0.0015551846076274918 ], linf=[ - 0.0029301131365355726, - 0.0034427051471457304, - 0.0028721569841545502, - 0.011125365074589944 + 0.002930113136531798, + 0.003442705146861069, + 0.002872156984277563, + 0.011125365075300486 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -313,8 +321,8 @@ end 4.128314378397532, 4.081366752807379 ], rtol = 0.05) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -334,8 +342,8 @@ end 0.0015060064614331736, 0.0019371156800773726, 0.0019371156800769285, 0.004742431684202408 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -354,8 +362,8 @@ end 2.509979394305084e-5, 2.2683711321080935e-5, 2.6180377720841363e-5, 5.575278031910713e-5 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -381,8 +389,8 @@ end 0.1235468309508845, 0.26911424973147735 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -409,8 +417,8 @@ end 0.12344140473946856, 0.26978428189564774 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -424,15 +432,19 @@ end "elixir_euler_rayleigh_taylor_instability.jl"), cells_per_dimension=(8, 8), tspan=(0.0, 0.2), l2=[ - 0.07097806723891838, 0.005168550941966817, - 0.013820912272220933, 0.03243357220022434 + 0.07097806924106471, + 0.005168545523460976, + 0.013820905434253445, + 0.03243358478653133 ], linf=[ - 0.4783395896753895, 0.02244629340135818, - 0.04023357731088538, 0.08515807256615027 + 0.4783395366569936, + 0.022446258588973853, + 0.04023354591166624, + 0.08515791118082117 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -456,8 +468,8 @@ end 0.021957154972646383, 0.33773439650806303 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -470,19 +482,19 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing.jl"), cells_per_dimension=4, tspan=(0.0, 0.1), l2=[ - 0.05685148333985476, - 0.04308122135907089, - 0.043081221359070915, - 0.21098131003847664 + 0.05685180852320552, + 0.04308097439005265, + 0.04308097439005263, + 0.21098250258804 ], linf=[ - 0.2360672306096051, - 0.16684417686971842, - 0.1668441768697189, - 0.8572572782118661 + 0.2360805191601203, + 0.16684117462697776, + 0.16684117462697767, + 0.8573034682049414 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -495,19 +507,19 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_shockcapturing_curved.jl"), cells_per_dimension=4, tspan=(0.0, 0.1), l2=[ - 0.05565849298766252, - 0.042322816017256494, - 0.042322816017256466, - 0.2064212098324083 + 0.055659339125865195, + 0.042323245380073364, + 0.042323245380073315, + 0.20642426004746467 ], linf=[ - 0.23633287875008924, - 0.16930148707515683, - 0.16930148707515688, - 0.8587706761131937 + 0.23633597150568753, + 0.16929779869845438, + 0.16929779869845438, + 0.8587814448153765 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -540,8 +552,8 @@ end 0.0016243096040242655, 0.004447503691245913 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -576,8 +588,8 @@ end 0.003452046522624652, 0.007677153211004928 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -596,8 +608,8 @@ end 2.7000151703315822e-6, 3.988595025372632e-6, 3.9885950240403645e-6, 8.848583036513702e-6 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -622,8 +634,8 @@ end 3.988595024928543e-6, 8.84858303740188e-6 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -654,8 +666,8 @@ end 0.5244468027020814, 1.2210130256735705 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -687,8 +699,8 @@ end 5.067812786885284e-5, 9.887976803746312e-5 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -710,8 +722,8 @@ end 0.7450328339751362, 0.06357382685763713, 0.0635738268576378, 0.1058830287485999, 0.005740591170062146]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -733,8 +745,8 @@ end 0.6906308908961734, 0.05349939593012487, 0.05349939593013042, 0.08830587480616725, 0.0029551359803035027]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -809,8 +821,8 @@ end 0.0, 0.0002992621756946904 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -824,19 +836,19 @@ end cells_per_dimension=8, element_type=Quad(), approximation_type=SBP(), l2=[ - 0.0020316462913319046, - 0.023669019044882247, - 0.03446194752754684, - 1.9333465252381796e-15 + 0.0020316463892983217, + 0.02366902012965938, + 0.03446194535725363, + 1.921676942941478e-15 ], linf=[ - 0.010385010095182778, - 0.08750628939565086, - 0.12088392994348407, + 0.010384996665098178, + 0.08750632767286826, + 0.12088391569555768, 9.325873406851315e-15 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -850,19 +862,19 @@ end cells_per_dimension=8, element_type=Tri(), approximation_type=SBP(), l2=[ - 0.004180680322490383, - 0.07026192411558974, - 0.11815151697006446, - 2.329788936151192e-15 + 0.004180679992535108, + 0.07026193567927695, + 0.11815151184746633, + 2.3786840926019625e-15 ], linf=[ - 0.02076003852980346, - 0.29169601664914424, - 0.5674183379872275, - 1.1546319456101628e-14 + 0.020760033097378283, + 0.29169608872805686, + 0.567418412384793, + 1.1102230246251565e-14 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -878,19 +890,19 @@ end # The last l2, linf error are the L2 projection error in approximating `b`, so they are not # zero for general non-collocated quadrature rules (e.g., for `element_type=Tri()`, `polydeg > 2`). l2=[ - 0.0008309356912456799, - 0.01522451288799231, - 0.016033969387208476, - 1.2820247308150876e-5 + 0.0008309358577296097, + 0.015224511207450263, + 0.016033971785878454, + 1.282024730815488e-5 ], linf=[ - 0.001888045014140971, - 0.05466838692127718, - 0.06345885709961152, - 3.3989933098554914e-5 + 0.0018880416154898327, + 0.05466845626696504, + 0.06345896594568323, + 3.398993309877696e-5 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] @@ -908,19 +920,19 @@ end # a `(polydeg + 1)`-point Gauss quadrature rule in each coordinate (in general, StartUpDG.jl defaults # to the quadrature rule with the fewest number of points which exactly integrates the mass matrix). l2=[ - 7.460461950323111e-5, - 0.003685589808444905, - 0.0039101604749887785, - 2.0636891126652983e-15 + 7.460473151203597e-5, + 0.0036855901000765463, + 0.003910160802530521, + 6.743418333559633e-15 ], linf=[ - 0.000259995400729629, - 0.0072236204211630906, - 0.010364675200833062, - 1.021405182655144e-14 + 0.0002599957400737374, + 0.007223608258381642, + 0.010364657535841815, + 2.042810365310288e-14 ]) - # Ensure that we do not have excessive memory allocations - # (e.g., from type instabilities) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) let t = sol.t[end] u_ode = sol.u[end] diff --git a/test/test_dgmulti_3d.jl b/test/test_dgmulti_3d.jl index 5af8c9ee911..824839f3d26 100644 --- a/test/test_dgmulti_3d.jl +++ b/test/test_dgmulti_3d.jl @@ -230,18 +230,18 @@ end @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_taylor_green_vortex.jl"), polydeg=3, tspan=(0.0, 1.0), cells_per_dimension=(2, 2, 2), l2=[ - 0.0003612827827560599, - 0.06219350883951729, - 0.062193508839503864, - 0.08121963221634831, - 0.07082703570808184 + 0.00036128264902931644, + 0.06219350570157671, + 0.062193505701565316, + 0.08121963725209637, + 0.0708269605813566 ], linf=[ - 0.0007893509649821162, - 0.1481953939988877, - 0.14819539399791176, - 0.14847291108358926, - 0.21313533492212855 + 0.0007893500666786846, + 0.14819541663164099, + 0.14819541663231595, + 0.148472950090691, + 0.2131352319423172 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -258,18 +258,18 @@ end polydeg=3, approximation_type=GaussSBP(), tspan=(0.0, 1.0), cells_per_dimension=(2, 2, 2), l2=[ - 0.00036128278275524326, - 0.062193508839511434, - 0.06219350883949677, - 0.08121963221635205, - 0.07082703570765223 + 0.0003612826490291416, + 0.06219350570157282, + 0.06219350570156088, + 0.08121963725209767, + 0.07082696058040763 ], linf=[ - 0.000789350964946367, - 0.14819539399525805, - 0.14819539399590542, - 0.14847291107658706, - 0.21313533492059378 + 0.0007893500666356079, + 0.14819541663140268, + 0.14819541663222624, + 0.14847295009030398, + 0.21313523192395678 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -354,18 +354,18 @@ end @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"), l2=[ - 7.561896970325353e-5, - 6.884047859361093e-5, - 6.884047859363204e-5, - 6.884047859361148e-5, - 0.000201107274617457 + 7.561468750241556e-5, + 6.882819932057486e-5, + 6.882819932056578e-5, + 6.882819932056221e-5, + 0.0002010869398143227 ], linf=[ - 0.0001337520020225913, - 0.00011571467799287305, - 0.0001157146779990903, - 0.00011571467799376123, - 0.0003446082308800058 + 0.00013375688549710496, + 0.00011568674101658516, + 0.00011568674101614107, + 0.00011568674101658516, + 0.0003446273444156489 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_parabolic_2d.jl b/test/test_parabolic_2d.jl index 0bf69fa85b3..7e80d0e1af9 100644 --- a/test/test_parabolic_2d.jl +++ b/test/test_parabolic_2d.jl @@ -50,22 +50,26 @@ isdir(outdir) && rm(outdir, recursive = true) fill!(gradients[dim], zero(eltype(gradients[dim]))) end + # unpack VectorOfArray + u0 = Base.parent(ode.u0) t = 0.0 # pass in `boundary_condition_periodic` to skip boundary flux/integral evaluation - Trixi.calc_gradient!(gradients, ode.u0, t, mesh, equations_parabolic, + Trixi.calc_gradient!(gradients, u0, t, mesh, equations_parabolic, boundary_condition_periodic, dg, cache, cache_parabolic) @unpack x, y, xq, yq = mesh.md @test getindex.(gradients[1], 1) ≈ 2 * xq .* yq @test getindex.(gradients[2], 1) ≈ xq .^ 2 u_flux = similar.(gradients) - Trixi.calc_viscous_fluxes!(u_flux, ode.u0, gradients, mesh, equations_parabolic, + Trixi.calc_viscous_fluxes!(u_flux, u0, gradients, mesh, + equations_parabolic, dg, cache, cache_parabolic) @test u_flux[1] ≈ gradients[1] @test u_flux[2] ≈ gradients[2] - du = similar(ode.u0) - Trixi.calc_divergence!(du, ode.u0, t, u_flux, mesh, equations_parabolic, + du = similar(u0) + Trixi.calc_divergence!(du, u0, t, u_flux, mesh, + equations_parabolic, boundary_condition_periodic, dg, semi.solver_parabolic, cache, cache_parabolic) @test getindex.(du, 1) ≈ 2 * y @@ -124,16 +128,16 @@ end "elixir_navierstokes_convergence.jl"), cells_per_dimension=(4, 4), tspan=(0.0, 0.1), l2=[ - 0.0015355076812512347, - 0.0033843168272690953, - 0.0036531858107444093, - 0.009948436427520873 + 0.0015355076237431118, + 0.003384316785885901, + 0.0036531858026850757, + 0.009948436101649498 ], linf=[ - 0.005522560467186466, - 0.013425258500736614, - 0.013962115643462503, - 0.0274831021205042 + 0.005522560543588462, + 0.013425258431728926, + 0.013962115936715924, + 0.027483099961148838 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -150,16 +154,16 @@ end "elixir_navierstokes_convergence_curved.jl"), cells_per_dimension=(4, 4), tspan=(0.0, 0.1), l2=[ - 0.004255101916146402, - 0.011118488923215394, - 0.011281831283462784, - 0.035736564473886 + 0.0042551020940351444, + 0.011118489080358264, + 0.011281831362358863, + 0.035736565778376306 ], linf=[ - 0.015071710669707805, - 0.04103132025860057, - 0.03990424085748012, - 0.13094017185987106 + 0.015071709836357083, + 0.04103131887989486, + 0.03990424032494211, + 0.13094018584692968 ],) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -176,16 +180,16 @@ end "elixir_navierstokes_lid_driven_cavity.jl"), cells_per_dimension=(4, 4), tspan=(0.0, 0.5), l2=[ - 0.00022156125227115747, - 0.028318325921401, - 0.009509168701070296, - 0.028267900513550506 + 0.0002215612357465129, + 0.028318325887331217, + 0.009509168805093485, + 0.028267893004691534 ], linf=[ - 0.001562278941298234, - 0.14886653390744856, - 0.0716323565533752, - 0.19472785105241996 + 0.0015622793960574644, + 0.1488665309341318, + 0.07163235778907852, + 0.19472797949052278 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_parabolic_3d.jl b/test/test_parabolic_3d.jl index fc14aa4db31..fe1aa7a0190 100644 --- a/test/test_parabolic_3d.jl +++ b/test/test_parabolic_3d.jl @@ -17,18 +17,18 @@ isdir(outdir) && rm(outdir, recursive = true) "elixir_navierstokes_convergence.jl"), cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1), l2=[ - 0.000553284711585015, - 0.0006592634909666629, - 0.0007776436127373607, - 0.0006592634909664286, - 0.0038073628897819524 + 0.0005532846479614563, + 0.000659263463988067, + 0.0007776436003494915, + 0.000659263463988129, + 0.0038073624941206956 ], linf=[ - 0.0017039861523628907, - 0.002628561703550747, - 0.0035310574250866367, - 0.002628561703585053, - 0.015587829540340437 + 0.001703986341275776, + 0.0026285618026252733, + 0.00353105737957371, + 0.002628561802588858, + 0.015587831432887 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -45,18 +45,18 @@ end "elixir_navierstokes_convergence_curved.jl"), cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.1), l2=[ - 0.001402722725120774, - 0.0021322235533272546, - 0.002787374144745514, - 0.002458747307062751, - 0.009978368180191861 + 0.0014027227340680359, + 0.0021322235583299425, + 0.002787374145873934, + 0.002458747307842109, + 0.009978368214450204 ], linf=[ - 0.0063417504028402405, - 0.010306014252245865, - 0.015207402509253565, - 0.010968264045485343, - 0.04745438983155026 + 0.006341750448945582, + 0.010306014425485621, + 0.015207402553448324, + 0.010968264060799426, + 0.04745438898236998 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) @@ -73,18 +73,18 @@ end "elixir_navierstokes_taylor_green_vortex.jl"), cells_per_dimension=(4, 4, 4), tspan=(0.0, 0.25), l2=[ - 0.0001825713444029892, - 0.015589736382772248, - 0.015589736382771884, - 0.021943924667273653, - 0.01927370280244222 + 0.00018257125088549987, + 0.015589736346235174, + 0.015589736346235415, + 0.021943924698669025, + 0.019273688367502154 ], linf=[ - 0.0006268463584697681, - 0.03218881662749007, - 0.03218881662697948, - 0.053872495395614256, - 0.05183822000984151 + 0.0006268461326666142, + 0.03218881686243058, + 0.03218881686357877, + 0.053872494644958, + 0.05183811394229565 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_threaded.jl b/test/test_threaded.jl index bcc5f21b8ea..1726bcd755f 100644 --- a/test/test_threaded.jl +++ b/test/test_threaded.jl @@ -391,16 +391,16 @@ end "elixir_euler_curved.jl"), alg=RDPK3SpFSAL49(thread = OrdinaryDiffEq.True()), l2=[ - 1.7204593127904542e-5, - 1.5921547179522804e-5, - 1.5921547180107928e-5, - 4.894071422525737e-5 + 1.720916434676505e-5, + 1.5928649356300228e-5, + 1.5928649356913923e-5, + 4.896339454587786e-5 ], linf=[ - 0.00010525416930584619, - 0.00010003778091061122, - 0.00010003778085621029, - 0.00036426282101720275 + 0.00010525404319960963, + 0.00010003768703459315, + 0.00010003768696797977, + 0.0003642622844073351 ]) # Ensure that we do not have excessive memory allocations @@ -469,18 +469,18 @@ end @test_trixi_include(joinpath(examples_dir(), "dgmulti_3d/elixir_euler_fdsbp_periodic.jl"), l2=[ - 7.561896970325353e-5, - 6.884047859361093e-5, - 6.884047859363204e-5, - 6.884047859361148e-5, - 0.000201107274617457 + 7.561468750241556e-5, + 6.882819932057486e-5, + 6.882819932056578e-5, + 6.882819932056221e-5, + 0.0002010869398143227 ], linf=[ - 0.0001337520020225913, - 0.00011571467799287305, - 0.0001157146779990903, - 0.00011571467799376123, - 0.0003446082308800058 + 0.00013375688549710496, + 0.00011568674101658516, + 0.00011568674101614107, + 0.00011568674101658516, + 0.0003446273444156489 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities) diff --git a/test/test_visualization.jl b/test/test_visualization.jl index 0803dae3649..6dcf789a08a 100644 --- a/test/test_visualization.jl +++ b/test/test_visualization.jl @@ -98,7 +98,12 @@ test_examples_2d = Dict("TreeMesh" => ("tree_2d_dgsem", semi = sol.prob.p if mesh == "DGMulti" - scalar_data = StructArrays.component(sol.u[end], 1) + if sol.u[end] isa Trixi.VectorOfArray + u = parent(sol.u[end]) + else + u = sol.u[end] + end + scalar_data = StructArrays.component(u, 1) @test_nowarn_mod Plots.plot(ScalarPlotData2D(scalar_data, semi)) else cache = semi.cache