Skip to content

Commit 5a11404

Browse files
SimonCanSimonCan
committed
Merge branch 'sc/p4est-view-coupled-enhanced' into sc/p4est-view-coupled-hanging-nodes
2 parents d7d67b2 + d38742c commit 5a11404

File tree

59 files changed

+1477
-305
lines changed

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

59 files changed

+1477
-305
lines changed

‎NEWS.md‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,15 @@ used in the Julia ecosystem. Notable changes will be documented in this file
66
for human readability.
77

88

9+
## Changes in the v0.13 lifecycle
10+
11+
#### Added
12+
- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582]).
13+
In the new version, IDP positivity limiting for conservative variables (using
14+
the keyword `positivity_variables_cons` in `SubcellLimiterIDP()`) is supported.
15+
`BoundsCheckCallback` is not supported in 3D yet.
16+
17+
918
## Changes when updating to v0.13 from v0.12.x
1019

1120
#### Changed

‎Project.toml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
name = "Trixi"
22
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
33
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>", "Andrés Rueda-Ramírez <[email protected]>"]
4-
version = "0.13.15-DEV"
4+
version = "0.13.16-DEV"
55

66
[deps]
77
Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

‎docs/Project.toml‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -26,9 +26,6 @@ Trixi = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
2626
Trixi2Vtk = "bc1476a1-1ca6-4cc3-950b-c312b255ff95"
2727
TrixiBase = "9a0f1c46-06d5-4909-a5a3-ce25d3fa3284"
2828

29-
[sources.Trixi]
30-
path = ".."
31-
3229
[compat]
3330
ADTypes = "1.11"
3431
Adapt = "4"
@@ -55,3 +52,6 @@ SparseMatrixColorings = "0.4.21"
5552
Test = "1"
5653
Trixi2Vtk = "0.3.16"
5754
TrixiBase = "0.1.1"
55+
56+
[sources.Trixi]
57+
path = ".."

‎docs/src/conventions.md‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -125,6 +125,16 @@ based on the following rules.
125125
be used when necessary, e.g., to avoid additional overhead when interfacing
126126
with external C libraries such as HDF5, MPI, or visualization.
127127

