Skip to content
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.

Sign up for GitHub

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

Merge branch from subcell limiting p4est PR #129

Merged
merged 2 commits into from
Jul 29, 2024
+362 −170
Merged

Merge branch from subcell limiting p4est PR #129

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
47d6d8e
Merge pestmesh pr branch
bennibolm Jul 29, 2024
1986068
Add parameter
bennibolm Jul 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 4 additions & 5 deletions examples/p4est_2d_dgsem/elixir_euler_double_mach.jl
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Loading