Skip to content

Commit e5b5534

Browse files
amruedaamrueda
committed
Added benchmark tests
1 parent 78fd1d0 commit e5b5534

File tree

3 files changed

+345
-0
lines changed

3 files changed

+345
-0
lines changed

‎examples/elixir_euler_gravity_equilibrium.jl‎

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
2+
using OrdinaryDiffEq
3+
using Trixi
4+
5+
###############################################################################
6+
# semidiscretization of the compressible Euler equations
7+
equations = CompressibleEulerEquationsWithGravity2D(1.4)
8+
9+
function initial_condition_constant(x, t,
10+
equations::CompressibleEulerEquationsWithGravity2D)
11+
rho = exp(-x[2])
12+
v1 = 0.0
13+
v2 = 0.0
14+
p = exp(-x[2])
15+
prim = SVector(rho, v1, v2, p, x[2])
16+
return prim2cons(prim, equations)
17+
end
18+
19+
initial_condition = initial_condition_constant
20+
21+
volume_flux = (flux_shima_etal, flux_nonconservative_waruszewski)
22+
surface_flux = (flux_lax_friedrichs, flux_nonconservative_waruszewski)
23+
24+
polydeg = 3
25+
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
26+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
27+
28+
coordinates_min = (0.0, 0.0)
29+
coordinates_max = (1.0, 1.0)
30+
mesh = TreeMesh(coordinates_min, coordinates_max,
31+
initial_refinement_level = 5,
32+
n_cells_max = 10_000,
33+
periodicity = false)
34+
35+
boundary_conditions = (; x_neg = boundary_condition_slip_wall,
36+
x_pos = boundary_condition_slip_wall,
37+
y_neg = boundary_condition_slip_wall,
38+
y_pos = boundary_condition_slip_wall)
39+
40+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
41+
boundary_conditions = boundary_conditions)
42+
43+
###############################################################################
44+
# ODE solvers, callbacks etc.
45+
46+
tspan = (0.0, 0.4)
47+
ode = semidiscretize(semi, tspan)
48+
49+
summary_callback = SummaryCallback()
50+
51+
analysis_interval = 1
52+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
53+
analysis_polydeg = polydeg)
54+
55+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
56+
57+
save_solution = SaveSolutionCallback(interval = 1,
58+
save_initial_solution = true,
59+
save_final_solution = true,
60+
solution_variables = cons2prim)
61+
62+
stepsize_callback = StepsizeCallback(cfl = 1.0)
63+
64+
callbacks = CallbackSet(summary_callback,
65+
analysis_callback, alive_callback,
66+
save_solution,
67+
stepsize_callback)
68+
69+
###############################################################################
70+
# run the simulation
71+
72+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
73+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
74+
save_everystep = false, callback = callbacks);

‎examples/elixir_euler_gravity_schaer_mountains.jl‎

