Skip to content

Commit

Permalink
tests and fmt
Browse files Browse the repository at this point in the history
  • Loading branch information
@DanielDoehring
DanielDoehring committed Jan 5, 2024
1 parent c783247 commit 06fb60c
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 55 deletions.
54 changes: 19 additions & 35 deletions examples/p4est_2d_dgsem/elixir_euler_airfoil_mach2.jl
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@

using Trixi
using OrdinaryDiffEq, Plots
using OrdinaryDiffEq
using Downloads: download

###############################################################################
Expand All @@ -9,14 +9,14 @@ using Downloads: download
equations = CompressibleEulerEquations2D(1.4)

@inline function initial_condition_mach2_flow(x, t, equations::CompressibleEulerEquations2D)
# set the freestream flow parameters
rho_freestream = 1.4
v1 = 2.0
v2 = 0.0
p_freestream = 1.0

prim = SVector(rho_freestream, v1, v2, p_freestream)
return prim2cons(prim, equations)
# set the freestream flow parameters
rho_freestream = 1.4
v1 = 2.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_mach2_flow
Expand All @@ -28,10 +28,10 @@ initial_condition = initial_condition_mach2_flow
normal_direction::AbstractVector,
x, t, surface_flux_function,
equations::CompressibleEulerEquations2D)
u_boundary = initial_condition_mach2_flow(x, t, equations)
flux = Trixi.flux(u_boundary, normal_direction, equations)
u_boundary = initial_condition_mach2_flow(x, t, equations)
flux = Trixi.flux(u_boundary, normal_direction, equations)

return flux
return flux
end

# Supersonic outflow boundary condition.
Expand All @@ -40,15 +40,15 @@ end
@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)
flux = Trixi.flux(u_inner, normal_direction, equations)

return flux
return flux
end

d = 3

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha
volume_flux = flux_ranocha

basis = LobattoLegendreBasis(d)
shock_indicator = IndicatorHennemannGassner(equations, basis,
Expand All @@ -67,17 +67,16 @@ solver = DGSEM(polydeg = d, surface_flux = surface_flux,
mesh_file = joinpath(@__DIR__, "NACA6412.inp")
isfile(mesh_file) ||
download("https://gist.github.com/DanielDoehring/b434005800a1b0c0b4b50a8772283019/raw/1f6649b3370163d75d268176b51775f8685dd81d/NACA6412.inp",
mesh_file)
mesh_file)

boundary_symbols = [:PhysicalLine10, :PhysicalLine20, :PhysicalLine30, :PhysicalLine40]

mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols=boundary_symbols)
mesh = P4estMesh{2}(mesh_file, polydeg = d, boundary_symbols = boundary_symbols)

boundary_conditions = Dict(:PhysicalLine10 => boundary_condition_supersonic_inflow,
:PhysicalLine20 => boundary_condition_outflow, # Right boundary
:PhysicalLine30 => boundary_condition_slip_wall, # Airfoil
:PhysicalLine40 => boundary_condition_outflow # Top/Bottom
)
:PhysicalLine40 => boundary_condition_outflow)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
boundary_conditions = boundary_conditions)
Expand All @@ -96,25 +95,10 @@ callbacks = CallbackSet(summary_callback,
analysis_callback,
stepsize_callback)

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

# Positivity limiter might be necessary for examples with strong shocks. Very sensitive
# to the order of the limiter variables, pressure must come first.
stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-7, 1.0e-6),
variables = (pressure, Trixi.density))

