gravity wave test, testing ec fluxes

trixi-framework · Sep 25, 2024 · 14c3820 · 14c3820
1 parent a820c97
commit 14c3820
Show file tree

Hide file tree

Showing 9 changed files with 218 additions and 8 deletions.
diff --git a/examples/elixir_euler_potential_temperature_dry_bubble.jl b/examples/elixir_euler_potential_temperature_dry_bubble.jl
@@ -1,7 +1,6 @@
 using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
-using TrixiAtmo: FluxLMARS, source_terms_gravity, flux_theta
 
 function initial_condition_warm_bubble(x, t,
                                        equations::CompressibleEulerPotentialTemperatureEquations2D)

diff --git a/examples/elixir_euler_potential_temperature_ec.jl b/examples/elixir_euler_potential_temperature_ec.jl
@@ -1,7 +1,7 @@
 using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
-using TrixiAtmo: flux_theta
+
 ###############################################################################
 # semidiscretization of the compressible Euler equations
 equations = CompressibleEulerPotentialTemperatureEquations2D()

diff --git a/examples/elixir_euler_potential_temperature_gravity_wave.jl b/examples/elixir_euler_potential_temperature_gravity_wave.jl
@@ -0,0 +1,102 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+
+function initial_condition_gravity_wave(x, t,
+                                        equations::CompressibleEulerPotentialTemperatureEquations2D)
+    g = equations.g
+    c_p = equations.c_p
+    c_v = equations.c_v
+    # center of perturbation
+    x_c = 100_000.0
+    a = 5_000
+    H = 10_000
+    # perturbation in potential temperature
+    potential_temperature_ref = 300.0 * exp(0.01^2 / g * x[2])
+    potential_temperature_perturbation = 0.01 * sinpi(x[2] / H) / (1 + (x[1] - x_c)^2 / a^2)
+    potential_temperature = potential_temperature_ref + potential_temperature_perturbation
+
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = 1 + g^2 / (c_p * 300.0 * 0.01^2) * (exp(-0.01^2 / g * x[2]) - 1)
+
+    # pressure
+    p_0 = 100_000.0  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
+    p = p_0 * exner^(c_p / R)
+    T = potential_temperature * exner
+
+    # density
+    rho = p / (R * T)
+    v1 = 20.0
+    v2 = 0.0
+
+    return SVector(rho, rho * v1, rho * v2, rho * potential_temperature)
+end
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerPotentialTemperatureEquations2D()
+
+boundary_conditions = (x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
+
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+
+surface_flux = FluxLMARS(340.0)
+
+solver = DGSEM(basis, surface_flux)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (300_000.0, 10_000.0)
+
+cells_per_dimension = (600, 20) # Delta x = Delta z = 1 km
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = (true, false))
+initial_condition = initial_condition_gravity_wave
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_gravity,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3000.0)  # 1000 seconds final time
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:entropy_conservation_error,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = analysis_interval,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     output_directory = "out",
+                                     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),
+            maxiters = 1.0e7,
+            dt = 1e-1, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback()
diff --git a/examples/elixir_moist_euler_dry_bubble.jl b/examples/elixir_moist_euler_dry_bubble.jl
@@ -1,7 +1,7 @@
 using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
-using TrixiAtmo: FluxLMARS, source_terms_geopotential, cons2drypot
+using TrixiAtmo: source_terms_geopotential, cons2drypot
 
 ###############################################################################
 # semidiscretization of the compressible moist Euler equations

diff --git a/examples/elixir_moist_euler_moist_bubble.jl b/examples/elixir_moist_euler_moist_bubble.jl
@@ -1,7 +1,6 @@
 using OrdinaryDiffEq
 using Trixi, TrixiAtmo
-using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble,
-                 FluxLMARS
+using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble
 using NLsolve: nlsolve
 
 ###############################################################################

diff --git a/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl b/examples/elixir_moist_euler_nonhydrostatic_gravity_waves.jl
@@ -3,7 +3,7 @@ using Trixi, TrixiAtmo
 using TrixiAtmo: CompressibleMoistEulerEquations2D, source_terms_geopotential,
                  source_terms_phase_change,
                  source_terms_nonhydrostatic_rayleigh_sponge,
-                 cons2drypot, FluxLMARS
+                 cons2drypot
 
 ###############################################################################
 # semidiscretization of the compressible moist Euler equation

diff --git a/src/TrixiAtmo.jl b/src/TrixiAtmo.jl
@@ -26,7 +26,8 @@ include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
 export CompressibleMoistEulerEquations2D
 export CompressibleEulerPotentialTemperatureEquations2D
 
-export flux_chandrashekar, FluxLMARS, flux_theta
+export flux_chandrashekar, flux_theta, flux_theta_rhoAM, flux_theta_physentropy,
+       flux_theta_physentropy2
 
 export examples_dir
 

