Skip to content

Commit

Permalink
Merge branch 'main' into helper
Browse files Browse the repository at this point in the history
  • Loading branch information
huiyuxie authored May 17, 2024
2 parents d7716cb + b98b6f0 commit fa2c9d8
Show file tree
Hide file tree
Showing 17 changed files with 330 additions and 82 deletions.
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "Trixi"
uuid = "a7f1ee26-1774-49b1-8366-f1abc58fbfcb"
authors = ["Michael Schlottke-Lakemper <[email protected]>", "Gregor Gassner <[email protected]>", "Hendrik Ranocha <[email protected]>", "Andrew R. Winters <[email protected]>", "Jesse Chan <[email protected]>"]
version = "0.7.13-pre"
version = "0.7.14-pre"

[deps]
CodeTracking = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
Expand Down
19 changes: 12 additions & 7 deletions docs/literate/src/files/subcell_shock_capturing.jl
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,8 @@ positivity_variables_nonlinear = [pressure]

# ### Local bounds
# Second, Trixi.jl supports the limiting with local bounds for conservative variables using a
# two-sided Zalesak-type limiter ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2)).
# two-sided Zalesak-type limiter ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2))
# and for general non-linear variables using a one-sided Newton-bisection algorithm.
# They allow to avoid spurious oscillations within the global bounds and to improve the
# shock-capturing capabilities of the method. The corresponding numerical admissibility conditions
# are frequently formulated as local maximum or minimum principles. The local bounds are computed
Expand All @@ -108,6 +109,11 @@ positivity_variables_nonlinear = [pressure]
# the following.
local_twosided_variables_cons = ["rho"]

# To limit non-linear variables locally, pass the variable function combined with the requested
# bound (`min` or `max`) as a tuple. For instance, to impose a lower local bound on the modified
# specific entropy [`Trixi.entropy_guermond_etal`](@ref), use
local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]

# ## Exemplary simulation
# How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary
# setup. This will be a simplified version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`.
Expand Down Expand Up @@ -277,11 +283,10 @@ plot(sol)
# timestep and simulation time.
# ````
# iter, simu_time, rho_min, rho_max
# 1, 0.0, 0.0, 0.0
# 101, 0.29394033217556337, 0.0, 0.0
# 201, 0.6012597465597065, 0.0, 0.0
# 301, 0.9559096690030839, 0.0, 0.0
# 401, 1.3674274981949077, 0.0, 0.0
# 501, 1.8395301696603052, 0.0, 0.0
# 100, 0.29103427131404924, 0.0, 0.0
# 200, 0.5980281923063808, 0.0, 0.0
# 300, 0.9520853560765293, 0.0, 0.0
# 400, 1.3630295622683186, 0.0, 0.0
# 500, 1.8344999624013498, 0.0, 0.0
# 532, 1.9974179806990118, 0.0, 0.0
# ````
71 changes: 67 additions & 4 deletions docs/src/conventions.md
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ following naming conventions:
use these indices.


# Keywords in elixirs
## Keywords in elixirs

