-
Notifications
You must be signed in to change notification settings - Fork 99
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
31 changed files
with
2,612 additions
and
73 deletions.
There are no files selected for viewing
114 changes: 114 additions & 0 deletions
114
examples/structured_2d_dgsem/elixir_shallowwater_conical_island.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,114 @@ | ||
|
||
using OrdinaryDiffEq | ||
using Trixi | ||
|
||
############################################################################### | ||
# Semidiscretization of the shallow water equations | ||
# | ||
# TODO: TrixiShallowWater: wet/dry example elixir | ||
|
||
equations = ShallowWaterEquations2D(gravity_constant=9.81, H0=1.4) | ||
|
||
""" | ||
initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D) | ||
Initial condition for the [`ShallowWaterEquations2D`](@ref) to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) | ||
and its handling of discontinuous water heights at the start in combination with wetting and | ||
drying. The bottom topography is given by a conical island in the middle of the domain. Around that | ||
island, there is a cylindrical water column at t=0 and the rest of the domain is dry. This | ||
discontinuous water height is smoothed by a logistic function. This simulation uses periodic | ||
boundary conditions. | ||
""" | ||
function initial_condition_conical_island(x, t, equations::ShallowWaterEquations2D) | ||
# Set the background values | ||
|
||
v1 = 0.0 | ||
v2 = 0.0 | ||
|
||
x1, x2 = x | ||
b = max(0.1, 1.0 - 4.0 * sqrt(x1^2 + x2^2)) | ||
|
||
# use a logistic function to transfer water height value smoothly | ||
L = equations.H0 # maximum of function | ||
x0 = 0.3 # center point of function | ||
k = -25.0 # sharpness of transfer | ||
|
||
H = max(b, L/(1.0 + exp(-k*(sqrt(x1^2+x2^2) - x0)))) | ||
|
||
# It is mandatory to shift the water level at dry areas to make sure the water height h | ||
# stays positive. The system would not be stable for h set to a hard 0 due to division by h in | ||
# the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold | ||
# with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above | ||
# for the ShallowWaterEquations and added to the initial condition if h = 0. | ||
# This default value can be changed within the constructor call depending on the simulation setup. | ||
H = max(H, b + equations.threshold_limiter) | ||
return prim2cons(SVector(H, v1, v2, b), equations) | ||
end | ||
|
||
initial_condition = initial_condition_conical_island | ||
|
||
############################################################################### | ||
# Get the DG approximation space | ||
|
||
volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) | ||
surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), | ||
flux_nonconservative_chen_noelle) | ||
|
||
basis = LobattoLegendreBasis(4) | ||
|
||
indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, | ||
alpha_max=0.5, | ||
alpha_min=0.001, | ||
alpha_smooth=true, | ||
variable=waterheight_pressure) | ||
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; | ||
volume_flux_dg=volume_flux, | ||
volume_flux_fv=surface_flux) | ||
|
||
solver = DGSEM(basis, surface_flux, volume_integral) | ||
|
||
############################################################################### | ||
# Get the StructuredMesh and setup a periodic mesh | ||
|
||
coordinates_min = (-1.0, -1.0) | ||
coordinates_max = (1.0, 1.0) | ||
|
||
cells_per_dimension = (16, 16) | ||
|
||
mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) | ||
|
||
# Create the semi discretization object | ||
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) | ||
|
||
############################################################################### | ||
# ODE solver | ||
|
||
tspan = (0.0, 10.0) | ||
ode = semidiscretize(semi, tspan) | ||
|
||
############################################################################### | ||
# Callbacks | ||
|
||
summary_callback = SummaryCallback() | ||
|
||
analysis_interval = 1000 | ||
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, | ||
extra_analysis_integrals=(lake_at_rest_error,)) | ||
|
||
alive_callback = AliveCallback(analysis_interval=analysis_interval) | ||
|
||
save_solution = SaveSolutionCallback(interval=100, | ||
save_initial_solution=true, | ||
save_final_solution=true) | ||
|
||
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) | ||
|
||
############################################################################### | ||
# run the simulation | ||
|
||
stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,)) | ||
|
||
sol = solve(ode, SSPRK43(stage_limiter!); | ||
ode_default_options()..., callback=callbacks); | ||
|
||
summary_callback() # print the timer summary |
120 changes: 120 additions & 0 deletions
120
examples/structured_2d_dgsem/elixir_shallowwater_parabolic_bowl.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,120 @@ | ||
|
||
using OrdinaryDiffEq | ||
using Trixi | ||
|
||
############################################################################### | ||
# Semidiscretization of the shallow water equations | ||
# | ||
# TODO: TrixiShallowWater: wet/dry example elixir | ||
|
||
equations = ShallowWaterEquations2D(gravity_constant=9.81) | ||
|
||
""" | ||
initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D) | ||
Well-known initial condition to test the [`hydrostatic_reconstruction_chen_noelle`](@ref) and its | ||
wet-dry mechanics. This test has an analytical solution. The initial condition is defined by the | ||
analytical solution at time t=0. The bottom topography defines a bowl and the water level is given | ||
by an oscillating lake. | ||
The original test and its analytical solution were first presented in | ||
- William C. Thacker (1981) | ||
Some exact solutions to the nonlinear shallow-water wave equations | ||
[DOI: 10.1017/S0022112081001882](https://doi.org/10.1017/S0022112081001882). | ||
The particular setup below is taken from Section 6.2 of | ||
- Niklas Wintermeyer, Andrew R. Winters, Gregor J. Gassner and Timothy Warburton (2018) | ||
An entropy stable discontinuous Galerkin method for the shallow water equations on | ||
curvilinear meshes with wet/dry fronts accelerated by GPUs | ||
[DOI: 10.1016/j.jcp.2018.08.038](https://doi.org/10.1016/j.jcp.2018.08.038). | ||
""" | ||
function initial_condition_parabolic_bowl(x, t, equations:: ShallowWaterEquations2D) | ||
a = 1.0 | ||
h_0 = 0.1 | ||
sigma = 0.5 | ||
ω = sqrt(2 * equations.gravity * h_0) / a | ||
|
||
v1 = -sigma * ω * sin(ω * t) | ||
v2 = sigma * ω * cos(ω * t) | ||
|
||
b = h_0 * ((x[1])^2 + (x[2])^2) / a^2 | ||
|
||
H = sigma * h_0 / a^2 * (2 * x[1] * cos(ω * t) + 2 * x[2] * sin(ω * t) - sigma) + h_0 | ||
|
||
# It is mandatory to shift the water level at dry areas to make sure the water height h | ||
# stays positive. The system would not be stable for h set to a hard 0 due to division by h in | ||
# the computation of velocity, e.g., (h v1) / h. Therefore, a small dry state threshold | ||
# with a default value of 500*eps() ≈ 1e-13 in double precision, is set in the constructor above | ||
# for the ShallowWaterEquations and added to the initial condition if h = 0. | ||
# This default value can be changed within the constructor call depending on the simulation setup. | ||
H = max(H, b + equations.threshold_limiter) | ||
return prim2cons(SVector(H, v1, v2, b), equations) | ||
end | ||
|
||
initial_condition = initial_condition_parabolic_bowl | ||
|
||
|
||
############################################################################### | ||
# Get the DG approximation space | ||
|
||
volume_flux = (flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal) | ||
surface_flux = (FluxHydrostaticReconstruction(flux_hll_chen_noelle, hydrostatic_reconstruction_chen_noelle), | ||
flux_nonconservative_chen_noelle) | ||
|
||
basis = LobattoLegendreBasis(4) | ||
|
||
indicator_sc = IndicatorHennemannGassnerShallowWater(equations, basis, | ||
alpha_max=0.6, | ||
alpha_min=0.001, | ||
alpha_smooth=true, | ||
variable=waterheight_pressure) | ||
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc; | ||
volume_flux_dg=volume_flux, | ||
volume_flux_fv=surface_flux) | ||
|
||
solver = DGSEM(basis, surface_flux, volume_integral) | ||
|
||
|
||
############################################################################### | ||
|
||
coordinates_min = (-2.0, -2.0) | ||
coordinates_max = (2.0, 2.0) | ||
|
||
cells_per_dimension = (150, 150) | ||
|
||
mesh = StructuredMesh(cells_per_dimension, coordinates_min, coordinates_max) | ||
|
||
# create the semi discretization object | ||
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) | ||
|
||
############################################################################### | ||
# ODE solvers, callbacks etc. | ||
|
||
tspan = (0.0, 1.0) | ||
ode = semidiscretize(semi, tspan) | ||
|
||
summary_callback = SummaryCallback() | ||
|
||
analysis_interval = 1000 | ||
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=false, | ||
extra_analysis_integrals=(energy_kinetic, | ||
energy_internal, | ||
lake_at_rest_error)) | ||
|
||
alive_callback = AliveCallback(analysis_interval=analysis_interval) | ||
|
||
save_solution = SaveSolutionCallback(interval=100, | ||
save_initial_solution=true, | ||
save_final_solution=true) | ||
|
||
callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution) | ||
|
||
stage_limiter! = PositivityPreservingLimiterShallowWater(variables=(Trixi.waterheight,)) | ||
|
||
############################################################################### | ||
# run the simulation | ||
|
||
sol = solve(ode, SSPRK43(stage_limiter!); | ||
ode_default_options()..., callback=callbacks); | ||
|
||
summary_callback() # print the timer summary |
Oops, something went wrong.