Skip to content

Commit f0c6807

Browse files
patrickersingandrewwinters5000
patrickersing
and
andrewwinters5000
authored
Add Shallow water moment equation models in 1D (#128)
* add swme * add tests * add well-balanced test, remove allocations from dissipation term * add unit tests * remove unused import of diagm * fix unit test * add news entry * add Symbolics dependency for testing * fix BCs * fix unit_test * fix tests * update tests * update test reference * update wb test with wall BC * add unit test * add unit test * apply formatter * clean up comments * Apply suggestions from code review Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se> * re-introduce the union for wb test * switch to SArrays * update SArray for linearized swme --------- Co-authored-by: Andrew Winters <andrew.ross.winters@liu.se>
1 parent 729b1a8 commit f0c6807

14 files changed

+1995
-10
lines changed

‎NEWS.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,7 @@ for human readability.
99

1010
#### Added
1111
- Introduce node-wise limiting functionality for the `ShallowWaterMultiLayerEquations2D`. ([#111])
12+
- New equation types `ShallowWaterMomentEquations1D` and `ShallowWaterLinearizedMomentEquations1D` have been added. ([#128])
1213

1314
#### Changed
1415
- Velocity desingularization procedure has been moved into a distinct `VelocityDesingularization`

‎examples/tree_1d_dgsem/elixir_shallowwater_linearized_moments_convergence.jl‎

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
using TrixiShallowWater
4+
using Symbolics
5+
6+
###############################################################################
7+
# Semidiscretization of the shallow water linearized moment equations to test convergence against a
8+
# manufactured solution.
9+
10+
equations = ShallowWaterLinearizedMomentEquations1D(gravity = 9.81, n_moments = 2)
11+
12+
### Create manufactured solution for method of manufactured solutions (MMS)
13+
14+
# Symbolic Variables
15+
@variables x_sym, t_sym, g
16+
17+
# Define Differentials
18+
Dt, Dx = Differential(t_sym), Differential(x_sym[1])
19+
20+
## Initial condition
21+
###############################################################################
22+
# Primitive Variables
23+
H = 7 + cos(sqrt(2) * 2 * pi * x_sym) * cos(2 * pi * t_sym)
24+
v = 0.5
25+
b = 2 + 0.5 * sinpi(sqrt(2) * x_sym)
26+
h = H - b
27+
a = [1 for i in 1:(equations.n_moments)]
28+
29+
init = [H, v, a..., b]
30+
31+
## PDE
32+
###############################################################################
33+
# precompute the sum term
34+
sum_moments = sum(h * a[j]^2 / (2j + 1) for j in 1:(equations.n_moments))
35+
36+
# additional moment equations
37+
mom_eqs = [Dt(h * a[i]) + Dx(2 * h * v * a[i]) - v * Dx(h * a[i])
38+
for i in 1:(equations.n_moments)]
39+
40+
# PDE Source Terms
41+
eqs = [
42+
Dt(h) + Dx(h * v),
43+
Dt(h * v) + Dx(h * v^2 + sum_moments) + g * h * Dx(h + b),
44+
mom_eqs...,
45+
0
46+
]
47+
48+
## Create the functions for the manufactured solution
49+
########################################################################################
50+
# Expand derivatives
51+
du_exprs = expand_derivatives.(eqs)
52+
53+
# Build functions
54+
du_funcs = build_function.(du_exprs, Ref(x_sym), t_sym, g, expression = Val(false))
55+
56+
init_funcs = build_function.(init, Ref(x_sym), t_sym, expression = Val(false))
57+
58+
# Trixi functions
59+
function initial_condition_convergence(x,
60+
t,
61+
equations::ShallowWaterLinearizedMomentEquations1D)
62+
prim = SVector{3 + equations.n_moments, Float64}([f(x[1], t) for f in init_funcs]...)
63+
return prim2cons(prim, equations)
64+
end
65+
66+
function source_terms_convergence(u,
67+
x,
68+
t,
69+
equations::ShallowWaterLinearizedMomentEquations1D)
70+
g = equations.gravity
71+
return SVector{3 + equations.n_moments, Float64}([f(x[1], t, g) for f in du_funcs]...)
72+
end
73+
74+
initial_condition = initial_condition_convergence
75+
76+
###############################################################################
77+
# Get the DG approximation space
78+
volume_flux = (flux_careaga_etal, flux_nonconservative_careaga_etal)
79+
surface_flux = (FluxPlusDissipation(flux_careaga_etal,
80+
DissipationLaxFriedrichsEntropyVariables(max_abs_speed)),
81+
flux_nonconservative_careaga_etal)
82+
83+
solver = DGSEM(polydeg = 3,
84+
surface_flux = surface_flux,
85+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
86+
87+
###############################################################################
88+
# Get the TreeMesh and setup a periodic mesh
89+
90+
coordinates_min = 0.0
91+
coordinates_max = sqrt(2.0)
92+
mesh = TreeMesh(coordinates_min,
93+
coordinates_max,
94+
initial_refinement_level = 4,
95+
n_cells_max = 10_000,
96+
periodicity = true)
97+
98+
# create the semi discretization object
99+
semi = SemidiscretizationHyperbolic(mesh,
100+
equations,
101+
initial_condition,
102+
solver,
103+
source_terms = source_terms_convergence;
104+
boundary_conditions = boundary_condition_periodic)
105+
106+
###############################################################################
107+
# ODE solvers, callbacks etc.
108+
109+
tspan = (0.0, 0.05)
110+
ode = semidiscretize(semi, tspan)
111+
112+
summary_callback = SummaryCallback()
113+
114+
analysis_interval = 500
115+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
116+
117+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
118+
119+
save_solution = SaveSolutionCallback(interval = 200,
120+
save_initial_solution = true,
121+
save_final_solution = true)
122+
123+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
124+
125+
###############################################################################
126+
# run the simulation
127+
128+
# use a Runge-Kutta method with fixed step size
129+
sol = solve(ode,
130+
CarpenterKennedy2N54(williamson_condition = false);
131+
dt = 1e-4,
132+
ode_default_options()...,
133+
callback = callbacks,);

‎examples/tree_1d_dgsem/elixir_shallowwater_moments_convergence.jl‎

Lines changed: 131 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,131 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
using TrixiShallowWater
4+
using Symbolics
5+
6+
###############################################################################
7+
# Semidiscretization of the shallow water moment equations to test convergence against a
8+
# manufactured solution.
9+
10+
equations = ShallowWaterMomentEquations1D(gravity = 9.81, n_moments = 2)
11+
12+
### Create manufactured solution for method of manufactured solutions (MMS)
13+
n = equations.n_moments
14+
15+
# Symbolic Variables
16+
@variables x_sym, t_sym, g
17+
18+
# Define Differentials
19+
Dt, Dx = Differential(t_sym), Differential(x_sym)
20+
21+
## Initial condition
22+
##################################################################################################
23+
H = 7 + cos(sqrt(2) * 2 * pi * x_sym) * cos(2 * pi * t_sym)
24+
v = 0.5
25+
b = 2 + 0.5 * sinpi(sqrt(2) * x_sym)
26+
h = H - b
27+
a = [0.5 for i in 1:n]
28+
29+
init = [H, v, a..., b]
30+
31+
## PDE
32+
###################################################################################################
33+
# precompute the sum term
34+
sum_moments = sum(h * a[j]^2 / (2j + 1) for j in 1:n)
35+
sum_A = [sum(h * equations.A[i, j, k] * a[j] * a[k] for j in 1:n, k in 1:n) for i in 1:n]
36+
sum_B = [sum(equations.B[i, j, k] * a[k] * Dx(h * a[j]) for j in 1:n, k in 1:n)
37+
for i in 1:n]
38+
39+
# additional moment equations
40+
mom_eqs = [Dt(h * a[i]) + Dx(2 * h * v * a[i] + sum_A[i]) - v * Dx(h * a[i]) + sum_B[i]
41+
for i in 1:n]
42+
43+
# PDE Source Terms
44+
eqs = [
45+
Dt(h) + Dx(h * v),
46+
Dt(h * v) + Dx(h * v^2 + sum_moments) + g * h * Dx(h + b),
47+
mom_eqs...,
48+
0
49+
]
50+
51+
## Create the functions for the manufactured solution
52+
###################################################################################################
53+
# Expand derivatives
54+
du_exprs = expand_derivatives.(eqs)
55+
56+
# Build functions
57+
du_funcs = build_function.(du_exprs, Ref(x_sym), Ref(t_sym), g, expression = Val(false))
58+
59+
init_funcs = build_function.(init, Ref(x_sym), t_sym, expression = Val(false))
60+
61+
# Trixi functions
62+
function initial_condition_convergence(x, t, equations::ShallowWaterMomentEquations1D)
63+
prim = SVector{3 + n, Float64}([f(x[1], t) for f in init_funcs]...)
64+
return prim2cons(prim, equations)
65+
end
66+
67+
function source_terms_convergence(u, x, t, equations::ShallowWaterMomentEquations1D)
68+
g = equations.gravity
69+
return SVector{3 + n, Float64}([f(x[1], t, g) for f in du_funcs]...)
70+
end
71+
72+
initial_condition = initial_condition_convergence
73+
74+
###############################################################################
75+
# Get the DG approximation space
76+
volume_flux = (flux_careaga_etal, flux_nonconservative_careaga_etal)
77+
surface_flux = (FluxPlusDissipation(flux_careaga_etal,
78+
DissipationLaxFriedrichsEntropyVariables(max_abs_speed)),
79+
flux_nonconservative_careaga_etal)
80+
81+
solver = DGSEM(polydeg = 3,
82+
surface_flux = surface_flux,
83+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
84+
85+
###############################################################################
86+
# Get the TreeMesh and setup a periodic mesh
87+
88+
coordinates_min = 0.0
89+
coordinates_max = sqrt(2.0)
90+
mesh = TreeMesh(coordinates_min,
91+
coordinates_max,
92+
initial_refinement_level = 4,
93+
n_cells_max = 10_000,
94+
periodicity = true)
95+
96+
# create the semi discretization object
97+
semi = SemidiscretizationHyperbolic(mesh,
98+
equations,
99+
initial_condition,
100+
solver,
101+
source_terms = source_terms_convergence;
102+
boundary_conditions = boundary_condition_periodic)
103+
104+
###############################################################################
105+
# ODE solvers, callbacks etc.
106+
107+
tspan = (0.0, 0.05)
108+
ode = semidiscretize(semi, tspan)
109+
110+
summary_callback = SummaryCallback()
111+
112+
analysis_interval = 500
113+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
114+
115+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
116+
117+
save_solution = SaveSolutionCallback(interval = 200,
118+
save_initial_solution = true,
119+
save_final_solution = true)
120+
121+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution)
122+
123+
###############################################################################
124+
# run the simulation
125+
126+
# use a Runge-Kutta method fixed step size
127+
sol = solve(ode,
128+
CarpenterKennedy2N54(williamson_condition = false);
129+
dt = 1e-4,
130+
ode_default_options()...,
131+
callback = callbacks,);

‎examples/tree_1d_dgsem/elixir_shallowwater_moments_smooth_wave.jl‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
2+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
3+
using Trixi
4+
using TrixiShallowWater
5+
6+
# Semidiscretization of the shallow water moment equations with two moments
7+
equations = ShallowWaterMomentEquations1D(gravity = 1.0, n_moments = 2)
8+
9+
# Initial condition with a smooth wave wave in a periodic domain.
10+
# See section 4.2 in the paper:
11+
# Julian Koellermeier and Marvin Rominger (2020)
12+
# "Analysis and Numerical Simulation of Hyperbolic Shallow Water Moment Equations"
13+
# [DOI: 10.4208/cicp.OA-2019-0065](https://doi.org/10.4208/cicp.OA-2019-0065)
14+
function initial_condition_smooth_periodic_wave(x, t,
15+
equations::Union{ShallowWaterMomentEquations1D,
16+
ShallowWaterLinearizedMomentEquations1D})
17+
H = 1.0 + exp(3.0 * cos(π * (x[1] + 0.5))) / exp(4.0)
18+
v = 0.25
19+
a = zeros(equations.n_moments)
20+
a[2] = -0.25
21+
22+
b = 0.0
23+
24+
return prim2cons(SVector(H, v, a..., b), equations)
25+
end
26+
27+
initial_condition = initial_condition_smooth_periodic_wave
28+
29+
###############################################################################
30+
# Get the DG approximation space
31+
32+
volume_flux = (flux_careaga_etal, flux_nonconservative_careaga_etal)
33+
surface_flux = (FluxPlusDissipation(flux_careaga_etal,
34+
DissipationLaxFriedrichsEntropyVariables(Trixi.max_abs_speed)),
35+
flux_nonconservative_careaga_etal)
36+
37+
indicator_var(u, equations) = u[2]^3
38+
39+
basis = LobattoLegendreBasis(3)
40+
indicator_sc = IndicatorHennemannGassner(equations, basis,
41+
alpha_max = 0.5,
42+
alpha_min = 0.001,
43+
alpha_smooth = true,
44+
variable = indicator_var)
45+
46+
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
47+
volume_flux_dg = volume_flux,
48+
volume_flux_fv = surface_flux,)
49+
50+
solver = DGSEM(basis, surface_flux, volume_integral)
51+
52+
###############################################################################
53+
# Create the TreeMesh for the domain [-1, 1]
54+
55+
coordinates_min = -1.0
56+
coordinates_max = 1.0
57+
58+
mesh = TreeMesh(coordinates_min,
59+
coordinates_max,
60+
initial_refinement_level = 6, # 2^refinement_level
61+
n_cells_max = 10_000,
62+
periodicity = true)
63+
64+
# create the semi discretization object
65+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
66+
source_terms = source_term_manning_friction;
67+
boundary_conditions = boundary_condition_periodic)
68+
69+
###############################################################################
70+
# ODE solver
71+
tspan = (0.0, 1.0)
72+
ode = semidiscretize(semi, tspan)
73+
74+
###############################################################################
75+
# Callbacks
76+
summary_callback = SummaryCallback()
77+
78+
analysis_interval = 200
79+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
80+
save_analysis = false)
81+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
82+
stepsize_callback = StepsizeCallback(cfl = 0.9)
83+
84+
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback,
85+
stepsize_callback)
86+
###############################################################################
87+
# run the simulation
88+
89+
sol = solve(ode,
90+
CarpenterKennedy2N54(williamson_condition = false);
91+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
92+
ode_default_options()...,
93+
callback = callbacks,);

0 commit comments

Comments
 (0)