Trixi.jl is distributed with several examples in the form of elixirs, small
Julia scripts containing everything to set up and run a simulation. Working
Expand All @@ -61,9 +61,9 @@ can only perform simple replacements. Some standard variables names are
Moreover, [`convergence_test`](@ref) requires that the spatial resolution is
set via the keywords
- `initial_refinement_level`
(an integer, e.g. for the [`TreeMesh`](@ref) and the [`P4estMesh`](@ref)) or
(an integer, e.g., for the [`TreeMesh`](@ref) and the [`P4estMesh`](@ref)) or
- `cells_per_dimension`
(a tuple of integers, one per spatial dimension, e.g. for the [`StructuredMesh`](@ref)
(a tuple of integers, one per spatial dimension, e.g., for the [`StructuredMesh`](@ref)
and the [`DGMultiMesh`](@ref)).


Expand Down Expand Up @@ -101,8 +101,71 @@ based on the following rules.
(or `wrap_array_native(u_ode, semi)`) for further processing.
- When some solution is passed together with the `mesh, equations, solver, cache, ...`,
it is already wrapped via `wrap_array` (or `wrap_array_native`).
- Exceptions of this rule are possible, e.g. for AMR, but must be documented in
- Exceptions of this rule are possible, e.g., for AMR, but must be documented in
the code.
- `wrap_array` should be used as default option. `wrap_array_native` should only
be used when necessary, e.g., to avoid additional overhead when interfacing
with external C libraries such as HDF5, MPI, or visualization.

## Numeric types and type stability

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.
Specifically, `Float32` and `Float64` types are fully supported, including the ability to run Trixi.jl on hardware that only supports `Float32` types.
We ensure the type stability of these numeric types throughout the development process.
Below are some guidelines to apply in various scenarios.

### Exact floating-point numbers

Some real numbers can be represented exactly as both `Float64` and `Float32` values (e.g., `0.25`, `0.5`, `1/2`). We prefer to use `Float32` type for these numbers to achieve a concise way of possible type promotion. For example,
```julia
# Assume we have `0.25`, `0.5`, `1/2` in function
0.25f0, 0.5f0, 0.5f0 # corresponding numbers
```
Generally, this equivalence is true for integer multiples of powers of two. That is, numbers that can be written as ``m 2^n``, where ``m, n \in \mathbb{Z}``, and where ``m`` and ``n`` are such that the result is representable as a [single precision floating point](https://en.wikipedia.org/wiki/Single-precision_floating-point_format) value. If a decimal value `v` is exactly representable in `Float32`, the expression
```julia
Float32(v) == v
```
will evaluate to `true`.

### Non-exact floating-point numbers

For real numbers that cannot be exactly represented in machine precision (e.g., `0.1`, `1/3`, `pi`), use the `convert` function to make them consistent with the type of the function input. For example,
```julia
# Assume we are handling `pi` in function
function foo(..., input, ...)
RealT = eltype(input) # see **notes** below
# ...
c1 = convert(RealT, pi) * c2 # sample operation
# ...
end
```

### Integer numbers

Integers need special consideration. Using functions like `convert` (as mentioned above), `zero`, and `one` to change integers to a specific type or keeping them in integer format is feasible in most cases. We aim to balance code readability and consistency while maintaining type stability. Here are some examples to guide our developers.
```julia
# The first example - `SVector`, keep code consistency within `SVector`
SVector(0, 0, 0)
SVector(zero(RealT), zero(RealT), zero(RealT))
Svector(one(RealT), one(RealT), one(RealT))

# The second example - inner functions, keep them type-stable as well
function foo(..., input, ...)
RealT = eltype(input) # see **notes** below
# ...
c1 = c2 > 0.5f0 ? one(RealT) : convert(RealT, 0.1) # make type-stable
# ...
end

# The third example - some operations (e.g., `/`, `sqrt`, `inv`), convert them definitely
c1 = convert(RealT, 4) # suppose we get RealT before
c2 = 1 / c1
c3 = sqrt(c1)
c4 = inv(c1)
```
In general, in the case of integer numbers, our developers should apply a case-by-case strategy to maintain type stability.

### Notes
1. If the function gets a local pointwise vector of the solution variables `u` such as `flux(u, equations)`, use `u` to determine the real type `eltype(u)`.
2. If `u` is not passed as an argument but a vector of coordinates `x` such as `initial_condition(x, t, equations)`, use `eltype(x)` instead.
3. Choose an appropriate argument to determine the real type otherwise.
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
using Trixi
using OrdinaryDiffEq

gamma_gas = 1.4
equations = CompressibleEulerEquations1D(gamma_gas)

###############################################################################
# setup the GSBP DG discretization that uses the Gauss operators from
# Chan, Del Rey Fernandez, Carpenter (2019).
# [https://doi.org/10.1137/18M1209234](https://doi.org/10.1137/18M1209234)

# Shu-Osher initial condition for 1D compressible Euler equations
# Example 8 from Shu, Osher (1989).
# [https://doi.org/10.1016/0021-9991(89)90222-2](https://doi.org/10.1016/0021-9991(89)90222-2)
function initial_condition_shu_osher(x, t, equations::CompressibleEulerEquations1D)
x0 = -4

rho_left = 27 / 7
v_left = 4 * sqrt(35) / 9
p_left = 31 / 3

# Replaced v_right = 0 to v_right = 0.1 to avoid positivity issues.
v_right = 0.1
p_right = 1.0

rho = ifelse(x[1] > x0, 1 + 1 / 5 * sin(5 * x[1]), rho_left)
v = ifelse(x[1] > x0, v_right, v_left)
p = ifelse(x[1] > x0, p_right, p_left)

return prim2cons(SVector(rho, v, p),
equations)
end

initial_condition = initial_condition_shu_osher

surface_flux = flux_lax_friedrichs
volume_flux = flux_ranocha

polydeg = 3
basis = DGMultiBasis(Line(), polydeg, approximation_type = GaussSBP())

indicator_sc = IndicatorHennemannGassner(equations, basis,
alpha_max = 0.5,
alpha_min = 0.001,
alpha_smooth = true,
variable = density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
volume_flux_dg = volume_flux,
volume_flux_fv = surface_flux)

dg = DGMulti(basis,
surface_integral = SurfaceIntegralWeakForm(surface_flux),
volume_integral = volume_integral)

boundary_condition = BoundaryConditionDirichlet(initial_condition)
boundary_conditions = (; :entire_boundary => boundary_condition)

###############################################################################
# setup the 1D mesh

cells_per_dimension = (64,)
mesh = DGMultiMesh(dg, cells_per_dimension,
coordinates_min = (-5.0,), coordinates_max = (5.0,),
periodicity = false)

###############################################################################
# setup the semidiscretization and ODE problem

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
dg, boundary_conditions = boundary_conditions)

tspan = (0.0, 1.0)
ode = semidiscretize(semi, tspan)

###############################################################################
# setup the callbacks

# prints a summary of the simulation setup and resets the timers
summary_callback = SummaryCallback()

# analyse the solution in regular intervals and prints the results
analysis_callback = AnalysisCallback(semi, interval = 100, uEltype = real(dg))

# handles the re-calculation of the maximum Δt after each time step
stepsize_callback = StepsizeCallback(cfl = 0.1)

# collect all callbacks such that they can be passed to the ODE solver
callbacks = CallbackSet(summary_callback, analysis_callback, stepsize_callback)

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

sol = solve(ode, SSPRK43(), adaptive = true, callback = callbacks, save_everystep = false)
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,9 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
stage_callbacks = (SubcellLimiterIDPCorrection(),
BoundsCheckCallback(save_errors = false, interval = 100))
# `interval` is used when calling this elixir in the tests with `save_errors=true`.

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
Original file line number Diff line number Diff line change
Expand Up @@ -85,7 +85,9 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
stage_callbacks = (SubcellLimiterIDPCorrection(),
BoundsCheckCallback(save_errors = false, interval = 100))
# `interval` is used when calling this elixir in the tests with `save_errors=true`.

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
Original file line number Diff line number Diff line change
Expand Up @@ -133,10 +133,9 @@ callbacks = CallbackSet(summary_callback,
###############################################################################
# run the simulation

output_directory = "out"
stage_callbacks = (SubcellLimiterIDPCorrection(),
BoundsCheckCallback(save_errors = true, interval = 100,
output_directory = output_directory))
BoundsCheckCallback(save_errors = false, interval = 100))
# `interval` is used when calling this elixir in the tests with `save_errors=true`.

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
11 changes: 8 additions & 3 deletions src/callbacks_stage/positivity_zhang_shu_dg2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
function limiter_zhang_shu!(u, threshold::Real, variable,
mesh::AbstractMesh{2}, equations, dg::DGSEM, cache)
@unpack weights = dg.basis
@unpack inverse_jacobian = cache.elements

@threaded for element in eachelement(dg, cache)
# determine minimum value
Expand All @@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable,

# compute mean value
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, element))
total_volume = zero(eltype(u))
for j in eachnode(dg), i in eachnode(dg)
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
i, j, element)))
u_node = get_node_vars(u, equations, dg, i, j, element)
u_mean += u_node * weights[i] * weights[j]
u_mean += u_node * weights[i] * weights[j] * volume_jacobian
total_volume += weights[i] * weights[j] * volume_jacobian
end
# note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
u_mean = u_mean / 2^ndims(mesh)
# normalize with the total volume
u_mean = u_mean / total_volume