128+
### Usage of statically sized arrays
129+
130+
Trixi.jl employs statically sized vectors and arrays from the
131+
[StaticArrays.jl](https://github.com/JuliaArrays/StaticArrays.jl) package, i.e.,
132+
`SVector, SMatrix` for the immutable versions and `MVector, MMatrix` for the mutable variants.
133+
These types are preferred for "small" (StaticArrays.jl [recommends as a rule of thumb about 100 entries](https://juliaarrays.github.io/StaticArrays.jl/stable/#When-Static-Arrays-may-be-useful)) arrays with a size known at compile time.
134+
Inspired by this recommendation, Trixi.jl uses the static size `SVector` most notably for numerical fluxes and initial/boundary conditions.
135+
The mutable statically sized `MArray` type is used for **arrays up to two dimensions** which show up in mortars, indicators, and certain volume integrals.
136+
For arrays of higher dimensions (mostly three and four) we use standard Julia `Array` types.
137+
128138
## Numeric types and type stability
129139

130140
In Trixi.jl, we use generic programming to support custom data types to store the numerical simulation data, including standard floating point types and automatic differentiation types.

‎examples/p4est_2d_dgsem/elixir_euler_sedov.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,8 +9,8 @@ equations = CompressibleEulerEquations2D(1.4)
99
"""
1010
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
1111
12-
The Sedov blast wave setup based on Flash
13-
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
12+
The Sedov blast wave setup based on example 35.1.4 from Flash
13+
- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
1414
"""
1515
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
1616
# Set up polar coordinates
@@ -19,7 +19,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq
1919
y_norm = x[2] - inicenter[2]
2020
r = sqrt(x_norm^2 + y_norm^2)
2121

22-
# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
22+
# Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
2323
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
2424
E = 1.0
2525
p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
@@ -43,7 +43,7 @@ initial_condition = initial_condition_sedov_blast_wave
4343
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
4444
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
4545
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
46-
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
46+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
4747
# `StepsizeCallback` (CFL-Condition) and less diffusion.
4848
surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
4949
volume_flux = flux_ranocha

‎examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -8,8 +8,8 @@ equations = CompressibleEulerEquations2D(1.4)
88
"""
99
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
1010
11-
The Sedov blast wave setup based on Flash
12-
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
11+
The Sedov blast wave setup based on example 35.1.4 from Flash
12+
- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
1313
"""
1414
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
1515
# Set up polar coordinates
@@ -18,7 +18,7 @@ function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEq
1818
y_norm = x[2] - inicenter[2]
1919
r = sqrt(x_norm^2 + y_norm^2)
2020

21-
# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
21+
# Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
2222
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
2323
E = 1.0
2424
p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
@@ -42,7 +42,7 @@ initial_condition = initial_condition_sedov_blast_wave
4242
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
4343
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
4444
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
45-
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
45+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
4646
# `StepsizeCallback` (CFL-Condition) and less diffusion.
4747
surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
4848
volume_flux = flux_ranocha

‎examples/p4est_2d_dgsem/elixir_eulermulti_shock.jl‎

Lines changed: 125 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,125 @@
1+
using OrdinaryDiffEqSSPRK, OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
# 1) Dry Air 2) SF_6
5+
equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
6+
gas_constants = (0.287, 1.578))
7+
8+
"""
9+
initial_condition_shock(coordinates, t, equations::CompressibleEulerEquations2D)
10+
11+
Shock traveling from left to right where it interacts with a Perturbed interface.
12+
"""
13+
@inline function initial_condition_shock(x, t,
14+
equations::CompressibleEulerMulticomponentEquations2D)
15+
rho_0 = 1.25 # kg/m^3
16+
p_0 = 101325 # Pa
17+
T_0 = 293 # K
18+
u_0 = 352 # m/s
19+
d = 20
20+
w = 40
21+
22+
if x[1] < 25
23+
# Shock region.
24+
v1 = 0.35
25+
v2 = 0.0
26+
rho1 = 1.72
27+
rho2 = 0.03
28+
p = 1.57
29+
elseif (x[1] <= 30) || (x[1] <= 70 && abs(125 - x[2]) > w / 2)
30+
# Intermediate region.
31+
v1 = 0.0
32+
v2 = 0.0
33+
rho1 = 1.25
34+
rho2 = 0.03
35+
p = 1.0
36+
else
37+
(x[1] <= 70 + d)
38+
# SF_6 region.
39+
v1 = 0.0
40+
v2 = 0.0
41+
rho1 = 0.03
42+
rho2 = 6.03 #SF_6
43+
p = 1.0
44+
end
45+
46+
return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
47+
end
48+
49+
# Define the simulation.
50+
initial_condition = initial_condition_shock
51+
52+
surface_flux = flux_lax_friedrichs
53+
volume_flux = flux_ranocha
54+
basis = LobattoLegendreBasis(3)
55+
56+
limiter_idp = SubcellLimiterIDP(equations, basis;
57+
positivity_variables_cons = ["rho" * string(i)
58+
for i in eachcomponent(equations)])
59+
60+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
61+
volume_flux_dg = volume_flux,
62+
volume_flux_fv = surface_flux)
63+
64+
solver = DGSEM(basis, surface_flux, volume_integral)
65+
66+
coordinates_min = (0.0, 0.0)
67+
coordinates_max = (250.0, 250.0)
68+
mesh = P4estMesh((32, 32), polydeg = 3, coordinates_min = (0.0, 0.0),
69+
coordinates_max = (250.0, 250.0), periodicity = (false, true),
70+
initial_refinement_level = 0)
71+
72+
# Completely free outflow
73+
function boundary_condition_outflow(u_inner, normal_direction::AbstractVector,
74+
x, t,
75+
surface_flux_function,
76+
equations)
77+
# Calculate the boundary flux entirely from the internal solution state
78+
return flux(u_inner, normal_direction, equations)
79+
end
80+
81+
boundary_conditions = Dict(:x_neg => BoundaryConditionDirichlet(initial_condition),
82+
:x_pos => boundary_condition_outflow)
83+
84+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
85+
boundary_conditions = boundary_conditions)
86+
87+
###############################################################################
88+
# ODE solvers, callbacks etc.
89+
90+
tspan = (0.0, 5.0)
91+
ode = semidiscretize(semi, tspan)
92+
93+
summary_callback = SummaryCallback()
94+
95+
analysis_interval = 10
96+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
97+
save_analysis = true)
98+
99+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
100+
101+
save_solution = SaveSolutionCallback(dt = 2.0,
102+
save_initial_solution = true,
103+
save_final_solution = true,
104+
solution_variables = cons2prim)
105+
106+
stepsize_callback = StepsizeCallback(cfl = 0.2)
107+
108+
callbacks = CallbackSet(summary_callback,
109+
analysis_callback,
110+
alive_callback,
111+
save_solution,
112+
stepsize_callback)
113+
114+
###############################################################################
115+
# run the simulation
116+
117+
stage_callbacks = (SubcellLimiterIDPCorrection(),
118+
BoundsCheckCallback(save_errors = false, interval = 100))
119+
# `interval` is used when calling this elixir in the tests with `save_errors=true`.
120+
121+
@time begin
122+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
123+
dt = 0.1, # solve needs some value here but it will be overwritten by the stepsize_callback
124+
ode_default_options()..., callback = callbacks)
125+
end

