TrixiAtmo.jl/examples/elixir_euler_potential_temperature_inertia_gravity_waves.jl at 0011d93ce3f552ae2e7e56ae1bb50b555eee6158 · trixi-framework/TrixiAtmo.jl · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
using OrdinaryDiffEqSSPRK
using Trixi, TrixiAtmo

"""
	initial_condition_gravity_waves(x, t,
                                        equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D)

Test cases for linearized analytical solution by
-  Baldauf, Michael and Brdar, Slavko (2013)
   An analytic solution for linear gravity waves in a channel as a test
   for numerical models using the non-hydrostatic, compressible {E}uler equations
   [DOI: 10.1002/qj.2105] (https://doi.org/10.1002/qj.2105)
"""
function initial_condition_gravity_waves(x, t,
                                         equations::CompressibleEulerPotentialTemperatureEquationsWithGravity2D)
    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
    R = c_p - c_v    # gas constant (dry air)
    T0 = 250
    delta = g / (R * T0)
    DeltaT = 0.001
    Tb = DeltaT * sinpi(x[2] / H) * exp(-(x[1] - x_c)^2 / a^2)
    ps = 100_000  # reference pressure
    rhos = ps / (T0 * R)
    rho_b = rhos * (-Tb / T0)
    p = ps * exp(-delta * x[2])
    rho = rhos * exp(-delta * x[2]) + rho_b * exp(-0.5 * delta * x[2])
    v1 = 20
    v2 = 0

    return prim2cons(SVector(rho, v1, v2, p, g * x[2]), equations)
end

equations = CompressibleEulerPotentialTemperatureEquationsWithGravity2D(c_p = 1004,
                                                                        c_v = 717,
                                                                        gravity = EARTH_GRAVITATIONAL_ACCELERATION)

# We have an isothermal background state with T0 = 250 K.
# The reference speed of sound can be computed as:
# cs = sqrt(gamma * R * T0)
cs = sqrt(equations.gamma * equations.R * 250)
surface_flux = (FluxLMARS(cs), flux_zero)
volume_flux = (flux_tec, flux_nonconservative_waruszewski_etal)
polydeg = 3
solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

boundary_conditions = (x_neg = boundary_condition_periodic,
                       x_pos = boundary_condition_periodic,
                       y_neg = boundary_condition_slip_wall,
                       y_pos = boundary_condition_slip_wall)

coordinates_min = (0.0, 0.0)
coordinates_max = (300_000.0, 10_000.0)
cells_per_dimension = (60, 8)
mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max,
                      periodicity = (true, false))
source_terms = nothing
initial_condition = initial_condition_gravity_waves
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms = source_terms,
                                    boundary_conditions = boundary_conditions)
tspan = (0.0, 1800.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 10000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                     extra_analysis_integrals = (entropy,))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

stepsize_callback = StepsizeCallback(cfl = 1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        stepsize_callback)

sol = solve(ode,
            SSPRK43(thread = Trixi.True());
            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, adaptive = false)