Using the an implicit solver based on Ariadne with Trixi.jl
using Trixi
using Theseus
using CairoMakie
using LinearAlgebra
import Ariadne: JacobianOperatorNotes: Must disable both Polyester and LoopVectorization for Enzyme to be able to differentiate Trixi.jl
LocalPreferences.jl
[Trixi]
loop_vectorization = false
backend = "static"@assert Trixi._PREFERENCE_THREADING !== :polyester
@assert !Trixi._PREFERENCE_LOOPVECTORIZATIONLoad Trixi Example
trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl"), sol = nothing);
u = copy(ode.u0)
du = zero(ode.u0)
res = zero(ode.u0)
F! = Theseus.nonlinear_problem(Theseus.ImplicitEuler(), ode.f)
J = JacobianOperator(F!, res, u, (ode.u0, 1.0, du, ode.p, 0.0, (), 1))
out = zero(u)
v = zero(u)
# precompile
mul!(u, J, v)
F!(res, u, (ode.u0, 1.0, du, ode.p, 0.0, (), 1))
@time mul!(u, J, v)
@time F!(res, u, (ode.u0, 1.0, du, ode.p, 0.0, (), 1))[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.
0.000206 seconds (79 allocations: 7.844 KiB)
0.000142 seconds (7 allocations: 352 bytes)Cost of time(mul!) ≈ 2 * time(F!)
Solve using ODE interface
sol_trbdf2 = solve(
ode, Theseus.TRBDF2();
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks,
# verbose=1,
krylov_algo = :gmres,
# krylov_kwargs=(;verbose=1)
);
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ……………………………………… 2 │
│ mesh: ……………………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 341 │
│ equations: ………………………………………………………………… LinearScalarAdvectionEquation2D │
│ initial condition: …………………………………………… initial_condition_convergence_test │
│ boundary conditions: ……………………………………… Trixi.BoundaryConditionPeriodic │
│ source terms: ………………………………………………………… nothing │
│ solver: ………………………………………………………………………… DG │
│ total #DOFs per field: ………………………………… 4096 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2, Float64}} │
│ ═════════════════════════════════════════ │
│ center: ………………………………………………………………………… [0.0, 0.0] │
│ length: ………………………………………………………………………… 2.0 │
│ periodicity: …………………………………………………………… (true, true) │
│ current #cells: …………………………………………………… 341 │
│ #leaf-cells: …………………………………………………………… 256 │
│ current capacity: ……………………………………………… 341 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D │
│ ═══════════════════════════════ │
│ #variables: ……………………………………………………………… 1 │
│ │ variable 1: ………………………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: …………………………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: ………………………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ……………………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: …………………………………………………… FluxLaxFriedrichs(max_abs_speed) │
│ 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 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveSolutionCallback │
│ ════════════════════ │
│ interval: …………………………………………………………………… 100 │
│ solution variables: ………………………………………… cons2prim │
│ save initial solution: ………………………………… yes │
│ save final solution: ……………………………………… yes │
│ output directory: ……………………………………………… /home/runner/work/Ariadne.jl/…e.jl/docs/build/generated/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Hyperbolic: …………………………………………………… Returns{Float64}(1.6) │
│ CFL Parabolic: ……………………………………………………… Returns{Float64}(0.0) │
│ Interval: …………………………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ……………………………………………………………… 0.0 │
│ Final time: ……………………………………………………………… 1.0 │
│ time integrator: ………………………………………………… TRBDF2 │
│ adaptive: …………………………………………………………………… false │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: …………………………………………………………………… 1 │
│ threading backend: …………………………………………… static │
│ LoopVectorization: …………………………………………… disabled │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 1.06333980e+00 s
Δt: 1.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 2.83805542e-08 s
PID: Inf s
#DOFs per field: 4096 alloc'd memory: 572.010 MiB
#elements: 256 device memory: 0.000 MiB
Variable: scalar
L2 error: 3.94120825e-06
Linf error: 1.23670425e-05
∑∂S/∂U ⋅ Uₜ : -2.19124035e-16
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 18 run time: 2.95675672e+01 s
Δt: 5.55555556e-02 └── GC time: 3.71339549e-01 s (1.256%)
sim. time: 1.00000000e+00 (100.000%) time/DOF/rhs!: 2.10521483e-08 s
PID: 1.29613340e-05 s
#DOFs per field: 4096 alloc'd memory: 438.460 MiB
#elements: 256 device memory: 0.000 MiB
Variable: scalar
L2 error: 1.71224434e-04
Linf error: 2.52822142e-04
∑∂S/∂U ⋅ Uₜ : -1.97551223e-09
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 30.0s / 4.7% 1.61GiB / 2.2%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────────
I/O 3 1.16s 82.1% 388ms 16.9MiB 47.5% 5.62MiB
~I/O~ 3 968ms 68.3% 323ms 12.6MiB 35.5% 4.20MiB
save solution 2 188ms 13.3% 93.9ms 3.67MiB 10.3% 1.84MiB
get element variables 2 7.36ms 0.5% 3.68ms 622KiB 1.7% 311KiB
get node variables 2 931ns 0.0% 466ns 0.00B 0.0% 0.00B
save mesh 2 531ns 0.0% 266ns 0.00B 0.0% 0.00B
analyze solution 2 208ms 14.7% 104ms 18.6MiB 52.5% 9.31MiB
rhs! 533 45.8ms 3.2% 86.0μs 4.78KiB 0.0% 9.19B
volume integral 533 29.7ms 2.1% 55.8μs 0.00B 0.0% 0.00B
interface flux 533 4.76ms 0.3% 8.93μs 0.00B 0.0% 0.00B
surface integral 533 3.76ms 0.3% 7.06μs 0.00B 0.0% 0.00B
prolong2interfaces 533 3.53ms 0.2% 6.63μs 0.00B 0.0% 0.00B
Jacobian 533 2.76ms 0.2% 5.17μs 0.00B 0.0% 0.00B
~rhs!~ 533 594μs 0.0% 1.11μs 4.78KiB 0.0% 9.19B
reset ∂u/∂t 533 564μs 0.0% 1.06μs 0.00B 0.0% 0.00B
prolong2mortars 533 51.4μs 0.0% 96.4ns 0.00B 0.0% 0.00B
prolong2boundaries 533 39.6μs 0.0% 74.4ns 0.00B 0.0% 0.00B
mortar flux 533 37.5μs 0.0% 70.4ns 0.00B 0.0% 0.00B
boundary flux 533 17.1μs 0.0% 32.0ns 0.00B 0.0% 0.00B
source terms 533 16.6μs 0.0% 31.2ns 0.00B 0.0% 0.00B
calculate dt 19 15.5μs 0.0% 815ns 2.67KiB 0.0% 144B
────────────────────────────────────────────────────────────────────────────────────Plot the solution
We have to manually convert the sol since Theseus has it's own lightweight solution type.
plot(Trixi.PlotData2DTriangulated(sol_trbdf2.u[end], sol_trbdf2.prob.p))Solve using OrdinaryDiffEqSDIRK
import OrdinaryDiffEqSDIRK
import DifferentiationInterface: AutoFiniteDiff
sol_sdrik = solve(
ode, OrdinaryDiffEqSDIRK.TRBDF2(autodiff = AutoFiniteDiff());
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks,
adaptive = false
);
████████╗██████╗ ██╗██╗ ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
██║ ██████╔╝██║ ╚███╔╝ ██║
██║ ██╔══██╗██║ ██╔██╗ ██║
██║ ██║ ██║██║██╔╝ ██╗██║
╚═╝ ╚═╝ ╚═╝╚═╝╚═╝ ╚═╝╚═╝
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolic │
│ ════════════════════════════ │
│ #spatial dimensions: ……………………………………… 2 │
│ mesh: ……………………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 341 │
│ equations: ………………………………………………………………… LinearScalarAdvectionEquation2D │
│ initial condition: …………………………………………… initial_condition_convergence_test │
│ boundary conditions: ……………………………………… Trixi.BoundaryConditionPeriodic │
│ source terms: ………………………………………………………… nothing │
│ solver: ………………………………………………………………………… DG │
│ total #DOFs per field: ………………………………… 4096 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ TreeMesh{2, Trixi.SerialTree{2, Float64}} │
│ ═════════════════════════════════════════ │
│ center: ………………………………………………………………………… [0.0, 0.0] │
│ length: ………………………………………………………………………… 2.0 │
│ periodicity: …………………………………………………………… (true, true) │
│ current #cells: …………………………………………………… 341 │
│ #leaf-cells: …………………………………………………………… 256 │
│ current capacity: ……………………………………………… 341 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ LinearScalarAdvectionEquation2D │
│ ═══════════════════════════════ │
│ #variables: ……………………………………………………………… 1 │
│ │ variable 1: ………………………………………………………… scalar │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ DG{Float64} │
│ ═══════════ │
│ basis: …………………………………………………………………………… LobattoLegendreBasis{Float64}(polydeg=3) │
│ mortar: ………………………………………………………………………… LobattoLegendreMortarL2{Float64}(polydeg=3) │
│ surface integral: ……………………………………………… SurfaceIntegralWeakForm │
│ │ surface flux: …………………………………………………… FluxLaxFriedrichs(max_abs_speed) │
│ 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 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SaveSolutionCallback │
│ ════════════════════ │
│ interval: …………………………………………………………………… 100 │
│ solution variables: ………………………………………… cons2prim │
│ save initial solution: ………………………………… yes │
│ save final solution: ……………………………………… yes │
│ output directory: ……………………………………………… /home/runner/work/Ariadne.jl/…e.jl/docs/build/generated/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Hyperbolic: …………………………………………………… Returns{Float64}(1.6) │
│ CFL Parabolic: ……………………………………………………… Returns{Float64}(0.0) │
│ Interval: …………………………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration │
│ ════════════════ │
│ Start time: ……………………………………………………………… 0.0 │
│ Final time: ……………………………………………………………… 1.0 │
│ time integrator: ………………………………………………… TRBDF2 │
│ adaptive: …………………………………………………………………… false │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information │
│ ═══════════════════════ │
│ #threads: …………………………………………………………………… 1 │
│ threading backend: …………………………………………… static │
│ LoopVectorization: …………………………………………… disabled │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 0 run time: 6.71000000e-07 s
Δt: 1.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 2.21640625e-08 s
PID: Inf s
#DOFs per field: 4096 alloc'd memory: 763.582 MiB
#elements: 256 device memory: 0.000 MiB
Variable: scalar
L2 error: 3.94120825e-06
Linf error: 1.23670425e-05
∑∂S/∂U ⋅ Uₜ : -2.19124035e-16
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 18 run time: 3.89770212e+01 s
Δt: 5.55555556e-02 └── GC time: 1.04168967e+00 s (2.673%)
sim. time: 1.00000000e+00 (100.000%) time/DOF/rhs!: 2.21131036e-08 s
PID: 6.44979058e-08 s
#DOFs per field: 4096 alloc'd memory: 1025.216 MiB
#elements: 256 device memory: 0.000 MiB
Variable: scalar
L2 error: 1.71161806e-04
Linf error: 2.52798215e-04
∑∂S/∂U ⋅ Uₜ : -1.97752414e-09
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 39.0s / 34.3% 4.67GiB / 0.1%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────────
rhs! 148k 13.3s 99.9% 90.5μs 4.78KiB 0.1% 0.03B
volume integral 148k 8.59s 64.3% 58.2μs 0.00B 0.0% 0.00B
interface flux 148k 1.27s 9.5% 8.61μs 0.00B 0.0% 0.00B
surface integral 148k 1.24s 9.3% 8.41μs 0.00B 0.0% 0.00B
prolong2interfaces 148k 1.05s 7.8% 7.09μs 0.00B 0.0% 0.00B
Jacobian 148k 939ms 7.0% 6.36μs 0.00B 0.0% 0.00B
~rhs!~ 148k 140ms 1.0% 948ns 4.78KiB 0.1% 0.03B
reset ∂u/∂t 148k 98.2ms 0.7% 666ns 0.00B 0.0% 0.00B
prolong2boundaries 148k 8.04ms 0.1% 54.5ns 0.00B 0.0% 0.00B
prolong2mortars 148k 6.51ms 0.0% 44.1ns 0.00B 0.0% 0.00B
mortar flux 148k 6.39ms 0.0% 43.3ns 0.00B 0.0% 0.00B
boundary flux 148k 4.70ms 0.0% 31.9ns 0.00B 0.0% 0.00B
source terms 148k 4.60ms 0.0% 31.2ns 0.00B 0.0% 0.00B
analyze solution 2 10.2ms 0.1% 5.11ms 6.06MiB 96.7% 3.03MiB
I/O 3 2.29ms 0.0% 763μs 203KiB 3.2% 67.8KiB
save solution 2 2.28ms 0.0% 1.14ms 200KiB 3.1% 100KiB
~I/O~ 3 4.63μs 0.0% 1.54μs 1.70KiB 0.0% 581B
get element variables 2 1.91μs 0.0% 956ns 1.41KiB 0.0% 720B
get node variables 2 1.03μs 0.0% 516ns 0.00B 0.0% 0.00B
save mesh 2 201ns 0.0% 100ns 0.00B 0.0% 0.00B
calculate dt 19 75.7μs 0.0% 3.98μs 2.67KiB 0.0% 144B
────────────────────────────────────────────────────────────────────────────────────Plot the solution
plot(Trixi.PlotData2DTriangulated(sol_sdrik.u[end], sol_sdrik.prob.p))Increase CFL numbers
trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_basic.jl"), cfl = 10, sol = nothing);
sol = solve(
ode, Theseus.ImplicitEuler();
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
ode_default_options()..., callback = callbacks,
# verbose=1,
krylov_algo = :gmres,
# krylov_kwargs=(;verbose=1)
);
@show callbacks.discrete_callbacks[4]┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Hyperbolic: …………………………………………………… Returns{Int64}(10) │
│ CFL Parabolic: ……………………………………………………… Returns{Float64}(0.0) │
│ Interval: …………………………………………………………………… 1 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘Plot the solution
plot(Trixi.PlotData2DTriangulated(sol.u[end], sol.prob.p))This page was generated using Literate.jl.