WIP: Euler with gravity and a total energy that includes the geopotential

examples/p4est_2d_dgsem/elixir_euler_gravity_schaer_mountains.jl

@@ -0,0 +1,135 @@
+using OrdinaryDiffEq
+using Trixi
+
+#schaer mountain test case 
+
+# Initial condition
+function initial_condition_schaer_mountain(x, t,
+                                           equations::CompressibleEulerEquationsWithGravity2D)
+    g = 9.81
+    c_p = 1004.0
+    c_v = 717.0
+    gamma = c_p / c_v
+
+    # Exner pressure from hydrostatic balance for x[2]
+    potential_temperature_int = 280.0 #constant of integration 
+    bvfrequency = 0.01 #Brunt-Väisälä frequency
+
+    exner = 1 +
+            g^2 / (c_p * potential_temperature_int * bvfrequency^2) *
+            (exp(-bvfrequency^2 / g * x[2]) - 1)
+
+    # mean potential temperature
+    potential_temperature = potential_temperature_int * exp(bvfrequency^2 / g * x[2])
+
+    # temperature
+    T = potential_temperature * exner
+
+    # pressure
+    p_0 = 100_000.0  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
+    p = p_0 * exner^(c_p / R)
+
+    # density
+    rho = p / (R * T)
+
+    #Geopotential
+    phi = g * x[2]
+
+    v1 = 10.0
+    v2 = 0.0
+    E = c_v * T + 0.5 * (v1^2 + v2^2) + phi
+    return SVector(rho, rho * v1, rho * v2, rho * E, phi)
+end
+
+###############################################################################
+
+function mapping(xi_, eta_)
+    xi = xi_ * 25_000  #xi_ * 10_000.0 
+    eta = eta_ * 10_500 + 10_500# eta_ * 500.0 + 500.0 
+    #upper boundary 
+    H = 21_000.0
+
+    #topography
+    h_c = 250.0
+    lambda_c = 4000.0
+    a_c = 5000.0
+
+    topo = -h_c * exp(-(xi / a_c)^2) * cos(pi * xi / lambda_c)^2
+
+    x = xi
+    y = H * (eta - topo) / (H - topo)
+    return SVector(x, y)
+end
+
+# Create curved mesh with 200 x 100 elements
+polydeg = 3
+cells_per_dimension = (60, 30)
+mesh = P4estMesh(cells_per_dimension; polydeg = polydeg, mapping = mapping,
+                 periodicity = false)
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquationsWithGravity2D(1004.0 / 717.0)
+
+initial_condition = initial_condition_schaer_mountain
+
+boundary_conditions_dirichlet = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_schaer_mountain),
+                                     :x_pos => BoundaryConditionDirichlet(initial_condition_schaer_mountain),
+                                     :y_neg => boundary_condition_slip_wall,
+                                     :y_pos => boundary_condition_slip_wall)
+
+basis = LobattoLegendreBasis(polydeg)
+
+volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski)
+surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski)
+
+volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
+
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-25_000.0, 0.0)
+coordinates_max = (25_000.0, 21_000.0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions_dirichlet)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 60 * 60)# * 10)  # 10h = 36000 s
+
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+solution_variables = cons2prim
+
+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 = solution_variables)
+
+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/tree_2d_dgsem/elixir_euler_gravity_equilibrium.jl

@@ -0,0 +1,77 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+equations = CompressibleEulerEquationsWithGravity2D(1.4)
+
+function initial_condition_constant(x, t,
+                                    equations::CompressibleEulerEquationsWithGravity2D)
+    rho = exp(-x[2])
+    v1 = 0.0
+    v2 = 0.0
+    p = exp(-x[2])
+    prim = SVector(rho, v1, v2, p, x[2])
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_constant
+
+volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski)
+#surface_flux = (flux_shima_etal, flux_nonconservative_waruszewski)
+#surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski) # bas
+surface_flux = (FluxPlusDissipation(flux_shima_etal,
+                                    DissipationLocalLaxFriedrichs(max_abs_speed_naive)),
+                flux_nonconservative_waruszewski)
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, 0.0)
+coordinates_max = (1.0, 1.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 5,
+                n_cells_max = 10_000,
+                periodicity = false)
+
+boundary_conditions = (; x_neg = boundary_condition_slip_wall,
+                       x_pos = boundary_condition_slip_wall,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.4)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_polydeg = polydeg)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1,
+                                     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);

examples/tree_2d_dgsem/elixir_euler_gravity_warmbubble.jl

