Add 2D Potential Temperature Formulation for the Euler Eqs.

examples/elixir_euler_potential_temperature_dry_bubble.jl

@@ -0,0 +1,112 @@
+using OrdinaryDiffEq
+using Trixi, TrixiAtmo
+using TrixiAtmo: source_terms_gravity
+
+function initial_condition_warm_bubble(x, t,
+                                       equations::CompressibleEulerPotentialTemperatureEquations2D)
+    g = equations.g
+    c_p = equations.c_p
+    c_v = equations.c_v
+    # center of perturbation
+    center_x = 10000.0
+    center_z = 2000.0
+    # radius of perturbation
+    radius = 2000.0
+    # distance of current x to center of perturbation
+    r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
+
+    # perturbation in potential temperature
+    potential_temperature_ref = 300.0
+    potential_temperature_perturbation = 0.0
+    if r <= radius
+        potential_temperature_perturbation = 2 * cospi(0.5 * r / radius)^2
+    end
+    potential_temperature = potential_temperature_ref + potential_temperature_perturbation
+
+    # Exner pressure, solves hydrostatic equation for x[2]
+    exner = 1 - g / (c_p * potential_temperature_ref) * x[2]
+
+    # 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)
+
+volume_flux = flux_theta
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (20_000.0, 10_000.0)
+
+cells_per_dimension = (64, 32)
+mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
+                      periodicity = (true, false))
+initial_condition = initial_condition_warm_bubble
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms_gravity,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1000.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 = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);
+
+summary_callback()

examples/elixir_euler_potential_temperature_ec.jl

@@ -0,0 +1,77 @@
+using OrdinaryDiffEq
+using Trixi
+using TrixiAtmo
+
+###############################################################################
+# 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 TrixiAtmo.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

examples/elixir_euler_potential_temperature_gravity_wave.jl

@@ -0,0 +1,102 @@
+using OrdinaryDiffEq
+using Trixi, TrixiAtmo
+using TrixiAtmo: source_terms_gravity
+
+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()

examples/elixir_moist_euler_dry_bubble.jl

@@ -1,7 +1,7 @@
 using OrdinaryDiffEq
 using Trixi
 using TrixiAtmo
-using TrixiAtmo: flux_LMARS, source_terms_geopotential, cons2drypot
+using TrixiAtmo: source_terms_geopotential, cons2drypot
 
 ###############################################################################
 # semidiscretization of the compressible moist Euler equations
@@ -61,7 +61,7 @@ source_term = source_terms_geopotential
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
 
-surface_flux = flux_LMARS
+surface_flux = FluxLMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

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,
-                 flux_LMARS
+using TrixiAtmo: cons2aeqpot, saturation_pressure, source_terms_moist_bubble
 using NLsolve: nlsolve
 
 ###############################################################################
@@ -242,7 +241,7 @@ source_term = source_terms_moist_bubble
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
 
-surface_flux = flux_LMARS
+surface_flux = FluxLMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

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, flux_LMARS
+                 cons2drypot
 
 ###############################################################################
 # semidiscretization of the compressible moist Euler equation
@@ -55,7 +55,7 @@ boundary_conditions = (x_neg = boundary_condition_periodic,
 
 polydeg = 4
 basis = LobattoLegendreBasis(polydeg)
-surface_flux = flux_LMARS
+surface_flux = FluxLMARS(360.0)
 volume_flux = flux_chandrashekar
 
 volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

src/TrixiAtmo.jl

@@ -24,8 +24,10 @@ include("solvers/solvers.jl")
 include("semidiscretization/semidiscretization_hyperbolic_2d_manifold_in_3d.jl")
 
 export CompressibleMoistEulerEquations2D
+export CompressibleEulerPotentialTemperatureEquations2D
 
-export flux_chandrashekar, flux_LMARS
+export flux_chandrashekar, flux_theta, flux_theta_rhoAM, flux_theta_physentropy,
+       flux_theta_physentropy2
 
 export examples_dir