From e382140be0f9a5dcd4b9fd072c9b420f7793771d Mon Sep 17 00:00:00 2001 From: Marco Artiano Date: Mon, 23 Sep 2024 15:40:49 +0200 Subject: [PATCH] added ec test --- .../elixir_euler_potential_temperature_ec.jl | 77 +++++++++++++++++++ src/TrixiAtmo.jl | 2 +- test/test_2d_euler_potential_temperature.jl | 29 ++++++- test/test_2d_moist_euler.jl | 28 +++---- test/test_spherical_advection.jl | 4 +- 5 files changed, 121 insertions(+), 19 deletions(-) create mode 100644 examples/elixir_euler_potential_temperature_ec.jl diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl new file mode 100644 index 0000000..eb27c5e --- /dev/null +++ b/examples/elixir_euler_potential_temperature_ec.jl @@ -0,0 +1,77 @@ +using OrdinaryDiffEq +using Trixi +using TrixiAtmo +using TrixiAtmo: flux_theta, prim2cons +############################################################################### +# semidiscretization of the compressible Euler equations +equations = CompressibleEulerPotentialTemperatureEquations2D() + +function initial_condition_weak_blast_wave(x, t, + equations::CompressibleEulerPotentialTemperatureEquations2D) + # From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) + # Set up polar coordinates + inicenter = SVector(0, 0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + RealT = eltype(x) + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos_phi + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin_phi + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) + + return prim2cons(SVector(rho, v1, v2, p), equations) +end + +initial_condition = initial_condition_weak_blast_wave + +volume_flux = flux_theta +solver = DGSEM(polydeg = 3, surface_flux = flux_theta, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) + +cells_per_dimension = (32, 32) +mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max, + periodicity = (true, true)) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + boundary_conditions = boundary_condition_periodic) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 1.0) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl index cdce0a3..da2cc5b 100644 --- a/src/TrixiAtmo.jl +++ b/src/TrixiAtmo.jl @@ -26,7 +26,7 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl") export CompressibleMoistEulerEquations2D export CompressibleEulerPotentialTemperatureEquations2D -export flux_chandrashekar, flux_LMARS, flux_theta +export flux_chandrashekar, flux_LMARS, flux_theta, prim2cons export examples_dir diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl index 0b96e69..69629af 100644 --- a/test/test_2d_euler_potential_temperature.jl +++ b/test/test_2d_euler_potential_temperature.jl @@ -14,13 +14,13 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 1.300427456987984e-6, 2.6010873806861595e-5, 0.0006660120005093007, - 9.51074191163579e-6, + 9.51074191163579e-6 ], linf=[ 1.031131432183141e-5, 0.00020623855042956052, 0.006392070867001616, - 7.841424786647622e-5, + 7.841424786647622e-5 ], polydeg=3, tspan=(0.0, 0.1)) @@ -34,4 +34,29 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory end end +@trixiatmo_testset "elixir_euler_potential_temperature_ec" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, + "elixir_euler_potential_temperature_ec.jl"), + l2=[ + 0.06174365254988615, + 0.05018008812040643, + 0.05018774429827882, + 0.0057820937341387935 + ], + linf=[ + 0.2932352942992503, + 0.3108954213591686, + 0.31082193857647905, + 0.027465217019157606 + ]) + # 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 diff --git a/test/test_2d_moist_euler.jl b/test/test_2d_moist_euler.jl index d669e96..646ea0e 100644 --- a/test/test_2d_moist_euler.jl +++ b/test/test_2d_moist_euler.jl @@ -16,7 +16,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.0006660124630171347, 0.008969786054960861, 0.0, - 0.0, + 0.0 ], linf=[ 1.0312042909910168e-5, @@ -24,7 +24,7 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory 0.006392107590872236, 0.07612038028310053, 0.0, - 0.0, + 0.0 ], polydeg=3, tspan=(0.0, 0.1)) @@ -46,7 +46,7 @@ end 7.938812668709457, 4500.747616411578, 0.00015592413050814787, - 0.00014163475049532796, + 0.00014163475049532796 ], linf=[ 0.1427479052298466, @@ -54,7 +54,7 @@ end 91.56822550162855, 49528.383866247605, 0.0019364397602254623, - 0.0013259689889851285, + 0.0013259689889851285 ], polydeg=3, cells_per_dimension=(16, 16), @@ -77,7 +77,7 @@ end 0.0006974588377288118, 1.715668353329522, 8.831269143134121e-7, - 1.025579538944668e-6, + 1.025579538944668e-6 ], linf=[ 8.055695643149896e-5, @@ -85,7 +85,7 @@ end 0.005897639251024437, 19.24776030163048, 1.0043133039065386e-5, - 1.1439046776775402e-5, + 1.1439046776775402e-5 ], polydeg=3, cells_per_dimension=(16, 8), @@ -109,7 +109,7 @@ end 0.01172830034500581, 9.898560584459009, 0.0, - 0.0, + 0.0 ], linf=[ 0.0017602202683439927, @@ -117,7 +117,7 @@ end 0.5945327351674782, 489.89171406268724, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -140,7 +140,7 @@ end 2.609465609739115e-5, 6.323484379066432e-5, 0.0, - 0.0, + 0.0 ], linf=[ 7.549984224430872e-5, @@ -148,7 +148,7 @@ end 0.00015964938767742964, 0.0005425860570698049, 0.0, - 0.0, + 0.0 ], polydeg=3, cells_per_dimension=(10, 8), @@ -171,7 +171,7 @@ end 0.015104690877128945, 0.5242098451814421, 5.474006801215573e-10, - 1.1103883907226752e-10, + 1.1103883907226752e-10 ], linf=[ 0.00013219484616722177, @@ -179,7 +179,7 @@ end 0.03789645369775574, 3.90888311638264, 3.938382289041286e-9, - 6.892033377287209e-10, + 6.892033377287209e-10 ], polydeg=3, cells_per_dimension=(10, 8), @@ -203,7 +203,7 @@ end 0.07460345535073129, 5.943049264963717, 4.471792794168725e-9, - 7.10320253652373e-10, + 7.10320253652373e-10 ], linf=[ 0.0007084183215528839, @@ -211,7 +211,7 @@ end 0.3697160082709827, 27.843155286573165, 2.1168438904322837e-8, - 3.691699932047233e-9, + 3.691699932047233e-9 ], polydeg=3, cells_per_dimension=(10, 8), diff --git a/test/test_spherical_advection.jl b/test/test_spherical_advection.jl index 31f7a29..c826091 100644 --- a/test/test_spherical_advection.jl +++ b/test/test_spherical_advection.jl @@ -15,14 +15,14 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") 8.23414117e-03, 1.84210648e-03, 0.00000000e+00, - 6.44302430e-02, + 6.44302430e-02 ], linf=[ 9.48950488e-02, 9.64811952e-02, 1.37453400e-02, 0.00000000e+00, - 4.09322999e-01, + 4.09322999e-01 ]) # Ensure that we do not have excessive memory allocations # (e.g., from type instabilities)