‎examples/p4est_3d_dgsem/elixir_euler_sedov.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,8 +9,8 @@ equations = CompressibleEulerEquations3D(1.4)
99
"""
1010
initial_condition_medium_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
1111
12-
The Sedov blast wave setup based on Flash
13-
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
12+
The Sedov blast wave setup based on example 35.1.4 from Flash
13+
- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
1414
with smaller strength of the initial discontinuity.
1515
"""
1616
function initial_condition_medium_sedov_blast_wave(x, t,
@@ -22,7 +22,7 @@ function initial_condition_medium_sedov_blast_wave(x, t,
2222
z_norm = x[3] - inicenter[3]
2323
r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
2424

25-
# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
25+
# Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
2626
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
2727
E = 1.0
2828
p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2)
@@ -45,7 +45,7 @@ initial_condition = initial_condition_medium_sedov_blast_wave
4545
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
4646
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
4747
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
48-
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
48+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
4949
# `StepsizeCallback` (CFL-Condition) and less diffusion.
5050
surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
5151
volume_flux = flux_ranocha

‎examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl‎

Lines changed: 98 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,98 @@
1+
using Trixi
2+
3+
###############################################################################
4+
# semidiscretization of the compressible Euler equations
5+
6+
equations = CompressibleEulerEquations3D(1.4)
7+
8+
"""
9+
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
10+
11+
The Sedov blast wave setup based on example 35.1.4 from Flash
12+
- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
13+
with smaller strength of the initial discontinuity.
14+
"""
15+
function initial_condition_sedov_blast_wave(x, t,
16+
equations::CompressibleEulerEquations3D)
17+
# Set up polar coordinates
18+
inicenter = SVector(0.0, 0.0, 0.0)
19+
x_norm = x[1] - inicenter[1]
20+
y_norm = x[2] - inicenter[2]
21+
z_norm = x[3] - inicenter[3]
22+
r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
23+
24+
# Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
25+
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
26+
E = 1.0
27+
p0_inner = 3 * (equations.gamma - 1) * E / (4 * pi * r0^2)
28+
p0_outer = 1.0e-1 # "simpler" setup since positivity limiter for pressure is not yet supported in 3D
29+
30+
# Calculate primitive variables
31+
rho = 1.0
32+
v1 = 0.0
33+
v2 = 0.0
34+
v3 = 0.0
35+
p = r > r0 ? p0_outer : p0_inner
36+
37+
return prim2cons(SVector(rho, v1, v2, v3, p), equations)
38+
end
39+
40+
initial_condition = initial_condition_sedov_blast_wave
41+
42+
surface_flux = flux_lax_friedrichs
43+
volume_flux = flux_ranocha
44+
polydeg = 3
45+
basis = LobattoLegendreBasis(polydeg)
46+
limiter_idp = SubcellLimiterIDP(equations, basis;
47+
positivity_variables_cons = ["rho"])
48+
volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
49+
volume_flux_dg = volume_flux,
50+
volume_flux_fv = surface_flux)
51+
solver = DGSEM(basis, surface_flux, volume_integral)
52+
53+
coordinates_min = (-1.0, -1.0, -1.0)
54+
coordinates_max = (1.0, 1.0, 1.0)
55+
56+
trees_per_dimension = (8, 8, 8)
57+
mesh = P4estMesh(trees_per_dimension,
58+
polydeg = 1, initial_refinement_level = 0,
59+
coordinates_min = coordinates_min, coordinates_max = coordinates_max,
60+
periodicity = true)
61+
62+
# create the semi discretization object
63+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
64+
65+
###############################################################################
66+
# ODE solvers, callbacks etc.
67+
68+
tspan = (0.0, 3.0)
69+
ode = semidiscretize(semi, tspan)
70+
71+
summary_callback = SummaryCallback()
72+
73+
analysis_interval = 100
74+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
75+
76+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
77+
78+
save_solution = SaveSolutionCallback(interval = 10,
79+
save_initial_solution = true,
80+
save_final_solution = true,
81+
extra_node_variables = (:limiting_coefficient,))
82+
83+
stepsize_callback = StepsizeCallback(cfl = 0.5)
84+
85+
callbacks = CallbackSet(summary_callback,
86+
analysis_callback,
87+
alive_callback,
88+
save_solution,
89+
stepsize_callback)
90+
91+
###############################################################################
92+
# run the simulation
93+
94+
stage_callbacks = (SubcellLimiterIDPCorrection(),)
95+
96+
sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);
97+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
98+
callback = callbacks);

