Skip to content

Commit

Permalink
Merge branch 'main' into HLL_Davis_Polytropic_Euler
Browse files Browse the repository at this point in the history
  • Loading branch information
DanielDoehring authored Nov 13, 2023
2 parents 1bd1b56 + 7bb3b46 commit 945c336
Show file tree
Hide file tree
Showing 17 changed files with 744 additions and 30 deletions.
3 changes: 2 additions & 1 deletion NEWS.md
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,8 @@ for human readability.
- Wetting and drying feature and examples for 1D and 2D shallow water equations
- Implementation of the polytropic Euler equations in 2D
- Implementation of the quasi-1D shallow water equations
- Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
- Subcell (positivity and local min/max) limiting support for conservative variables
in 2D for `TreeMesh`
- AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`
- Added `GradientVariables` type parameter to `AbstractEquationsParabolic`

Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

"""
initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
A medium blast wave taken from
- Sebastian Hennemann, Gregor J. Gassner (2020)
A provably entropy stable subcell shock capturing approach for high order split form DG
[arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)
# Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave"
# Set up polar coordinates
inicenter = SVector(0.0, 0.0)
x_norm = x[1] - inicenter[1]
y_norm = x[2] - inicenter[2]
r = sqrt(x_norm^2 + y_norm^2)
phi = atan(y_norm, x_norm)
sin_phi, cos_phi = sincos(phi)

# Calculate primitive variables
rho = r > 0.5 ? 1.0 : 1.1691
v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi
v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi
p = r > 0.5 ? 1.0E-3 : 1.245

return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_blast_wave

boundary_condition = BoundaryConditionDirichlet(initial_condition)

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
basis = LobattoLegendreBasis(3)
limiter_idp = SubcellLimiterIDP(equations, basis;
local_minmax_variables_cons = ["rho"])
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
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)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 6,
n_cells_max = 10_000,
periodicity = false)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_condition)

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

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 100,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 0.3)

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

###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))

sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
91 changes: 91 additions & 0 deletions examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,91 @@

using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler equations
gamma = 1.4
equations = CompressibleEulerEquations2D(gamma)

"""
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
The Sedov blast wave setup based on Flash
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
"""
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
# Set up polar coordinates
inicenter = SVector(0.0, 0.0)
x_norm = x[1] - inicenter[1]
y_norm = x[2] - inicenter[2]
r = sqrt(x_norm^2 + y_norm^2)

# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
# r0 = 0.5 # = more reasonable setup
E = 1.0
p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
p0_outer = 1.0e-5 # = true Sedov setup
# p0_outer = 1.0e-3 # = more reasonable setup

# Calculate primitive variables
rho = 1.0
v1 = 0.0
v2 = 0.0
p = r > r0 ? p0_outer : p0_inner

return prim2cons(SVector(rho, v1, v2, p), equations)
end
initial_condition = initial_condition_sedov_blast_wave

surface_flux = flux_lax_friedrichs
volume_flux = flux_chandrashekar
basis = LobattoLegendreBasis(3)
limiter_idp = SubcellLimiterIDP(equations, basis;
local_minmax_variables_cons = ["rho"])
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
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)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 3,
n_cells_max = 100_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

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

tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 1000
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

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

stepsize_callback = StepsizeCallback(cfl = 0.6)

callbacks = CallbackSet(summary_callback,
analysis_callback, alive_callback,
stepsize_callback,
save_solution)
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))

sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
Original file line number Diff line number Diff line change
Expand Up @@ -39,7 +39,7 @@ surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
basis = LobattoLegendreBasis(3)
limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = [1],
positivity_variables_cons = ["rho"],
positivity_correction_factor = 0.5)
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible Euler multicomponent equations

# 1) Dry Air 2) Helium + 28% Air
equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
gas_constants = (0.287, 1.578))

"""
initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
A shock-bubble testcase for multicomponent Euler equations
- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
[arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
"""
function initial_condition_shock_bubble(x, t,
equations::CompressibleEulerMulticomponentEquations2D{
5,
2
})
# bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
# other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
# typical domain is rectangular, we change it to a square, as Trixi can only do squares
@unpack gas_constants = equations

# Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
delta = 0.03

# Region I
rho1_1 = delta
rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta
v1_1 = zero(delta)
v2_1 = zero(delta)
p_1 = 101325

# Region II
rho1_2 = 1.225 - delta
rho2_2 = delta
v1_2 = zero(delta)
v2_2 = zero(delta)
p_2 = 101325

# Region III
rho1_3 = 1.6861 - delta
rho2_3 = delta
v1_3 = -113.5243
v2_3 = zero(delta)
p_3 = 159060

# Set up Region I & II:
inicenter = SVector(zero(delta), zero(delta))
x_norm = x[1] - inicenter[1]
y_norm = x[2] - inicenter[2]
r = sqrt(x_norm^2 + y_norm^2)

if (x[1] > 0.50)
# Set up Region III
rho1 = rho1_3
rho2 = rho2_3
v1 = v1_3
v2 = v2_3
p = p_3
elseif (r < 0.25)
# Set up Region I
rho1 = rho1_1
rho2 = rho2_1
v1 = v1_1
v2 = v2_1
p = p_1
else
# Set up Region II
rho1 = rho1_2
rho2 = rho2_2
v1 = v1_2
v2 = v2_2
p = p_2
end

return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
end
initial_condition = initial_condition_shock_bubble

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
basis = LobattoLegendreBasis(3)

limiter_idp = SubcellLimiterIDP(equations, basis;
local_minmax_variables_cons = ["rho" * string(i)
for i in eachcomponent(equations)])
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)

solver = DGSEM(basis, surface_flux, volume_integral)

coordinates_min = (-2.25, -2.225)
coordinates_max = (2.20, 2.225)
mesh = TreeMesh(coordinates_min, coordinates_max,
initial_refinement_level = 3,
n_cells_max = 1_000_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

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

tspan = (0.0, 0.01)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 300
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
extra_analysis_integrals = (Trixi.density,))

alive_callback = AliveCallback(analysis_interval = analysis_interval)

save_solution = SaveSolutionCallback(interval = 600,
save_initial_solution = true,
save_final_solution = true,
solution_variables = cons2prim)

stepsize_callback = StepsizeCallback(cfl = 0.5)

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

###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))

sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
save_everystep = false, callback = callbacks);
summary_callback() # print the timer summary
Original file line number Diff line number Diff line change
Expand Up @@ -88,9 +88,8 @@ volume_flux = flux_ranocha
basis = LobattoLegendreBasis(3)

limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = [
(i + 3 for i in eachcomponent(equations))...,
])
positivity_variables_cons = ["rho" * string(i)
for i in eachcomponent(equations)])

volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@ volume_flux = (flux_derigs_etal, flux_nonconservative_powell_local_symmetric)
basis = LobattoLegendreBasis(3)

limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = [1],
positivity_variables_cons = ["rho"],
positivity_correction_factor = 0.5)
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
volume_flux_dg = volume_flux,
Expand Down
Loading

0 comments on commit 945c336

Please sign in to comment.