# Run the simulation
###############################################################################
sol = solve(ode, #SSPRK104(; thread = OrdinaryDiffEq.True(), stage_limiter!);
SSPRK104(; thread = OrdinaryDiffEq.True());
sol = solve(ode, SSPRK104(; thread = OrdinaryDiffEq.True());
dt = 42.0, # overwritten
callback = callbacks);

summary_callback() # print the timer summary

plot(sol)

pd = PlotData2D(sol)
plot(pd["rho"])
plot!(getmesh(pd))
61 changes: 41 additions & 20 deletions src/meshes/p4est_mesh.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -482,7 +482,7 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
node_set_dict = parse_node_sets(meshfile, boundary_symbols)
# Read in all elements with associated nodes to specify the boundaries
element_node_matrix = parse_elements(meshfile, n_trees, n_dimensions)

# Initialize boundary information matrix with symbol for no boundary / internal connection
boundary_names = fill(Symbol("---"), 2 * n_dimensions, n_trees)

Expand All @@ -494,22 +494,26 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
# and search for "Node ordering and face numbering on elements"
for boundary in keys(node_set_dict)
# Check bottom edge
if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary]
if tree_nodes[1] in node_set_dict[boundary] &&
tree_nodes[2] in node_set_dict[boundary]
# Bottom boundary is position 3 in p4est indexing
boundary_names[3, tree] = boundary
end
# Check right edge
if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary]
if tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[3] in node_set_dict[boundary]
# Right boundary is position 2 in p4est indexing
boundary_names[2, tree] = boundary
end
# Check top edge
if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary]
if tree_nodes[3] in node_set_dict[boundary] &&
tree_nodes[4] in node_set_dict[boundary]
# Top boundary is position 4 in p4est indexing
boundary_names[4, tree] = boundary
end
# Check left edge
if tree_nodes[4] in node_set_dict[boundary] && tree_nodes[1] in node_set_dict[boundary]
if tree_nodes[4] in node_set_dict[boundary] &&
tree_nodes[1] in node_set_dict[boundary]
# Left boundary is position 1 in p4est indexing
boundary_names[1, tree] = boundary
end
Expand All @@ -523,38 +527,50 @@ function p4est_mesh_from_standard_abaqus(meshfile, mapping, polydeg,
# https://web.mit.edu/calculix_v2.7/CalculiX/ccx_2.7/doc/ccx/node26.html
for boundary in keys(node_set_dict)
# Check "front face" (y_min)
if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary]
if tree_nodes[1] in node_set_dict[boundary] &&
tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[5] in node_set_dict[boundary] &&
tree_nodes[6] in node_set_dict[boundary]
# Front face is position 3 in p4est indexing
boundary_names[3, tree] = boundary
end
# Check "back face" (y_max)
if tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] &&
tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary]
if tree_nodes[3] in node_set_dict[boundary] &&
tree_nodes[4] in node_set_dict[boundary] &&
tree_nodes[7] in node_set_dict[boundary] &&
tree_nodes[8] in node_set_dict[boundary]
# Front face is position 4 in p4est indexing
boundary_names[4, tree] = boundary
end
# Check "left face" (x_min)
if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary] &&
tree_nodes[5] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary]
if tree_nodes[1] in node_set_dict[boundary] &&
tree_nodes[4] in node_set_dict[boundary] &&
tree_nodes[5] in node_set_dict[boundary] &&
tree_nodes[8] in node_set_dict[boundary]
# Left face is position 1 in p4est indexing
boundary_names[1, tree] = boundary
end
# Check "right face" (x_max)
if tree_nodes[2] in node_set_dict[boundary] && tree_nodes[3] in node_set_dict[boundary] &&
tree_nodes[6] in node_set_dict[boundary] && tree_nodes[7] in node_set_dict[boundary]
if tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[3] in node_set_dict[boundary] &&
tree_nodes[6] in node_set_dict[boundary] &&
tree_nodes[7] in node_set_dict[boundary]
# Right face is position 2 in p4est indexing
boundary_names[2, tree] = boundary
end
# Check "bottom face" (z_min)
if tree_nodes[1] in node_set_dict[boundary] && tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[3] in node_set_dict[boundary] && tree_nodes[4] in node_set_dict[boundary]
if tree_nodes[1] in node_set_dict[boundary] &&
tree_nodes[2] in node_set_dict[boundary] &&
tree_nodes[3] in node_set_dict[boundary] &&
tree_nodes[4] in node_set_dict[boundary]
# Bottom face is position 5 in p4est indexing
boundary_names[5, tree] = boundary
end
# Check "top face" (z_max)
if tree_nodes[5] in node_set_dict[boundary] && tree_nodes[6] in node_set_dict[boundary] &&
tree_nodes[7] in node_set_dict[boundary] && tree_nodes[8] in node_set_dict[boundary]
if tree_nodes[5] in node_set_dict[boundary] &&
tree_nodes[6] in node_set_dict[boundary] &&
tree_nodes[7] in node_set_dict[boundary] &&
tree_nodes[8] in node_set_dict[boundary]
# Top face is position 6 in p4est indexing
boundary_names[6, tree] = boundary
end
Expand All @@ -568,8 +584,13 @@ end