‎examples/structured_1d_dgsem/elixir_euler_sedov.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -9,16 +9,16 @@ equations = CompressibleEulerEquations1D(1.4)
99
"""
1010
initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
1111
12-
The Sedov blast wave setup based on Flash
13-
- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
12+
The Sedov blast wave setup based on example 35.1.4 from Flash
13+
- https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
1414
"""
1515
function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations1D)
1616
# Set up polar coordinates
1717
inicenter = SVector(0.0)
1818
x_norm = x[1] - inicenter[1]
1919
r = abs(x_norm)
2020

21-
# Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
21+
# Setup based on example 35.1.4 in https://flash.rochester.edu/site/flashcode/user_support/flash4_ug_4p8.pdf
2222
r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
2323
# r0 = 0.5 # = more reasonable setup
2424
E = 1.0
@@ -40,7 +40,7 @@ initial_condition = initial_condition_sedov_blast_wave
4040
# In the `StepsizeCallback`, though, the less diffusive `max_abs_speeds` is employed which is consistent with `max_abs_speed`.
4141
# Thus, we exchanged in PR#2458 the default wave speed used in the LLF flux to `max_abs_speed`.
4242
# To ensure that every example still runs we specify explicitly `FluxLaxFriedrichs(max_abs_speed_naive)`.
43-
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
43+
# We remark, however, that the now default `max_abs_speed` is in general recommended due to compliance with the
4444
# `StepsizeCallback` (CFL-Condition) and less diffusion.
4545
surface_flux = FluxLaxFriedrichs(max_abs_speed_naive)
4646
volume_flux = flux_ranocha

0 commit comments

Comments
 (0)