Lines changed: 135 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,135 @@
1+
using OrdinaryDiffEq
2+
using Trixi
3+
4+
#schaer mountain test case
5+
6+
# Initial condition
7+
function initial_condition_schaer_mountain(x, t,
8+
equations::CompressibleEulerEquationsWithGravity2D)
9+
g = 9.81
10+
c_p = 1004.0
11+
c_v = 717.0
12+
gamma = c_p / c_v
13+
14+
# Exner pressure from hydrostatic balance for x[2]
15+
potential_temperature_int = 280.0 #constant of integration
16+
bvfrequency = 0.01 #Brunt-Väisälä frequency
17+
18+
exner = 1 +
19+
g^2 / (c_p * potential_temperature_int * bvfrequency^2) *
20+
(exp(-bvfrequency^2 / g * x[2]) - 1)
21+
22+
# mean potential temperature
23+
potential_temperature = potential_temperature_int * exp(bvfrequency^2 / g * x[2])
24+
25+
# temperature
26+
T = potential_temperature * exner
27+
28+
# pressure
29+
p_0 = 100_000.0 # reference pressure
30+
R = c_p - c_v # gas constant (dry air)
31+
p = p_0 * exner^(c_p / R)
32+
33+
# density
34+
rho = p / (R * T)
35+
36+
#Geopotential
37+
phi = g * x[2]
38+
39+
v1 = 10.0
40+
v2 = 0.0
41+
E = c_v * T + 0.5 * (v1^2 + v2^2) + phi
42+
return SVector(rho, rho * v1, rho * v2, rho * E, phi)
43+
end
44+
45+
###############################################################################
46+
47+
function mapping(xi_, eta_)
48+
xi = xi_ * 25_000 #xi_ * 10_000.0
49+
eta = eta_ * 10_500 + 10_500# eta_ * 500.0 + 500.0
50+
#upper boundary
51+
H = 21_000.0
52+
53+
#topography
54+
h_c = 250.0
55+
lambda_c = 4000.0
56+
a_c = 5000.0
57+
58+
topo = -h_c * exp(-(xi / a_c)^2) * cos(pi * xi / lambda_c)^2
59+
60+
x = xi
61+
y = H * (eta - topo) / (H - topo)
62+
return SVector(x, y)
63+
end
64+
65+
# Create curved mesh with 200 x 100 elements
66+
polydeg = 3
67+
cells_per_dimension = (60, 30)
68+
mesh = P4estMesh(cells_per_dimension; polydeg = polydeg, mapping = mapping,
69+
periodicity = false)
70+
71+
###############################################################################
72+
# semidiscretization of the compressible Euler equations
73+
equations = CompressibleEulerEquationsWithGravity2D(1004.0 / 717.0)
74+
75+
initial_condition = initial_condition_schaer_mountain
76+
77+
boundary_conditions_dirichlet = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition_schaer_mountain),
78+
:x_pos => BoundaryConditionDirichlet(initial_condition_schaer_mountain),
79+
:y_neg => boundary_condition_slip_wall,
80+
:y_pos => boundary_condition_slip_wall)
81+
82+
basis = LobattoLegendreBasis(polydeg)
83+
84+
volume_flux = (flux_ranocha, flux_nonconservative_waruszewski)
85+
surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski)
86+
87+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
88+
89+
solver = DGSEM(basis, surface_flux, volume_integral)
90+
91+
coordinates_min = (-25_000.0, 0.0)
92+
coordinates_max = (25_000.0, 21_000.0)
93+
94+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
95+
boundary_conditions = boundary_conditions_dirichlet)
96+
97+
###############################################################################
98+
# ODE solvers, callbacks etc.
99+
100+
tspan = (0.0, 60 * 60)# * 10) # 10h = 36000 s
101+
102+
ode = semidiscretize(semi, tspan)
103+
104+
summary_callback = SummaryCallback()
105+
106+
analysis_interval = 1000
107+
solution_variables = cons2prim
108+
109+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
110+
extra_analysis_errors = (:entropy_conservation_error,))
111+
112+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
113+
114+
save_solution = SaveSolutionCallback(interval = analysis_interval,
115+
save_initial_solution = true,
116+
save_final_solution = true,
117+
output_directory = "out",
118+
solution_variables = solution_variables)
119+
120+
stepsize_callback = StepsizeCallback(cfl = 1.0)
121+
122+
callbacks = CallbackSet(summary_callback,
123+
analysis_callback,
124+
alive_callback,
125+
save_solution,
126+
stepsize_callback)
127+
128+
###############################################################################
129+
# run the simulation
130+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
131+
maxiters = 1.0e7,
132+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
133+
save_everystep = false, callback = callbacks);
134+
135+
summary_callback()

‎examples/elixir_euler_gravity_warmbubble.jl‎

