-
Notifications
You must be signed in to change notification settings - Fork 98
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add subcell limiting support for P4estMesh #1954
base: main
Are you sure you want to change the base?
Changes from all commits
71fceb8
ea4c3e0
4b9d8b1
a0206b8
99e50e6
1dd6122
2641da3
f0d92e1
1250a94
e81d950
8a14a3d
29aa331
446c21c
a52c174
eadcf6f
ccbb71f
2bd75d3
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,101 @@ | ||
using OrdinaryDiffEq | ||
using Trixi | ||
|
||
############################################################################### | ||
# semidiscretization of the compressible Euler equations | ||
|
||
equations = CompressibleEulerEquations2D(1.4) | ||
|
||
""" | ||
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) | ||
E = 1.0 | ||
p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2) | ||
p0_outer = 1.0e-5 # = true Sedov 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 | ||
|
||
# Get the DG approximation space | ||
surface_flux = flux_lax_friedrichs | ||
volume_flux = flux_ranocha | ||
polydeg = 3 | ||
basis = LobattoLegendreBasis(polydeg) | ||
limiter_idp = SubcellLimiterIDP(equations, basis; | ||
local_twosided_variables_cons = ["rho"], | ||
local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, | ||
min)], | ||
max_iterations_newton = 40) # Default value of 10 iterations is too low to fulfill bounds. | ||
|
||
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; | ||
volume_flux_dg = volume_flux, | ||
volume_flux_fv = surface_flux) | ||
solver = DGSEM(basis, surface_flux, volume_integral) | ||
|
||
############################################################################### | ||
|
||
coordinates_min = (-1.0, -1.0) | ||
coordinates_max = (1.0, 1.0) | ||
|
||
trees_per_dimension = (4, 4) | ||
mesh = P4estMesh(trees_per_dimension, | ||
polydeg = polydeg, initial_refinement_level = 2, | ||
coordinates_min = coordinates_min, coordinates_max = coordinates_max, | ||
periodicity = true) | ||
|
||
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 = 300 | ||
analysis_callback = AnalysisCallback(semi, interval = analysis_interval) | ||
|
||
alive_callback = AliveCallback(analysis_interval = analysis_interval) | ||
|
||
save_solution = SaveSolutionCallback(interval = 300, | ||
save_initial_solution = true, | ||
save_final_solution = true) | ||
|
||
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()) | ||
|
||
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 | ||
callback = callbacks); | ||
summary_callback() # print the timer summary |
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,142 @@ | ||
# Channel flow around a cylinder at Mach 3 | ||
# | ||
# Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain | ||
# and supersonic outflow at the right portion of the domain. The top and bottom of the | ||
# channel as well as the cylinder are treated as Euler slip wall boundaries. | ||
# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz | ||
# instabilities at later times as two Mach stems form above and below the cylinder. | ||
# | ||
# For complete details on the problem setup see Section 5.7 of the paper: | ||
# - Jean-Luc Guermond, Murtazo Nazarov, Bojan Popov, and Ignacio Tomas (2018) | ||
# Second-Order Invariant Domain Preserving Approximation of the Euler Equations using Convex Limiting. | ||
# [DOI: 10.1137/17M1149961](https://doi.org/10.1137/17M1149961) | ||
# | ||
# Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D | ||
|
||
using OrdinaryDiffEq | ||
using Trixi | ||
|
||
############################################################################### | ||
# semidiscretization of the compressible Euler equations | ||
|
||
equations = CompressibleEulerEquations2D(1.4) | ||
|
||
@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D) | ||
# set the freestream flow parameters | ||
rho_freestream = 1.4 | ||
v1 = 3.0 | ||
v2 = 0.0 | ||
p_freestream = 1.0 | ||
|
||
prim = SVector(rho_freestream, v1, v2, p_freestream) | ||
return prim2cons(prim, equations) | ||
end | ||
|
||
initial_condition = initial_condition_mach3_flow | ||
|
||
# Supersonic inflow boundary condition. | ||
# Calculate the boundary flux entirely from the external solution state, i.e., set | ||
# external solution state values for everything entering the domain. | ||
@inline function boundary_condition_supersonic_inflow(u_inner, | ||
normal_direction::AbstractVector, | ||
x, t, surface_flux_function, | ||
equations::CompressibleEulerEquations2D) | ||
u_boundary = initial_condition_mach3_flow(x, t, equations) | ||
flux = Trixi.flux(u_boundary, normal_direction, equations) | ||
|
||
return flux | ||
end | ||
|
||
@inline function Trixi.get_boundary_outer_state(u_inner, t, | ||
boundary_condition::typeof(boundary_condition_supersonic_inflow), | ||
normal_direction::AbstractVector, direction, | ||
mesh::P4estMesh{2}, equations, dg, cache, | ||
indices...) | ||
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...) | ||
|
||
return initial_condition_mach3_flow(x, t, equations) | ||
end | ||
|
||
# Supersonic outflow boundary condition. | ||
# Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow | ||
# except all the solution state values are set from the internal solution as everything leaves the domain | ||
@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t, | ||
surface_flux_function, | ||
equations::CompressibleEulerEquations2D) | ||
flux = Trixi.flux(u_inner, normal_direction, equations) | ||
|
||
return flux | ||
end | ||
|
||
@inline function Trixi.get_boundary_outer_state(u_inner, t, | ||
boundary_condition::typeof(boundary_condition_outflow), | ||
orientation_or_normal, direction, | ||
mesh::P4estMesh{2}, equations, dg, cache, | ||
indices...) | ||
return u_inner | ||
end | ||
|
||
boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall, | ||
:Circle => boundary_condition_slip_wall, | ||
:Top => boundary_condition_slip_wall, | ||
:Right => boundary_condition_outflow, | ||
:Left => boundary_condition_supersonic_inflow) | ||
|
||
volume_flux = flux_ranocha_turbo | ||
surface_flux = flux_lax_friedrichs | ||
polydeg = 3 | ||
basis = LobattoLegendreBasis(polydeg) | ||
limiter_idp = SubcellLimiterIDP(equations, basis; | ||
local_twosided_variables_cons = ["rho"], | ||
positivity_variables_nonlinear = [pressure], | ||
local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, | ||
min)], | ||
max_iterations_newton = 50) # Default value of 10 iterations is too low to fulfill bounds. | ||
|
||
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; | ||
volume_flux_dg = volume_flux, | ||
volume_flux_fv = surface_flux) | ||
solver = DGSEM(basis, surface_flux, volume_integral) | ||
|
||
# Get the unstructured quad mesh from a file (downloads the file if not available locally) | ||
mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp", | ||
joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")) | ||
|
||
mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 0) | ||
|
||
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, | ||
boundary_conditions = boundary_conditions) | ||
|
||
############################################################################### | ||
# ODE solvers | ||
|
||
tspan = (0.0, 2.0) | ||
ode = semidiscretize(semi, tspan) | ||
|
||
# Callbacks | ||
|
||
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.8) | ||
|
||
callbacks = CallbackSet(summary_callback, | ||
analysis_callback, alive_callback, | ||
save_solution, | ||
stepsize_callback) | ||
|
||
stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) | ||
|
||
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 | ||
callback = callbacks); | ||
summary_callback() # print the timer summary |
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -5,7 +5,8 @@ | |
@muladd begin | ||
#! format: noindent | ||
|
||
@inline function check_bounds(u, equations, solver, cache, | ||
@inline function check_bounds(u::AbstractArray{<:Any, 4}, | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Why is this dispatch needed if you didn't add a new There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. I just wanted to make sure that an error is thrown when calling this routine with another dimension than 2. In the past, I had |
||
equations, solver, cache, | ||
limiter::SubcellLimiterIDP) | ||
(; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter | ||
(; variable_bounds) = limiter.cache.subcell_limiter_coefficients | ||
|
@@ -102,47 +103,4 @@ | |
|
||
return nothing | ||
end | ||
|
||
@inline function save_bounds_check_errors(output_directory, time, iter, equations, | ||
limiter::SubcellLimiterIDP) | ||
Comment on lines
-106
to
-107
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Note: Relocated this routine to |
||
(; local_twosided, positivity, local_onesided) = limiter | ||
(; idp_bounds_delta_local) = limiter.cache | ||
|
||
# Print to output file | ||
open(joinpath(output_directory, "deviations.txt"), "a") do f | ||
print(f, iter, ", ", time) | ||
if local_twosided | ||
for v in limiter.local_twosided_variables_cons | ||
v_string = string(v) | ||
print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")], | ||
", ", idp_bounds_delta_local[Symbol(v_string, "_max")]) | ||
end | ||
end | ||
if local_onesided | ||
for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear | ||
key = Symbol(string(variable), "_", string(min_or_max)) | ||
print(f, ", ", idp_bounds_delta_local[key]) | ||
end | ||
end | ||
if positivity | ||
for v in limiter.positivity_variables_cons | ||
if v in limiter.local_twosided_variables_cons | ||
continue | ||
end | ||
print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")]) | ||
end | ||
for variable in limiter.positivity_variables_nonlinear | ||
print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")]) | ||
end | ||
end | ||
println(f) | ||
end | ||
|
||
# Reset local maximum deviations | ||
for (key, _) in idp_bounds_delta_local | ||
idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key])) | ||
end | ||
|
||
return nothing | ||
end | ||
end # @muladd |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Note: Relocated this routine from
subcell_bounds_check_2d.jl
since it is not dimension-dependent.