# We compute the value directly with the mean values, as we assume that
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).
Expand Down
11 changes: 8 additions & 3 deletions src/callbacks_stage/positivity_zhang_shu_dg3d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
function limiter_zhang_shu!(u, threshold::Real, variable,
mesh::AbstractMesh{3}, equations, dg::DGSEM, cache)
@unpack weights = dg.basis
@unpack inverse_jacobian = cache.elements

@threaded for element in eachelement(dg, cache)
# determine minimum value
Expand All @@ -22,12 +23,16 @@ function limiter_zhang_shu!(u, threshold::Real, variable,

# compute mean value
u_mean = zero(get_node_vars(u, equations, dg, 1, 1, 1, element))
total_volume = zero(eltype(u))
for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
volume_jacobian = abs(inv(get_inverse_jacobian(inverse_jacobian, mesh,
i, j, k, element)))
u_node = get_node_vars(u, equations, dg, i, j, k, element)
u_mean += u_node * weights[i] * weights[j] * weights[k]
u_mean += u_node * weights[i] * weights[j] * weights[k] * volume_jacobian
total_volume += weights[i] * weights[j] * weights[k] * volume_jacobian
end
# note that the reference element is [-1,1]^ndims(dg), thus the weights sum to 2
u_mean = u_mean / 2^ndims(mesh)
# normalize with the total volume
u_mean = u_mean / total_volume

# We compute the value directly with the mean values, as we assume that
# Jensen's inequality holds (e.g. pressure for compressible Euler equations).
Expand Down
38 changes: 19 additions & 19 deletions src/callbacks_stage/subcell_bounds_check.jl
Original file line number Diff line number Diff line change
Expand Up @@ -38,35 +38,37 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
(; t, iter, alg) = integrator
u = wrap_array(u_ode, mesh, equations, solver, cache)

@trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
solver.volume_integral)