function parse_elements(meshfile, n_trees, n_dims)
@assert n_dims ∈ [2, 3] "Only 2D and 3D meshes are supported"
element_types = n_dims == 2 ? ["*ELEMENT, type=CPS4", "*ELEMENT, type=C2D4", "*ELEMENT, type=S4"] :
["*ELEMENT, type=C3D8"]
element_types = n_dims == 2 ?
[
"*ELEMENT, type=CPS4",
"*ELEMENT, type=C2D4",
"*ELEMENT, type=S4",
] :
["*ELEMENT, type=C3D8"]
# 2D quads: 4 nodes + element index, 3D hexes: 8 nodes + element index
expected_content_length = n_dims == 2 ? 5 : 9

Expand Down Expand Up @@ -621,7 +642,7 @@ function parse_node_sets(meshfile, boundary_symbols)
end
elseif current_symbol !== nothing # Read only if there was already a nodeset specified
# There is always a trailing comma, remove the corresponding empty string
append!(current_nodes, parse.(Int64, split(line, ",")[1:end-1]))
append!(current_nodes, parse.(Int64, split(line, ",")[1:(end - 1)]))
end
end
if current_symbol !== nothing
Expand Down
21 changes: 21 additions & 0 deletions test/test_p4est_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -391,6 +391,27 @@ end
end
end

@trixi_testset "elixir_euler_airfoil_mach2.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_airfoil_mach2.jl"),
l2=[
1.682112049215034e-10, 3.439606229463632e-10,
1.973815936404324e-10, 7.666245265417042e-10,
],
linf=[
2.636251394960709e-8, 4.4272578936244145e-8,
2.8049958258193942e-8, 1.0125492977408612e-7,
],
tspan=(0.0, 0.1))
# Ensure that we do not have excessive memory allocations
# (e.g., from type instabilities)
let
t = sol.t[end]
u_ode = sol.u[end]
du_ode = similar(u_ode)
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
end
end

@trixi_testset "elixir_eulergravity_convergence.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_convergence.jl"),
l2=[
Expand Down
23 changes: 23 additions & 0 deletions test/test_p4est_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -234,6 +234,29 @@ end
end
end

@trixi_testset "elixir_euler_free_stream_boundaries.jl" begin
@test_trixi_include(joinpath(EXAMPLES_DIR,
"elixir_euler_free_stream_boundaries.jl"),
l2=[
3.3824205900551404e-16, 5.761375062854377e-16,
1.0712466724440564e-15, 1.2802903179864802e-15,
2.310819174832704e-15,
],
linf=[
2.1094237467877974e-15, 3.3584246494910985e-15,
9.2148511043888e-15, 1.2101430968414206e-14,
1.9539925233402755e-14,
])
# Ensure that we do not have excessive memory allocations
# (e.g., from type instabilities)
let
t = sol.t[end]
u_ode = sol.u[end]
du_ode = similar(u_ode)
@test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
end
end

@trixi_testset "elixir_euler_free_stream_extruded.jl with HLLC FLux" begin
@test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_extruded.jl"),
l2=[
Expand Down

0 comments on commit 06fb60c

Please sign in to comment.