21: 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, Float64}} │
│ ═════════════════════════════════════════ │
│ 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, Float64}} 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, Float64}} 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, Float64}} │
│ ═════════════════════════════════════════ │
│ 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: 8.12000000e-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.43852840e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 1376.083 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.3061e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 5.0355e-04 s
#timesteps: 30 │ Δt: 6.9722e-02 │ sim. time: 1.9558e+00 (65.192%) │ run time: 7.3560e-04 s
#timesteps: 40 │ Δt: 6.8292e-02 │ sim. time: 2.6511e+00 (88.370%) │ run time: 9.6540e-04 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 1.50910900e-03 s
Δt: 7.36971224e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.44852721e-08 s
PID: 4.09569095e-08 s
#DOFs per field: 64 alloc'd memory: 1380.619 MiB
#elements: 16
Variable: scalar
L2 error: 4.46838745e-05
Linf error: 2.31363584e-04
∑∂S/∂U ⋅ Uₜ : -1.23948223e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 2.50ms / 65.9% 9.67MiB / 77.3%
Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────
rhs! 417 881μs 53.4% 2.11μs 6.61KiB 0.1% 16.2B
source terms 417 361μs 21.9% 866ns 0.00B 0.0% 0.00B
~rhs!~ 417 277μs 16.8% 665ns 6.61KiB 0.1% 16.2B
volume integral 417 99.0μs 6.0% 237ns 0.00B 0.0% 0.00B
interface flux 417 39.9μs 2.4% 95.6ns 0.00B 0.0% 0.00B
prolong2interfaces 417 22.3μs 1.3% 53.4ns 0.00B 0.0% 0.00B
surface integral 417 22.1μs 1.3% 53.0ns 0.00B 0.0% 0.00B
Jacobian 417 18.6μs 1.1% 44.6ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 14.7μs 0.9% 35.1ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.6μs 0.8% 32.6ns 0.00B 0.0% 0.00B
boundary flux 417 12.6μs 0.8% 30.3ns 0.00B 0.0% 0.00B
analyze solution 2 770μs 46.6% 385μs 7.46MiB 99.9% 3.73MiB
─────────────────────────────────────────────────────────────────────────────────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, Float64}} 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, Float64}} │
│ ═════════════════════════════════════════ │
│ 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: 7.52000000e-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.43182932e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 1738.271 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.4232e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 6.1590e-04 s
#timesteps: 30 │ Δt: 6.9722e-02 │ sim. time: 1.9558e+00 (65.192%) │ run time: 9.2375e-04 s
#timesteps: 40 │ Δt: 6.8292e-02 │ sim. time: 2.6511e+00 (88.370%) │ run time: 1.2285e-03 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 2.75427300e-03 s
Δt: 7.36971224e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.55346516e-08 s
PID: 5.33589628e-08 s
#DOFs per field: 64 alloc'd memory: 1742.807 MiB
#elements: 16
Variable: scalar
L2 error: 4.46838745e-05
Linf error: 2.31363584e-04
∑∂S/∂U ⋅ Uₜ : -1.23948223e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 5.48ms / 70.4% 9.67MiB / 77.3%
Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────
analyze solution 2 2.95ms 76.5% 1.48ms 7.46MiB 99.9% 3.73MiB
rhs! 417 908μs 23.5% 2.18μs 6.61KiB 0.1% 16.2B
source terms 417 377μs 9.8% 905ns 0.00B 0.0% 0.00B
~rhs!~ 417 285μs 7.4% 683ns 6.61KiB 0.1% 16.2B
volume integral 417 100μs 2.6% 239ns 0.00B 0.0% 0.00B
interface flux 417 40.6μs 1.1% 97.3ns 0.00B 0.0% 0.00B
prolong2interfaces 417 24.1μs 0.6% 57.9ns 0.00B 0.0% 0.00B
surface integral 417 21.8μs 0.6% 52.4ns 0.00B 0.0% 0.00B
Jacobian 417 18.0μs 0.5% 43.1ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 15.1μs 0.4% 36.1ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.8μs 0.4% 33.0ns 0.00B 0.0% 0.00B
boundary flux 417 12.6μs 0.3% 30.3ns 0.00B 0.0% 0.00B
─────────────────────────────────────────────────────────────────────────────────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, Float64}, Float64}, 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::Trixi.ElementContainer1D{Float64, Float64}, interfaces::Trixi.InterfaceContainer1D{Float64}, boundaries::Trixi.BoundaryContainer1D{Float64, Float64}}}, Base.RefValue{Float64}}(SemidiscretizationHyperbolic(TreeMesh{1, Trixi.SerialTree{1, Float64}} 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.00039566016129114303, 0.005595530597144552, 0.01420823700056966, 0.019099085166061614, 0.01933615033096657, 0.02399145283812463, 0.03128767555719924, 0.03512318003345417, 0.035313794068257194, 0.0387373345584359 … -0.038403718625379946, -0.034873676580076514, -0.034746550097046656, -0.030828914911087947, -0.02350468228582584, -0.018776524261631816, -0.01859887840964364, -0.013657386582008402, -0.0050293813323077285, 0.00016360783023488103]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, Float64}, Float64}, 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::Trixi.ElementContainer1D{Float64, Float64}, interfaces::Trixi.InterfaceContainer1D{Float64}, boundaries::Trixi.BoundaryContainer1D{Float64, Float64}}}, Base.RefValue{Float64}}(SemidiscretizationHyperbolic(TreeMesh{1, Trixi.SerialTree{1, Float64}} 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, Float64}} │
│ ═════════════════════════════════════════ │
│ 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: 8.01000000e-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.44428297e-08 s
PID: Inf s
#DOFs per field: 64 alloc'd memory: 1179.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.4206e-04 s
#timesteps: 20 │ Δt: 8.5503e-02 │ sim. time: 1.2059e+00 (40.198%) │ run time: 5.5566e-04 s
#timesteps: 30 │ Δt: 6.9722e-02 │ sim. time: 1.9558e+00 (65.192%) │ run time: 8.0691e-04 s
#timesteps: 40 │ Δt: 6.8292e-02 │ sim. time: 2.6511e+00 (88.370%) │ run time: 1.0453e-03 s
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation1D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 45 run time: 5.46132800e-03 s
Δt: 7.36971224e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 3.00000000e+00 (100.000%) time/DOF/rhs!: 3.58297697e-08 s
PID: 4.41342551e-08 s
#DOFs per field: 64 alloc'd memory: 1184.468 MiB
#elements: 16
Variable: scalar
L2 error: 4.46838745e-05
Linf error: 2.31363584e-04
∑∂S/∂U ⋅ Uₜ : -1.23948223e-03
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished. Final time: 3.0 Time steps: 45 (accepted), 46 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────
─────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 6.49ms / 85.7% 9.91MiB / 77.6%
Section ncalls time %tot avg alloc %tot avg
─────────────────────────────────────────────────────────────────────────────────
analyze solution 2 4.67ms 84.0% 2.33ms 7.68MiB 99.9% 3.84MiB
rhs! 417 892μs 16.0% 2.14μs 6.61KiB 0.1% 16.2B
source terms 417 366μs 6.6% 877ns 0.00B 0.0% 0.00B
~rhs!~ 417 281μs 5.0% 673ns 6.61KiB 0.1% 16.2B
volume integral 417 100μs 1.8% 239ns 0.00B 0.0% 0.00B
interface flux 417 39.6μs 0.7% 95.0ns 0.00B 0.0% 0.00B
prolong2interfaces 417 24.4μs 0.4% 58.4ns 0.00B 0.0% 0.00B
surface integral 417 22.3μs 0.4% 53.4ns 0.00B 0.0% 0.00B
Jacobian 417 18.3μs 0.3% 43.9ns 0.00B 0.0% 0.00B
reset ∂u/∂t 417 14.8μs 0.3% 35.5ns 0.00B 0.0% 0.00B
prolong2boundaries 417 13.6μs 0.2% 32.6ns 0.00B 0.0% 0.00B
boundary flux 417 12.6μs 0.2% 30.2ns 0.00B 0.0% 0.00B
─────────────────────────────────────────────────────────────────────────────────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.10.7
Commit 4976d05258e (2024-11-26 15:57 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-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (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.9
[a7f1ee26] Trixi v0.9.12 `~/work/Trixi.jl/Trixi.jl`
Info Packages marked with ⌃ have new versions available and may be upgradable.This page was generated using Literate.jl.