@@ -0,0 +1,144 @@
+
+using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
+using Trixi
+
+# Warm bubble test case from
+# - Wicker, L. J., and Skamarock, W. C. (1998)
+#   A time-splitting scheme for the elastic equations incorporating
+#   second-order Runge–Kutta time differencing
+#   [DOI: 10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2](https://doi.org/10.1175/1520-0493(1998)126%3C1992:ATSSFT%3E2.0.CO;2)
+# See also
+# - Bryan and Fritsch (2002)
+#   A Benchmark Simulation for Moist Nonhydrostatic Numerical Models
+#   [DOI: 10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2](https://doi.org/10.1175/1520-0493(2002)130<2917:ABSFMN>2.0.CO;2)
+# - Carpenter, Droegemeier, Woodward, Hane (1990)
+#   Application of the Piecewise Parabolic Method (PPM) to
+#   Meteorological Modeling
+#   [DOI: 10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2](https://doi.org/10.1175/1520-0493(1990)118<0586:AOTPPM>2.0.CO;2)
+struct WarmBubbleSetup
+    # Physical constants
+    g::Float64       # gravity of earth
+    c_p::Float64     # heat capacity for constant pressure (dry air)
+    c_v::Float64     # heat capacity for constant volume (dry air)
+    gamma::Float64   # heat capacity ratio (dry air)
+
+    function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v)
+        new(g, c_p, c_v, gamma)
+    end
+end
+
+# Initial condition
+function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquationsWithGravity2D)
+    RealT = eltype(x)
+    @unpack g, c_p, c_v = setup
+
+    # center of perturbation
+    center_x = 10000
+    center_z = 2000
+    # radius of perturbation
+    radius = 2000
+    # 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
+    potential_temperature_perturbation = zero(RealT)
+    if r <= radius
+        potential_temperature_perturbation = 2 * cospi(0.5f0 * 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) * x[2]
+
+    # pressure
+    p_0 = 100_000  # reference pressure
+    R = c_p - c_v    # gas constant (dry air)
+    p = p_0 * exner^(c_p / R)
+
+    # temperature
+    T = potential_temperature * exner
+    # T = potential_temperature - g / (c_p) * x[2]
+
+    # density
+    rho = p / (R * T)
+
+    # Geopotential
+    phi = g * x[2]
+
+    v1 = 20
+    v2 = 0
+    E = c_v * T + 0.5f0 * (v1^2 + v2^2) + phi
+    return SVector(rho, rho * v1, rho * v2, rho * E, phi)
+end
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+warm_bubble_setup = WarmBubbleSetup()
+
+equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma)
+
+initial_condition = warm_bubble_setup
+
+volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski)
+#surface_flux = (flux_shima_etal, flux_nonconservative_waruszewski)
+# surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski) # bas
+surface_flux = (flux_hllc, flux_nonconservative_waruszewski)
+
+# surface_flux = (FluxPlusDissipation(FluxLMARS(340.0),
+#                                     DissipationLocalLaxFriedrichs(max_abs_speed_naive)),
+#                 flux_nonconservative_waruszewski)
+# volume_flux = (flux_central, flux_nonconservative_waruszewski)
+# surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski)
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (0.0, -5000.0)
+coordinates_max = (20_000.0, 15_000.0)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 6,
+                n_cells_max = 10_000,
+                periodicity = (true, false))
+
+boundary_conditions = (; x_neg = boundary_condition_periodic,
+                       x_pos = boundary_condition_periodic,
+                       y_neg = boundary_condition_slip_wall,
+                       y_pos = boundary_condition_slip_wall)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    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 = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     analysis_polydeg = polydeg)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(dt = 10.0, #interval = 1, #dt = 10.0,
+                                     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);

src/Trixi.jl

@@ -153,6 +153,7 @@ include("visualization/visualization.jl")
 export AcousticPerturbationEquations2D,
        CompressibleEulerEquations1D, CompressibleEulerEquations2D,
        CompressibleEulerEquations3D,
+       CompressibleEulerEquationsWithGravity2D,
        CompressibleEulerMulticomponentEquations1D,
        CompressibleEulerMulticomponentEquations2D,
        CompressibleEulerEquationsQuasi1D,
@@ -188,6 +189,7 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_fjordholm_etal, flux_nonconservative_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,
        flux_chan_etal, flux_nonconservative_chan_etal, flux_winters_etal,
+       flux_nonconservative_waruszewski,
        hydrostatic_reconstruction_audusse_etal, flux_nonconservative_audusse_etal,
        FluxPlusDissipation, DissipationGlobalLaxFriedrichs, DissipationLocalLaxFriedrichs,
        DissipationLaxFriedrichsEntropyVariables, DissipationMatrixWintersEtal,