save_errors = callback.save_errors && (callback.interval > 0) &&
(stage == length(alg.c)) &&
(iter % callback.interval == 0 || integrator.finalstep)
@trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
solver.volume_integral, t,
iter + 1,
callback.output_directory,
save_errors)
((iter + 1) % callback.interval == 0 || # Every `interval` time steps
integrator.finalstep || # Planned last time step
(iter + 1) >= integrator.opts.maxiters) # Maximum iterations reached
if save_errors
@trixi_timeit timer() "save_errors" save_bounds_check_errors(callback.output_directory,
u, t, iter + 1,
equations,
solver.volume_integral)
end
end

function check_bounds(u, mesh, equations, solver, cache,
volume_integral::AbstractVolumeIntegral, t, iter,
output_directory, save_errors)
return nothing
@inline function check_bounds(u, mesh, equations, solver, cache,
volume_integral::VolumeIntegralSubcellLimiting)
check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter)
end

function check_bounds(u, mesh, equations, solver, cache,
volume_integral::VolumeIntegralSubcellLimiting, t, iter,
output_directory, save_errors)
check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter, t, iter,
output_directory, save_errors)
@inline function save_bounds_check_errors(output_directory, u, t, iter, equations,
volume_integral::VolumeIntegralSubcellLimiting)
save_bounds_check_errors(output_directory, u, t, iter, equations,
volume_integral.limiter)
end

function init_callback(callback::BoundsCheckCallback, semi)
init_callback(callback, semi, semi.solver.volume_integral)
end

init_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing

function init_callback(callback::BoundsCheckCallback, semi,
volume_integral::VolumeIntegralSubcellLimiting)
init_callback(callback, semi, volume_integral.limiter)
Expand Down Expand Up @@ -116,8 +118,6 @@ function finalize_callback(callback::BoundsCheckCallback, semi)
finalize_callback(callback, semi, semi.solver.volume_integral)
end

finalize_callback(callback::BoundsCheckCallback, semi, volume_integral::AbstractVolumeIntegral) = nothing

function finalize_callback(callback::BoundsCheckCallback, semi,
volume_integral::VolumeIntegralSubcellLimiting)
finalize_callback(callback, semi, volume_integral.limiter)
Expand Down
Loading

0 comments on commit fa2c9d8

Please sign in to comment.