From 15062dcd11d13e712220c1ccbd1a8384f890e02a Mon Sep 17 00:00:00 2001 From: Tristan Montoya Date: Tue, 3 Sep 2024 10:17:15 +0200 Subject: [PATCH] Add test for spherical advection in Cartesian coordinates (#30) * add tests for spherical advection in Cartesian coords * revert change to moist bubble case * revert change to moist bubble case * hard-code dimension of node coordinates to avoid type instability * metric terms for shell of radius != 1 --- ...xir_euler_spherical_advection_cartesian.jl | 2 + .../containers_2d_manifold_in_3d.jl | 19 +++++++--- test/runtests.jl | 4 ++ test/test_spherical_advection.jl | 37 +++++++++++++++++++ 4 files changed, 56 insertions(+), 6 deletions(-) create mode 100644 test/test_spherical_advection.jl diff --git a/examples/elixir_euler_spherical_advection_cartesian.jl b/examples/elixir_euler_spherical_advection_cartesian.jl index 4a3bb4c..9a0e7b5 100644 --- a/examples/elixir_euler_spherical_advection_cartesian.jl +++ b/examples/elixir_euler_spherical_advection_cartesian.jl @@ -6,6 +6,8 @@ using TrixiAtmo ############################################################################### # semidiscretization of the linear advection equation +# We use the Euler equations structure but modify the rhs! function to convert it to a +# variable-coefficient advection equation equations = CompressibleEulerEquations3D(1.4) # Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux diff --git a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl index 538b895..1fb976b 100644 --- a/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl +++ b/src/solvers/dgsem_p4est/containers_2d_manifold_in_3d.jl @@ -150,7 +150,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates, # Hanging nodes will cause holes in the mesh if its polydeg is higher # than the polydeg of the solver. @assert length(basis.nodes)>=length(mesh.nodes) "The solver can't have a lower polydeg than the mesh" - + @assert size(mesh.tree_node_coordinates, 1)==3 "Shell must use 3D node coordinates" calc_node_coordinates_2d_shell!(node_coordinates, mesh, basis.nodes) end @@ -162,7 +162,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates, # places and the additional information passed to the compiler makes them faster # than native `Array`s. tmp1 = Trixi.StrideArray(undef, real(mesh), - Trixi.StaticInt(size(mesh.tree_node_coordinates, 1)), + Trixi.StaticInt(3), Trixi.static_length(nodes), Trixi.static_length(mesh.nodes)) matrix1 = Trixi.StrideArray(undef, real(mesh), @@ -173,7 +173,6 @@ function calc_node_coordinates_2d_shell!(node_coordinates, # Macros from `p4est` p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL - p4est_root_len = 1 << Trixi.P4EST_MAXLEVEL p4est_quadrant_len(l) = 1 << (Trixi.P4EST_MAXLEVEL - l) trees = Trixi.unsafe_wrap_sc(Trixi.p4est_tree_t, mesh.p4est.trees) @@ -208,6 +207,7 @@ function calc_node_coordinates_2d_shell!(node_coordinates, return node_coordinates end +# This only works for a sphere function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArray{<:Any, 5}, element, @@ -216,6 +216,10 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr @unpack derivative_matrix = basis for j in eachnode(basis), i in eachnode(basis) + radius = sqrt(node_coordinates[1, i, j, element]^2 + + node_coordinates[2, i, j, element]^2 + + node_coordinates[3, i, j, element]^2) + for n in 1:3 # (n, m, l) cyclic m = (n % 3) + 1 @@ -229,7 +233,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr jacobian_matrix[l, 2, i, j, element] * node_coordinates[m, i, j, - element]) + element]) / + radius contravariant_vectors[n, 2, i, j, element] = (jacobian_matrix[l, 1, i, j, element] * @@ -239,7 +244,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr jacobian_matrix[m, 1, i, j, element] * node_coordinates[l, i, j, - element]) + element]) / + radius contravariant_vectors[n, 3, i, j, element] = (jacobian_matrix[m, 1, i, j, element] * @@ -249,7 +255,8 @@ function calc_contravariant_vectors_2d_shell!(contravariant_vectors::AbstractArr jacobian_matrix[m, 2, i, j, element] * jacobian_matrix[l, 1, i, j, - element]) + element]) / + radius end end diff --git a/test/runtests.jl b/test/runtests.jl index c59a4fa..541cfd8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,4 +16,8 @@ const TRIXIATMO_NTHREADS = clamp(Sys.CPU_THREADS, 2, 3) @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "moist_euler" include("test_2d_moist_euler.jl") end + + @time if TRIXIATMO_TEST == "all" || TRIXIATMO_TEST == "spherical_advection" + include("test_spherical_advection.jl") + end end diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl new file mode 100644 index 0000000..31f7a29 --- /dev/null +++ b/test/test_spherical_advection.jl @@ -0,0 +1,37 @@ +module TestSphericalAdvection + +using Test +using TrixiAtmo + +include("test_trixiatmo.jl") + +EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") + +@trixiatmo_testset "elixir_euler_spherical_advection_cartesian" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_spherical_advection_cartesian.jl"), + l2=[ + 8.44505073e-03, + 8.23414117e-03, + 1.84210648e-03, + 0.00000000e+00, + 6.44302430e-02, + ], + linf=[ + 9.48950488e-02, + 9.64811952e-02, + 1.37453400e-02, + 0.00000000e+00, + 4.09322999e-01, + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated TrixiAtmo.Trixi.rhs!(du_ode, u_ode, semi, t)) < 100 + end +end + +end # module