TrixiAtmo.jl/examples/elixir_moist_euler_dry_bubble.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
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
using OrdinaryDiffEqLowStorageRK
using Trixi
using TrixiAtmo
using TrixiAtmo: source_terms_geopotential, cons2drypot

###############################################################################
# semidiscretization of the compressible moist Euler equations
c_pd = 1004 # specific heat at constant pressure for dry air
c_vd = 717  # specific heat at constant volume for dry air
c_pv = 1885 # specific heat at constant pressure for moist air
c_vv = 1424 # specific heat at constant volume for moist air
equations = CompressibleMoistEulerEquations2D(c_pd = c_pd, c_vd = c_vd, c_pv = c_pv,
                                              c_vv = c_vv,
                                              gravity = EARTH_GRAVITATIONAL_ACCELERATION)

# Warm bubble test from paper:
# Wicker, L. J., and W. C. Skamarock, 1998: A time-splitting scheme
# for the elastic equations incorporating second-order Runge–Kutta
# time differencing. Mon. Wea. Rev., 126, 1992–1999.
function initial_condition_warm_bubble(x, t, equations::CompressibleMoistEulerEquations2D)
    @unpack p_0, kappa, g, c_pd, c_vd, R_d, R_v = equations
    xc = 10000
    zc = 2000
    r = sqrt((x[1] - xc)^2 + (x[2] - zc)^2)
    rc = 2000
    θ_ref = 300
    Δθ = 0

    if r <= rc
        Δθ = 2 * cospi(0.5f0 * r / rc)^2
    end

    # Perturbed state:
    θ = θ_ref + Δθ # potential temperature
    # π_exner = 1 - g / (c_pd * θ) * x[2] # exner pressure
    # rho = p_0 / (R_d * θ) * (π_exner)^(c_vd / R_d) # density

    # calculate background pressure with assumption hydrostatic and neutral
    p = p_0 * (1 - kappa * g * x[2] / (R_d * θ_ref))^(c_pd / R_d)

    # calculate rho and T with p and theta (now perturbed) rho = p / R_d T, T = θ / π
    rho = p / ((p / p_0)^kappa * R_d * θ)
    T = p / (R_d * rho)

    v1 = 20
    v2 = 0
    rho_v1 = rho * v1
    rho_v2 = rho * v2
    rho_E = rho * c_vd * T + 1 / 2 * rho * (v1^2 + v2^2)
    return SVector(rho, rho_v1, rho_v2, rho_E, zero(eltype(g)), zero(eltype(g)))
end

initial_condition = initial_condition_warm_bubble

boundary_conditions = Dict(:y_neg => boundary_condition_slip_wall,
                           :y_pos => boundary_condition_slip_wall)

# Gravity source since Q_ph=0
source_term = source_terms_geopotential

###############################################################################
# Get the DG approximation space
polydeg = 4

surface_flux = FluxLMARS(360.0)
volume_flux = flux_chandrashekar

volume_integral = VolumeIntegralFluxDifferencing(volume_flux)

# Create DG solver with polynomial degree = 4 and LMARS flux as surface flux
# and the EC flux (chandrashekar) as volume flux
solver = DGSEM(polydeg, surface_flux, volume_integral)

coordinates_min = (0.0, -5000.0)
coordinates_max = (20000.0, 15000.0)

trees_per_dimension = (2, 2)
mesh = P4estMesh(trees_per_dimension, polydeg = 1,
                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
                 initial_refinement_level = 5,
                 periodicity = (true, false))

# A semidiscretization collects data structures and functions for the spatial discretization
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    boundary_conditions = boundary_conditions,
                                    source_terms = source_term)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 1000.0)

# Create ODE problem with time span from 0.0 to 1000.0
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
solution_variables = cons2drypot

analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
                                     extra_analysis_errors = (:entropy_conservation_error,))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 1000,
                                     save_initial_solution = true,
                                     save_final_solution = true,
                                     solution_variables = solution_variables)

stepsize_callback = StepsizeCallback(cfl = 0.2)

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)