20: Custom semidiscretizations
As described in the overview section, semidiscretizations are high-level descriptions of spatial discretizations in Trixi.jl. Trixi.jl's main focus is on hyperbolic conservation laws represented in a SemidiscretizationHyperbolic. Hyperbolic-parabolic problems based on the advection-diffusion equation or the compressible Navier-Stokes equations can be represented in a SemidiscretizationHyperbolicParabolic. This is described in the basic tutorial on parabolic terms and its extension to custom parabolic terms. In this tutorial, we will describe how these semidiscretizations work and how they can be used to create custom semidiscretizations involving also other tasks.
Overview of the right-hand side evaluation
The semidiscretizations provided by Trixi.jl are set up to create ODEProblems from the SciML ecosystem for ordinary differential equations. In particular, a spatial semidiscretization can be wrapped in an ODE problem using semidiscretize, which returns an ODEProblem. This ODEProblem bundles an initial condition, a right-hand side (RHS) function, the time span, and possible parameters. The ODEProblems created by Trixi.jl use the semidiscretization passed to semidiscretize as a parameter. For a SemidiscretizationHyperbolic, the ODEProblem wraps Trixi.rhs! as ODE RHS. For a SemidiscretizationHyperbolicParabolic, Trixi.jl uses a SplitODEProblem combining Trixi.rhs_parabolic! for the (potentially) stiff part and Trixi.rhs! for the other part.
Standard Trixi.jl setup
In this tutorial, we will consider the linear advection equation with source term
\[\partial_t u(t,x) + \partial_x u(t,x) = -\exp(-t) \sin\bigl(\pi (x - t) \bigr)\]
with periodic boundary conditions in the domain [-1, 1] as a model problem. The initial condition is
\[u(0,x) = \sin(\pi x).\]
The source term results in some damping and the analytical solution
\[u(t,x) = \exp(-t) \sin\bigl(\pi (x - t) \bigr).\]
First, we discretize this equation using the standard functionality of Trixi.jl.
using Trixi, OrdinaryDiffEq, PlotsThe linear scalar advection equation is already implemented in Trixi.jl as LinearScalarAdvectionEquation1D. We construct it with an advection velocity 1.0.
equations = LinearScalarAdvectionEquation1D(1.0)┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation1D │
│ ═══════════════════════════════ │
│ #variables: ………………………………………………… 1 │
│ │ variable 1: …………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘Next, we use a standard DGSEM solver.
solver = DGSEM(polydeg = 3)┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: ……………………………………… flux_central │
│ volume integral: …………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘We create a simple TreeMesh in 1D.
coordinates_min = (-1.0,)
coordinates_max = (+1.0,)
mesh = TreeMesh(coordinates_min, coordinates_max;
initial_refinement_level = 4,
n_cells_max = 10^4,
periodicity = true)┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{1, Trixi.SerialTree{1}} │
│ ════════════════════════════════ │
│ center: …………………………………………………………… [0.0] │
│ length: …………………………………………………………… 2.0 │
│ periodicity: ……………………………………………… (true,) │
│ current #cells: ……………………………………… 31 │
│ #leaf-cells: ……………………………………………… 16 │
│ maximum #cells: ……………………………………… 10000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘We wrap everything in in a semidiscretization and pass the source terms as a standard Julia function. Please note that Trixi.jl uses SVectors from StaticArrays.jl to store the conserved variables u. Thus, the return value of the source terms must be wrapped in an SVector - even if we consider just a scalar problem.
function initial_condition(x, t, equations)
return SVector(exp(-t) * sinpi(x[1] - t))
end
function source_terms_standard(u, x, t, equations)
return -initial_condition(x, t, equations)
end
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition,
solver;
source_terms = source_terms_standard)┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ………………………… 1 │
│ mesh: ………………………………………………………………… TreeMesh{1, Trixi.SerialTree{1}} with length 31 │
│ equations: …………………………………………………… LinearScalarAdvectionEquation1D │
│ initial condition: ……………………………… initial_condition │
│ boundary conditions: ………………………… Trixi.BoundaryConditionPeriodic │
│ source terms: …………………………………………… source_terms_standard │
│ solver: …………………………………………………………… DG │
│ total #DOFs per field: …………………… 64 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘Now, we can create the ODEProblem, solve the resulting ODE using a time integration method from OrdinaryDiffEq.jl, and visualize the numerical solution at the final time using Plots.jl.
tspan = (0.0, 3.0)
ode = semidiscretize(semi, tspan)
sol = solve(ode, RDPK3SpFSAL49(); ode_default_options()...)
plot(sol; label = "numerical sol.", legend = :topright)
We can also plot the analytical solution for comparison. Since Trixi.jl uses SVectors for the variables, we take their first (and only) component to get the scalar value for manual plotting.
let
x = range(-1.0, 1.0; length = 200)
plot!(x, first.(initial_condition.(x, sol.t[end], equations)),
label = "analytical sol.", linestyle = :dash, legend = :topright)
end
We can also add the initial condition to the plot.
plot!(sol.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
You can of course also use some callbacks provided by Trixi.jl as usual.
summary_callback = SummaryCallback()
analysis_interval = 100
analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
alive_callback = AliveCallback(; analysis_interval)
callbacks = CallbackSet(summary_callback,
analysis_callback,
alive_callback)
sol = solve(ode, RDPK3SpFSAL49();
ode_default_options()..., callback = callbacks)
summary_callback()
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ………………………… 1 │
│ mesh: ………………………………………………………………… TreeMesh{1, Trixi.SerialTree{1}} with length 31 │
│ equations: …………………………………………………… LinearScalarAdvectionEquation1D │
│ initial condition: ……………………………… initial_condition │
│ boundary conditions: ………………………… Trixi.BoundaryConditionPeriodic │
│ source terms: …………………………………………… source_terms_standard │
│ solver: …………………………………………………………… DG │
│ total #DOFs per field: …………………… 64 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{1, Trixi.SerialTree{1}} │
│ ════════════════════════════════ │
│ center: …………………………………………………………… [0.0] │
│ length: …………………………………………………………… 2.0 │
│ periodicity: ……………………………………………… (true,) │
│ current #cells: ……………………………………… 31 │
│ #leaf-cells: ……………………………………………… 16 │
│ maximum #cells: ……………………………………… 10000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation1D │
│ ═══════════════════════════════ │
│ #variables: ………………………………………………… 1 │
│ │ variable 1: …………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: ……………………………………… flux_central │
│ volume integral: …………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback │
│ ════════════════ │
│ interval: ……………………………………………………… 100 │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6) │
│ │ error 1: …………………………………………………… l2_error │
│ │ error 2: …………………………………………………… linf_error │
│ │ integral 1: …………………………………………… entropy_timederivative │
│ save analysis to file: …………………… no │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AliveCallback │
│ ═════════════ │
│ interval: ……………………………………………………… 10 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ………………………………………………… 0.0 │
│ Final time: ………………………………………………… 3.0 │
│ time integrator: …………………………………… RDPK3SpFSAL49 │
│ adaptive: ……………………………………………………… true │
│ abstol: …………………………………………………………… 1.0e-6 │
│ reltol: …………………………………………………………… 0.001 │
│ controller: ………………………………………………… PIDController(beta=[0.38, -0.18,…iter=default_dt_factor_limiter) │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: ……………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 6.82000000e-07 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 3.59991862e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 3550.696 MiB
#elements: 16
Variable: scalar
L2 error: 5.57368408e-06
Linf error: 1.21294882e-05
∑∂S/∂U ⋅ Uₜ : -5.00000000e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 9.9686e-02 │ sim. time: 3.4078e-01 (11.359%) │ run time: 2.3719e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 4.9618e-04 s
#timesteps: 30 │ Δt: 6.9464e-02 │ sim. time: 1.9556e+00 (65.188%) │ run time: 7.2470e-04 s
#timesteps: 40 │ Δt: 8.5031e-02 │ sim. time: 2.7005e+00 (90.016%) │ run time: 9.5178e-04 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 1.31045100e-03 s
Δt: 4.31148548e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.59317060e-08 s
PID: 4.01605216e-08 s
#DOFs per field: 64 alloc'd memory: 3550.737 MiB
#elements: 16
Variable: scalar
L2 error: 3.33421293e-05
Linf error: 1.31835569e-04
∑∂S/∂U ⋅ Uₜ : -1.23948350e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 2.02ms / 67.1% 98.0KiB / 48.2%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────
rhs! 417 920μs 67.9% 2.21μs 6.61KiB 14.0% 16.2B
source terms 417 387μs 28.6% 928ns 0.00B 0.0% 0.00B
~rhs!~ 417 275μs 20.3% 659ns 6.61KiB 14.0% 16.2B
volume integral 417 114μs 8.4% 274ns 0.00B 0.0% 0.00B
interface flux 417 44.8μs 3.3% 107ns 0.00B 0.0% 0.00B
prolong2interfaces 417 24.5μs 1.8% 58.7ns 0.00B 0.0% 0.00B
Jacobian 417 17.6μs 1.3% 42.2ns 0.00B 0.0% 0.00B
surface integral 417 17.6μs 1.3% 42.2ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 14.4μs 1.1% 34.6ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.2μs 1.0% 31.5ns 0.00B 0.0% 0.00B
boundary flux 417 12.1μs 0.9% 29.1ns 0.00B 0.0% 0.00B
analyze solution 2 435μs 32.1% 217μs 40.6KiB 86.0% 20.3KiB
────────────────────────────────────────────────────────────────────────────────Using a custom ODE right-hand side function
Next, we will solve the same problem but use our own ODE RHS function. To demonstrate this, we will artificially create a global variable containing the current time of the simulation.
const GLOBAL_TIME = Ref(0.0)
function source_terms_custom(u, x, t, equations)
t = GLOBAL_TIME[]
return -initial_condition(x, t, equations)
endsource_terms_custom (generic function with 1 method)Next, we create our own RHS function to update the global time of the simulation before calling the RHS function from Trixi.jl.
function rhs_source_custom!(du_ode, u_ode, semi, t)
GLOBAL_TIME[] = t
Trixi.rhs!(du_ode, u_ode, semi, t)
endrhs_source_custom! (generic function with 1 method)Next, we create an ODEProblem manually copying over the data from the one we got from semidiscretize earlier.
ode_source_custom = ODEProblem(rhs_source_custom!,
ode.u0,
ode.tspan,
ode.p #= semi =#)
sol_source_custom = solve(ode_source_custom, RDPK3SpFSAL49();
ode_default_options()...)
plot(sol_source_custom; label = "numerical sol.")
let
x = range(-1.0, 1.0; length = 200)
plot!(x, first.(initial_condition.(x, sol_source_custom.t[end], equations)),
label = "analytical sol.", linestyle = :dash, legend = :topleft)
end
plot!(sol_source_custom.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
This also works with callbacks as usual.
summary_callback = SummaryCallback()
analysis_interval = 100
analysis_callback = AnalysisCallback(semi; interval = analysis_interval)
alive_callback = AliveCallback(; analysis_interval)
callbacks = CallbackSet(summary_callback,
analysis_callback,
alive_callback)
sol = solve(ode_source_custom, RDPK3SpFSAL49();
ode_default_options()..., callback = callbacks)
summary_callback()
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ………………………… 1 │
│ mesh: ………………………………………………………………… TreeMesh{1, Trixi.SerialTree{1}} with length 31 │
│ equations: …………………………………………………… LinearScalarAdvectionEquation1D │
│ initial condition: ……………………………… initial_condition │
│ boundary conditions: ………………………… Trixi.BoundaryConditionPeriodic │
│ source terms: …………………………………………… source_terms_standard │
│ solver: …………………………………………………………… DG │
│ total #DOFs per field: …………………… 64 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{1, Trixi.SerialTree{1}} │
│ ════════════════════════════════ │
│ center: …………………………………………………………… [0.0] │
│ length: …………………………………………………………… 2.0 │
│ periodicity: ……………………………………………… (true,) │
│ current #cells: ……………………………………… 31 │
│ #leaf-cells: ……………………………………………… 16 │
│ maximum #cells: ……………………………………… 10000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation1D │
│ ═══════════════════════════════ │
│ #variables: ………………………………………………… 1 │
│ │ variable 1: …………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: ……………………………………… flux_central │
│ volume integral: …………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback │
│ ════════════════ │
│ interval: ……………………………………………………… 100 │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6) │
│ │ error 1: …………………………………………………… l2_error │
│ │ error 2: …………………………………………………… linf_error │
│ │ integral 1: …………………………………………… entropy_timederivative │
│ save analysis to file: …………………… no │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AliveCallback │
│ ═════════════ │
│ interval: ……………………………………………………… 10 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ………………………………………………… 0.0 │
│ Final time: ………………………………………………… 3.0 │
│ time integrator: …………………………………… RDPK3SpFSAL49 │
│ adaptive: ……………………………………………………… true │
│ abstol: …………………………………………………………… 1.0e-6 │
│ reltol: …………………………………………………………… 0.001 │
│ controller: ………………………………………………… PIDController(beta=[0.38, -0.18,…iter=default_dt_factor_limiter) │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: ……………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 5.31000000e-07 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 3.57776117e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 3523.832 MiB
#elements: 16
Variable: scalar
L2 error: 5.57368408e-06
Linf error: 1.21294882e-05
∑∂S/∂U ⋅ Uₜ : -5.00000000e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 9.9686e-02 │ sim. time: 3.4078e-01 (11.359%) │ run time: 2.3987e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 5.0123e-04 s
#timesteps: 30 │ Δt: 6.9464e-02 │ sim. time: 1.9556e+00 (65.188%) │ run time: 7.3164e-04 s
#timesteps: 40 │ Δt: 8.5031e-02 │ sim. time: 2.7005e+00 (90.016%) │ run time: 9.5995e-04 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 1.31167200e-03 s
Δt: 4.31148548e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.58799716e-08 s
PID: 4.05020983e-08 s
#DOFs per field: 64 alloc'd memory: 3523.873 MiB
#elements: 16
Variable: scalar
L2 error: 3.33421293e-05
Linf error: 1.31835569e-04
∑∂S/∂U ⋅ Uₜ : -1.23948350e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 2.05ms / 65.4% 98.1KiB / 48.1%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────
rhs! 417 919μs 68.4% 2.20μs 6.61KiB 14.0% 16.2B
source terms 417 386μs 28.8% 926ns 0.00B 0.0% 0.00B
~rhs!~ 417 276μs 20.5% 661ns 6.61KiB 14.0% 16.2B
volume integral 417 113μs 8.4% 271ns 0.00B 0.0% 0.00B
interface flux 417 44.5μs 3.3% 107ns 0.00B 0.0% 0.00B
prolong2interfaces 417 24.6μs 1.8% 58.9ns 0.00B 0.0% 0.00B
Jacobian 417 17.9μs 1.3% 42.9ns 0.00B 0.0% 0.00B
surface integral 417 17.6μs 1.3% 42.3ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 14.3μs 1.1% 34.2ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.0μs 1.0% 31.2ns 0.00B 0.0% 0.00B
boundary flux 417 12.2μs 0.9% 29.3ns 0.00B 0.0% 0.00B
analyze solution 2 424μs 31.6% 212μs 40.6KiB 86.0% 20.3KiB
────────────────────────────────────────────────────────────────────────────────Setting up a custom semidiscretization
Using a global constant is of course not really nice from a software engineering point of view. Thus, it can often be useful to collect additional data in the parameters of the ODEProblem. Thus, it is time to create our own semidiscretization. Here, we create a small wrapper of a standard semidiscretization of Trixi.jl and the current global time of the simulation.
struct CustomSemidiscretization{Semi, T} <: Trixi.AbstractSemidiscretization
semi::Semi
t::T
end
semi_custom = CustomSemidiscretization(semi, Ref(0.0))Main.CustomSemidiscretization{SemidiscretizationHyperbolic{TreeMesh{1, Trixi.SerialTree{1}}, LinearScalarAdvectionEquation1D{Float64}, typeof(Main.initial_condition), Trixi.BoundaryConditionPeriodic, typeof(Main.source_terms_standard), DGSEM{LobattoLegendreBasis{Float64, 4, SVector{4, Float64}, Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Matrix{Float64}, Matrix{Float64}}, SurfaceIntegralWeakForm{typeof(flux_central)}, VolumeIntegralWeakForm}, NamedTuple{(:elements, :interfaces, :boundaries), Tuple{Trixi.ElementContainer1D{Float64, Float64}, Trixi.InterfaceContainer1D{Float64}, Trixi.BoundaryContainer1D{Float64, Float64}}}}, Base.RefValue{Float64}}(SemidiscretizationHyperbolic(TreeMesh{1, Trixi.SerialTree{1}} with length 31, LinearScalarAdvectionEquation1D with one variable, initial_condition, boundary_condition_periodic, source_terms_standard, DG{Float64}(LobattoLegendreBasis{Float64}(polydeg=3), LobattoLegendreMortarL2{Float64}(polydeg=3), SurfaceIntegralWeakForm{typeof(flux_central)}(Trixi.flux_central), VolumeIntegralWeakForm()), cache(elements interfaces boundaries)), Base.RefValue{Float64}(0.0))To get pretty printing in the REPL, you can consider specializing
Base.show(io::IO, parameters::CustomSemidiscretization)Base.show(io::IO, ::MIME"text/plain", parameters::CustomSemidiscretization)
for your custom semidiscretiation.
Next, we create our own source terms that use the global time stored in the custom semidiscretiation.
source_terms_custom_semi = let semi_custom = semi_custom
function source_terms_custom_semi(u, x, t, equations)
t = semi_custom.t[]
return -initial_condition(x, t, equations)
end
endsource_terms_custom_semi (generic function with 1 method)We also create a custom ODE RHS to update the current global time stored in the custom semidiscretization. We unpack the standard semidiscretization created by Trixi.jl and pass it to Trixi.rhs!.
function rhs_semi_custom!(du_ode, u_ode, semi_custom, t)
semi_custom.t[] = t
Trixi.rhs!(du_ode, u_ode, semi_custom.semi, t)
endrhs_semi_custom! (generic function with 1 method)Finally, we set up an ODEProblem and solve it numerically.
ode_semi_custom = ODEProblem(rhs_semi_custom!,
ode.u0,
ode.tspan,
semi_custom)
sol_semi_custom = solve(ode_semi_custom, RDPK3SpFSAL49();
ode_default_options()...)retcode: Success
Interpolation: 1st order linear
t: 2-element Vector{Float64}:
0.0
3.0
u: 2-element Vector{Vector{Float64}}:
[-3.487868498008632e-16, -0.1083263689449018, -0.2803509724255764, -0.38268343236508945, -0.3826834323650901, -0.48051200262449845, -0.6263474196321541, -0.7071067811865472, -0.7071067811865478, -0.7795440397566659 … 0.7795440397566659, 0.7071067811865478, 0.7071067811865472, 0.6263474196321541, 0.48051200262449845, 0.3826834323650901, 0.38268343236508945, 0.2803509724255764, 0.1083263689449018, 3.487868498008632e-16]
[0.00039320882668261985, 0.005593037571101643, 0.014211415234637079, 0.019097988726944525, 0.01932495108774102, 0.02399181705537257, 0.03128895202056526, 0.0351157802915657, 0.03531125107864154, 0.03873867502331451 … -0.03840107811685022, -0.034871915240701996, -0.03475449432652092, -0.030829250470936792, -0.02350420904170509, -0.018777711432434757, -0.01859676365228886, -0.013658085806813441, -0.005028839697000977, 0.00017159954909493918]If we want to make use of additional functionality provided by Trixi.jl, e.g., for plotting, we need to implement a few additional specializations. In this case, we forward everything to the standard semidiscretization provided by Trixi.jl wrapped in our custom semidiscretization.
Base.ndims(semi::CustomSemidiscretization) = ndims(semi.semi)
function Trixi.mesh_equations_solver_cache(semi::CustomSemidiscretization)
Trixi.mesh_equations_solver_cache(semi.semi)
endNow, we can plot the numerical solution as usual.
plot(sol_semi_custom; label = "numerical sol.")
let
x = range(-1.0, 1.0; length = 200)
plot!(x, first.(initial_condition.(x, sol_semi_custom.t[end], equations)),
label = "analytical sol.", linestyle = :dash, legend = :topleft)
end
plot!(sol_semi_custom.u[1], semi, label = "u0", linestyle = :dot, legend = :topleft)
This also works with many callbacks as usual. However, the AnalysisCallback requires some special handling since it makes use of a performance counter contained in the standard semidiscretizations of Trixi.jl to report some performance metrics. Here, we forward all accesses to the performance counter to the wrapped semidiscretization.
function Base.getproperty(semi::CustomSemidiscretization, s::Symbol)
if s === :performance_counter
wrapped_semi = getfield(semi, :semi)
wrapped_semi.performance_counter
else
getfield(semi, s)
end
endMoreover, the AnalysisCallback also performs some error calculations. We also need to forward them to the wrapped semidiscretization.
function Trixi.calc_error_norms(func, u, t, analyzer,
semi::CustomSemidiscretization,
cache_analysis)
Trixi.calc_error_norms(func, u, t, analyzer,
semi.semi,
cache_analysis)
endNow, we can work with the callbacks used before as usual.
summary_callback = SummaryCallback()
analysis_interval = 100
analysis_callback = AnalysisCallback(semi_custom;
interval = analysis_interval)
alive_callback = AliveCallback(; analysis_interval)
callbacks = CallbackSet(summary_callback,
analysis_callback,
alive_callback)
sol = solve(ode_semi_custom, RDPK3SpFSAL49();
ode_default_options()..., callback = callbacks)
summary_callback()
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
Main.CustomSemidiscretization{SemidiscretizationHyperbolic{TreeMesh{1, Trixi.SerialTree{1}}, LinearScalarAdvectionEquation1D{Float64}, typeof(Main.initial_condition), Trixi.BoundaryConditionPeriodic, typeof(Main.source_terms_standard), DG{LobattoLegendreBasis{Float64, 4, StaticArraysCore.SArray{Tuple{4}, Float64, 1, 4}, Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}, Trixi.LobattoLegendreMortarL2{Float64, 4, Array{Float64, 2}, Array{Float64, 2}}, SurfaceIntegralWeakForm{typeof(flux_central)}, VolumeIntegralWeakForm}, NamedTuple{(:elements, :interfaces, :boundaries), Tuple{Trixi.ElementContainer1D{Float64, Float64}, Trixi.InterfaceContainer1D{Float64}, Trixi.BoundaryContainer1D{Float64, Float64}}}}, Base.RefValue{Float64}}(SemidiscretizationHyperbolic(TreeMesh{1, Trixi.SerialTree{1}} with length 31, LinearScalarAdvectionEquation1D with one variable, initial_condition, boundary_condition_periodic, source_terms_standard, DG{Float64}(LobattoLegendreBasis{Float64}(polydeg=3), LobattoLegendreMortarL2{Float64}(polydeg=3), SurfaceIntegralWeakForm{typeof(flux_central)}(Trixi.flux_central), VolumeIntegralWeakForm()), cache(elements interfaces boundaries)), Base.RefValue{Float64}(3.0))
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{1, Trixi.SerialTree{1}} │
│ ════════════════════════════════ │
│ center: …………………………………………………………… [0.0] │
│ length: …………………………………………………………… 2.0 │
│ periodicity: ……………………………………………… (true,) │
│ current #cells: ……………………………………… 31 │
│ #leaf-cells: ……………………………………………… 16 │
│ maximum #cells: ……………………………………… 10000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation1D │
│ ═══════════════════════════════ │
│ #variables: ………………………………………………… 1 │
│ │ variable 1: …………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: ……………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: …………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: ……………………………………… flux_central │
│ volume integral: …………………………………… VolumeIntegralWeakForm │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AnalysisCallback │
│ ════════════════ │
│ interval: ……………………………………………………… 100 │
│ analyzer: ……………………………………………………… LobattoLegendreAnalyzer{Float64}(polydeg=6) │
│ │ error 1: …………………………………………………… l2_error │
│ │ error 2: …………………………………………………… linf_error │
│ │ integral 1: …………………………………………… entropy_timederivative │
│ save analysis to file: …………………… no │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AliveCallback │
│ ═════════════ │
│ interval: ……………………………………………………… 10 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ………………………………………………… 0.0 │
│ Final time: ………………………………………………… 3.0 │
│ time integrator: …………………………………… RDPK3SpFSAL49 │
│ adaptive: ……………………………………………………… true │
│ abstol: …………………………………………………………… 1.0e-6 │
│ reltol: …………………………………………………………… 0.001 │
│ controller: ………………………………………………… PIDController(beta=[0.38, -0.18,…iter=default_dt_factor_limiter) │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: ……………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 6.41000000e-07 s
Δt: 0.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 3.56101871e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 3571.664 MiB
#elements: 16
Variable: scalar
L2 error: 5.57368408e-06
Linf error: 1.21294882e-05
∑∂S/∂U ⋅ Uₜ : -5.00000000e-01
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 10 │ Δt: 9.9686e-02 │ sim. time: 3.4078e-01 (11.359%) │ run time: 2.4325e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 5.0551e-04 s
#timesteps: 30 │ Δt: 6.9464e-02 │ sim. time: 1.9556e+00 (65.188%) │ run time: 7.3668e-04 s
#timesteps: 40 │ Δt: 8.5031e-02 │ sim. time: 2.7005e+00 (90.016%) │ run time: 9.6832e-04 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 1.33199100e-03 s
Δt: 4.31148548e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.59988038e-08 s
PID: 4.08516187e-08 s
#DOFs per field: 64 alloc'd memory: 3571.703 MiB
#elements: 16
Variable: scalar
L2 error: 3.33421293e-05
Linf error: 1.31835569e-04
∑∂S/∂U ⋅ Uₜ : -1.23948350e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 2.07ms / 65.2% 95.3KiB / 47.8%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────
rhs! 417 916μs 67.8% 2.20μs 6.61KiB 14.5% 16.2B
source terms 417 386μs 28.6% 926ns 0.00B 0.0% 0.00B
~rhs!~ 417 274μs 20.3% 656ns 6.61KiB 14.5% 16.2B
volume integral 417 112μs 8.3% 268ns 0.00B 0.0% 0.00B
interface flux 417 45.0μs 3.3% 108ns 0.00B 0.0% 0.00B
prolong2interfaces 417 24.7μs 1.8% 59.2ns 0.00B 0.0% 0.00B
Jacobian 417 18.0μs 1.3% 43.1ns 0.00B 0.0% 0.00B
surface integral 417 17.7μs 1.3% 42.3ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 14.2μs 1.1% 34.2ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.0μs 1.0% 31.3ns 0.00B 0.0% 0.00B
boundary flux 417 12.2μs 0.9% 29.2ns 0.00B 0.0% 0.00B
analyze solution 2 435μs 32.2% 217μs 38.9KiB 85.5% 19.4KiB
────────────────────────────────────────────────────────────────────────────────For even more advanced usage of custom semidiscretizations, you may look at the source code of the ones contained in Trixi.jl, e.g.,
SemidiscretizationHyperbolicParabolicSemidiscretizationEulerGravitySemidiscretizationEulerAcousticsSemidiscretizationCoupled
Package versions
These results were obtained using the following versions.
using InteractiveUtils
versioninfo()
using Pkg
Pkg.status(["Trixi", "OrdinaryDiffEq", "Plots"],
mode=PKGMODE_MANIFEST)Julia Version 1.9.4
Commit 8e5136fa297 (2023-11-14 08:46 UTC)
Build Info:
Official https://julialang.org/ release
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 4 × AMD EPYC 7763 64-Core Processor
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-14.0.6 (ORCJIT, znver3)
Threads: 1 on 4 virtual cores
Environment:
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
Status `~/work/Trixi.jl/Trixi.jl/docs/Manifest.toml`
⌅ [1dea7af3] OrdinaryDiffEq v6.66.0
[91a5bcdd] Plots v1.40.4
[a7f1ee26] Trixi v0.7.8 `~/work/Trixi.jl/Trixi.jl`
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated -m`This page was generated using Literate.jl.