Lines changed: 136 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,136 @@
1+
2+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
3+
using Trixi
4+
5+
# Warm bubble test case from
6+
# - Wicker, L. J., and Skamarock, W. C. (1998)
7+
# A time-splitting scheme for the elastic equations incorporating
8+
# second-order Runge–Kutta time differencing
9+
# [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)
10+
# See also
11+
# - Bryan and Fritsch (2002)
12+
# A Benchmark Simulation for Moist Nonhydrostatic Numerical Models
13+
# [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)
14+
# - Carpenter, Droegemeier, Woodward, Hane (1990)
15+
# Application of the Piecewise Parabolic Method (PPM) to
16+
# Meteorological Modeling
17+
# [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)
18+
struct WarmBubbleSetup
19+
# Physical constants
20+
g::Float64 # gravity of earth
21+
c_p::Float64 # heat capacity for constant pressure (dry air)
22+
c_v::Float64 # heat capacity for constant volume (dry air)
23+
gamma::Float64 # heat capacity ratio (dry air)
24+
25+
function WarmBubbleSetup(; g = 9.81, c_p = 1004.0, c_v = 717.0, gamma = c_p / c_v)
26+
new(g, c_p, c_v, gamma)
27+
end
28+
end
29+
30+
# Initial condition
31+
function (setup::WarmBubbleSetup)(x, t, equations::CompressibleEulerEquationsWithGravity2D)
32+
RealT = eltype(x)
33+
@unpack g, c_p, c_v = setup
34+
35+
# center of perturbation
36+
center_x = 10000
37+
center_z = 2000
38+
# radius of perturbation
39+
radius = 2000
40+
# distance of current x to center of perturbation
41+
r = sqrt((x[1] - center_x)^2 + (x[2] - center_z)^2)
42+
43+
# perturbation in potential temperature
44+
potential_temperature_ref = 300
45+
potential_temperature_perturbation = zero(RealT)
46+
if r <= radius
47+
potential_temperature_perturbation = 2 * cospi(0.5f0 * r / radius)^2
48+
end
49+
potential_temperature = potential_temperature_ref + potential_temperature_perturbation
50+
51+
# Exner pressure, solves hydrostatic equation for x[2]
52+
exner = 1 - g / (c_p * potential_temperature) * x[2]
53+
54+
# pressure
55+
p_0 = 100_000 # reference pressure
56+
R = c_p - c_v # gas constant (dry air)
57+
p = p_0 * exner^(c_p / R)
58+
59+
# temperature
60+
T = potential_temperature * exner
61+
# T = potential_temperature - g / (c_p) * x[2]
62+
63+
# density
64+
rho = p / (R * T)
65+
66+
# Geopotential
67+
phi = g * x[2]
68+
69+
v1 = 20
70+
v2 = 0
71+
E = c_v * T + 0.5f0 * (v1^2 + v2^2) + phi
72+
return SVector(rho, rho * v1, rho * v2, rho * E, phi)
73+
end
74+
75+
###############################################################################
76+
# semidiscretization of the compressible Euler equations
77+
warm_bubble_setup = WarmBubbleSetup()
78+
79+
equations = CompressibleEulerEquationsWithGravity2D(warm_bubble_setup.gamma)
80+
81+
initial_condition = warm_bubble_setup
82+
83+
volume_flux = (flux_kennedy_gruber, flux_nonconservative_waruszewski)
84+
surface_flux = (FluxLMARS(340.0), flux_nonconservative_waruszewski)
85+
86+
polydeg = 3
87+
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
88+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
89+
90+
coordinates_min = (0.0, -5000.0)
91+
coordinates_max = (20_000.0, 15_000.0)
92+
mesh = TreeMesh(coordinates_min, coordinates_max,
93+
initial_refinement_level = 6,
94+
n_cells_max = 10_000,
95+
periodicity = (true, false))
96+
97+
boundary_conditions = (; x_neg = boundary_condition_periodic,
98+
x_pos = boundary_condition_periodic,
99+
y_neg = boundary_condition_slip_wall,
100+
y_pos = boundary_condition_slip_wall)
101+
102+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
103+
boundary_conditions = boundary_conditions)
104+
105+
###############################################################################
106+
# ODE solvers, callbacks etc.
107+
108+
tspan = (0.0, 1000.0) # 1000 seconds final time
109+
ode = semidiscretize(semi, tspan)
110+
111+
summary_callback = SummaryCallback()
112+
113+
analysis_interval = 100
114+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
115+
analysis_polydeg = polydeg)
116+
117+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
118+
119+
save_solution = SaveSolutionCallback(dt = 10.0, #interval = 1, #dt = 10.0,
120+
save_initial_solution = true,
121+
save_final_solution = true,
122+
solution_variables = cons2prim)
123+
124+
stepsize_callback = StepsizeCallback(cfl = 1.0)
125+
126+
callbacks = CallbackSet(summary_callback,
127+
analysis_callback, alive_callback,
128+
save_solution,
129+
stepsize_callback)
130+
131+
###############################################################################
132+
# run the simulation
133+
134+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
135+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
136+
save_everystep = false, callback = callbacks);

0 commit comments

Comments
 (0)