diff --git a/src/equations/compressible_euler_potential_temperature_2d.jl b/src/equations/compressible_euler_potential_temperature_2d.jl
@@ -225,6 +225,88 @@ end
     return SVector(f1, f2, f3, f4)
 end
 
+@inline function flux_theta_rhoAM(u_ll, u_rr, normal_direction::AbstractVector,
+                                  equations::CompressibleEulerPotentialTemperatureEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+    _, _, _, rho_theta_ll = u_ll
+    _, _, _, rho_theta_rr = u_rr
+    # Compute the necessary mean values
+    rho_mean = 0.5f0 * (rho_ll + rho_rr)
+
+    gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma)
+
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+    f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    return SVector(f1, f2, f3, f4)
+end
+
+@inline function flux_theta_physentropy(u_ll, u_rr, normal_direction::AbstractVector,
+                                        equations::CompressibleEulerPotentialTemperatureEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+    _, _, _, rho_theta_ll = u_ll
+    _, _, _, rho_theta_rr = u_rr
+    # Compute the necessary mean values
+    #rho_mean = ln_mean(rho_ll, rho_rr)
+    #rho_mean = 0.5f0 * (rho_ll + rho_rr)
+    rho_mean = ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr)
+    gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma)
+
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
+
+    # Calculate fluxes depending on normal_direction
+    #f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    f4 = gammamean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    f1 = f4 * rho_mean
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+
+    return SVector(f1, f2, f3, f4)
+end
+
+@inline function flux_theta_physentropy2(u_ll, u_rr, normal_direction::AbstractVector,
+                                         equations::CompressibleEulerPotentialTemperatureEquations2D)
+    # Unpack left and right state
+    rho_ll, v1_ll, v2_ll, p_ll = cons2prim(u_ll, equations)
+    rho_rr, v1_rr, v2_rr, p_rr = cons2prim(u_rr, equations)
+    v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2]
+    v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2]
+    _, _, _, rho_theta_ll = u_ll
+    _, _, _, rho_theta_rr = u_rr
+    # Compute the necessary mean values
+    rho_mean = ln_mean(rho_ll, rho_rr)
+    #rho_mean = 0.5f0 * (rho_ll + rho_rr)
+    #rho_mean = ln_mean(rho_ll/rho_theta_ll, rho_rr/rho_theta_rr)
+    gammamean = stolarsky_mean(rho_theta_ll, rho_theta_rr, equations.gamma)
+
+    v1_avg = 0.5f0 * (v1_ll + v1_rr)
+    v2_avg = 0.5f0 * (v2_ll + v2_rr)
+    p_avg = 0.5f0 * (p_ll + p_rr)
+
+    # Calculate fluxes depending on normal_direction
+    f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
+    f2 = f1 * v1_avg + p_avg * normal_direction[1]
+    f3 = f1 * v2_avg + p_avg * normal_direction[2]
+    f4 = inv_ln_mean(rho_ll / rho_theta_ll, rho_rr / rho_theta_rr) * f1
+    return SVector(f1, f2, f3, f4)
+end
+
 @inline function prim2cons(prim,
                            equations::CompressibleEulerPotentialTemperatureEquations2D)
     rho, v1, v2, p = prim
@@ -267,7 +349,7 @@ end
     # Mathematical entropy
     p = equations.p_0 * (equations.R * cons[4] / equations.p_0)^equations.gamma
 
-    U = (p / (equations.gamma - 1) + 1 / 2 * (cons[2]^2 + cons[3]^2) / (cons[1]))
+    U = (p / (equations.gamma - 1) + 0.5f0 * (cons[2]^2 + cons[3]^2) / (cons[1]))
 
     return U
 end

diff --git a/test/test_2d_euler_potential_temperature.jl b/test/test_2d_euler_potential_temperature.jl
@@ -34,6 +34,33 @@ EXAMPLES_DIR = pkgdir(TrixiAtmo, "examples") # TODO - Do we need a subdirectory
     end
 end
 
+@trixiatmo_testset "elixir_euler_potential_temperature_gravity_wave" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_euler_potential_temperature_gravity_wave.jl"),
+                        l2=[
+                            8.92434879991671e-9,
+                            1.8131676850378482e-7,
+                            2.6374650049135436e-5,
+                            1.4620631924879524e-6,
+                        ],
+                        linf=[
+                            7.544232794032268e-8,
+                            1.654911628179434e-6,
+                            0.00023724858495745153,
+                            1.8884994915424613e-5,
+                        ],
+                        polydeg=3,
+                        tspan=(0.0, 1.0))
+    # 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
+
 @trixiatmo_testset "elixir_euler_potential_temperature_ec" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_potential_temperature_ec.jl"),