Skip to content

Commit

Permalink
Merge branch from subcell limiting p4est PR (#129)
Browse files Browse the repository at this point in the history 
Make sure I adapted everything correctly.
  • Loading branch information
@bennibolm
bennibolm authored Jul 29, 2024
1 parent 91b5db9 commit 0ac57af
Show file tree
Hide file tree
Showing 16 changed files with 362 additions and 170 deletions.
9 changes: 4 additions & 5 deletions examples/p4est_2d_dgsem/elixir_euler_double_mach.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond
return flux
end

# Note: Only for P4estMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@inline function Trixi.get_boundary_outer_state(u_inner, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
normal_direction::AbstractVector,
mesh::P4estMesh{2},
equations::CompressibleEulerEquations2D,
dg, indices...)
dg, cache, indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
if x[1] < 1 / 6 # BoundaryConditionCharacteristic
u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
Expand Down Expand Up @@ -160,7 +159,7 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
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
Expand Down
9 changes: 4 additions & 5 deletions examples/p4est_2d_dgsem/elixir_euler_double_mach_MCL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,13 +71,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond
return flux
end

# Note: Only for P4estMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@inline function Trixi.get_boundary_outer_state(u_inner, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
normal_direction::AbstractVector,
mesh::P4estMesh{2},
equations::CompressibleEulerEquations2D,
dg, indices...)
dg, cache, indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
if x[1] < 1 / 6 # BoundaryConditionCharacteristic
u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
Expand Down Expand Up @@ -162,7 +161,7 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (BoundsCheckCallback(save_errors = false),)
stage_callbacks = (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
Expand Down
102 changes: 102 additions & 0 deletions examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,102 @@
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.
bar_states = false)

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
12 changes: 6 additions & 6 deletions examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_MCL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,10 +48,10 @@ initial_condition = initial_condition_mach3_flow
return flux
end

@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@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,
normal_direction::AbstractVector,
mesh::P4estMesh{2}, equations, dg, cache,
indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)

Expand All @@ -69,10 +69,10 @@ end
return flux
end

@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@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,
orientation_or_normal,
mesh::P4estMesh{2}, equations, dg, cache,
indices...)
return u_inner
end
Expand Down
39 changes: 18 additions & 21 deletions examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@
#
# Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D

using Downloads: download
using OrdinaryDiffEq
using Trixi

Expand Down Expand Up @@ -48,10 +47,10 @@ initial_condition = initial_condition_mach3_flow
return flux
end

@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@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,
normal_direction::AbstractVector,
mesh::P4estMesh{2}, equations, dg, cache,
indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)

Expand All @@ -69,10 +68,10 @@ end
return flux
end

@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@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,
normal_direction::AbstractVector,
mesh::P4estMesh{2}, equations, dg, cache,
indices...)
return u_inner
end
Expand All @@ -91,28 +90,26 @@ boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,
:Right => boundary_condition_outflow,
:Left => boundary_condition_supersonic_inflow)

volume_flux = flux_ranocha_turbo
surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
polydeg = 3
basis = LobattoLegendreBasis(polydeg)
limiter_idp = SubcellLimiterIDP(equations, basis;
positivity_variables_cons = ["rho"],
local_twosided_variables_cons = ["rho"],
positivity_variables_nonlinear = [pressure],
local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
min)],
max_iterations_newton = 100,
max_iterations_newton = 50, # Default value of 10 iterations is too low to fulfill bounds.
bar_states = false)

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)
default_mesh_file = joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")
isfile(default_mesh_file) ||
download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp",
default_mesh_file)
mesh_file = default_mesh_file
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)

Expand All @@ -134,21 +131,21 @@ analysis_callback = AnalysisCallback(semi, interval = analysis_interval)

alive_callback = AliveCallback(analysis_interval = analysis_interval)

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

stepsize_callback = StepsizeCallback(cfl = 0.4)
stepsize_callback = StepsizeCallback(cfl = 0.8)

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

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
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
save_everystep = false, callback = callbacks);
callback = callbacks);
summary_callback() # print the timer summary
7 changes: 3 additions & 4 deletions examples/structured_2d_dgsem/elixir_euler_double_mach.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond
return flux
end

# Note: Only for StructuredMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@inline function Trixi.get_boundary_outer_state(u_inner, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
mesh::StructuredMesh{2},
equations::CompressibleEulerEquations2D,
dg, indices...)
dg, cache, indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
if x[1] < 1 / 6 # BoundaryConditionCharacteristic
u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
Expand Down Expand Up @@ -162,7 +161,7 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
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
Expand Down
7 changes: 3 additions & 4 deletions examples/structured_2d_dgsem/elixir_euler_double_mach_MCL.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,12 @@ boundary_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_cond
return flux
end

# Note: Only for StructuredMesh
@inline function Trixi.get_boundary_outer_state(u_inner, cache, t,
@inline function Trixi.get_boundary_outer_state(u_inner, t,
boundary_condition::typeof(boundary_condition_mixed_characteristic_wall),
normal_direction::AbstractVector, direction,
mesh::StructuredMesh{2},
equations::CompressibleEulerEquations2D,
dg, indices...)
dg, cache, indices...)
x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
if x[1] < 1 / 6 # BoundaryConditionCharacteristic
u_outer = Trixi.characteristic_boundary_value_function(initial_condition_double_mach_reflection,
Expand Down Expand Up @@ -164,7 +163,7 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (BoundsCheckCallback(save_errors = false),)
stage_callbacks = (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
Expand Down
Loading

0 comments on commit 0ac57af

Please sign in to comment.