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.000244 seconds (83 allocations: 8.906 KiB)
0.000141 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 │
│ maximum #cells: ……………………………………… 30000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ 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/Ariadne.jl/docs/build/generated/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Advective: ………………………………………… Returns{Float64}(1.6) │
│ CFL Diffusive: ………………………………………… 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: 9.78911368e-01 s
Δt: 1.00000000e+00 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 0.00000000e+00 (0.000%) time/DOF/rhs!: 2.80758057e-08 s
PID: Inf s
#DOFs per field: 4096 alloc'd memory: 721.600 MiB
#elements: 256
Variable: scalar
L2 error: 3.94120825e-06
Linf error: 1.23670425e-05
∑∂S/∂U ⋅ Uₜ : -2.19124035e-16
────────────────────────────────────────────────────────────────────────────────────────────────────
┌ Warning: TODO forward zero-set of memorycopy used memset rather than runtime type
│ Caused by:
│ Stacktrace:
│ [1] copy
│ @ ./array.jl:350
│ [2] unaliascopy
│ @ ./abstractarray.jl:1516
│ [3] unalias
│ @ ./abstractarray.jl:1500
│ [4] broadcast_unalias
│ @ ./broadcast.jl:946
│ [5] preprocess
│ @ ./broadcast.jl:953
│ [6] preprocess_args
│ @ ./broadcast.jl:955
│ [7] preprocess
│ @ ./broadcast.jl:952
│ [8] override_bc_copyto!
│ @ ~/.julia/packages/Enzyme/C5rSn/src/compiler/interpreter.jl:818
│ [9] copyto!
│ @ ./broadcast.jl:925
│ [10] materialize!
│ @ ./broadcast.jl:883
│ [11] materialize!
│ @ ./broadcast.jl:880
│ [12] TRBDF2
│ @ ~/work/Ariadne.jl/Ariadne.jl/libs/Theseus/src/Theseus.jl:139
└ @ Enzyme.Compiler ~/.julia/packages/Enzyme/C5rSn/src/rules/llvmrules.jl:812
────────────────────────────────────────────────────────────────────────────────────────────────────
Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
#timesteps: 18 run time: 3.36708565e+01 s
Δt: 5.55555556e-02 └── GC time: 1.44975603e+00 s (4.306%)
sim. time: 1.00000000e+00 (100.000%) time/DOF/rhs!: 2.25913360e-08 s
PID: 1.48843109e-05 s
#DOFs per field: 4096 alloc'd memory: 505.649 MiB
#elements: 256
Variable: scalar
L2 error: 1.71224434e-04
Linf error: 2.52822142e-04
∑∂S/∂U ⋅ Uₜ : -1.97551211e-09
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 34.0s / 4.1% 1.95GiB / 1.9%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────────
I/O 3 1.11s 79.6% 369ms 16.9MiB 45.5% 5.63MiB
~I/O~ 3 926ms 66.5% 309ms 12.6MiB 33.9% 4.20MiB
save solution 2 174ms 12.5% 87.2ms 3.67MiB 9.9% 1.84MiB
get element variables 2 6.88ms 0.5% 3.44ms 620KiB 1.6% 310KiB
get node variables 2 753ns 0.0% 376ns 0.00B 0.0% 0.00B
save mesh 2 453ns 0.0% 226ns 0.00B 0.0% 0.00B
analyze solution 2 235ms 16.9% 118ms 20.2MiB 54.5% 10.1MiB
rhs! 533 49.2ms 3.5% 92.3μs 4.78KiB 0.0% 9.19B
volume integral 533 30.9ms 2.2% 58.0μs 0.00B 0.0% 0.00B
interface flux 533 5.00ms 0.4% 9.38μs 0.00B 0.0% 0.00B
surface integral 533 4.58ms 0.3% 8.60μs 0.00B 0.0% 0.00B
prolong2interfaces 533 3.63ms 0.3% 6.81μs 0.00B 0.0% 0.00B
Jacobian 533 3.09ms 0.2% 5.81μs 0.00B 0.0% 0.00B
reset ∂u/∂t 533 1.09ms 0.1% 2.04μs 0.00B 0.0% 0.00B
~rhs!~ 533 694μs 0.0% 1.30μs 4.78KiB 0.0% 9.19B
prolong2mortars 533 66.1μs 0.0% 124ns 0.00B 0.0% 0.00B
prolong2boundaries 533 43.1μs 0.0% 80.9ns 0.00B 0.0% 0.00B
mortar flux 533 35.8μs 0.0% 67.1ns 0.00B 0.0% 0.00B
boundary flux 533 16.2μs 0.0% 30.5ns 0.00B 0.0% 0.00B
source terms 533 13.0μs 0.0% 24.4ns 0.00B 0.0% 0.00B
calculate dt 19 38.1μs 0.0% 2.01μs 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 │
│ maximum #cells: ……………………………………… 30000 │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ 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/Ariadne.jl/docs/build/generated/out │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback │
│ ════════════════ │
│ CFL Advective: ………………………………………… Returns{Float64}(1.6) │
│ CFL Diffusive: ………………………………………… 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: 2.84000000e-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.64396973e-08 s
PID: Inf s
#DOFs per field: 4096 alloc'd memory: 786.882 MiB
#elements: 256
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.78110389e+01 s
Δt: 5.55555556e-02 └── GC time: 0.00000000e+00 s (0.000%)
sim. time: 1.00000000e+00 (100.000%) time/DOF/rhs!: 1.95091534e-08 s
PID: 4.60205278e-08 s
#DOFs per field: 4096 alloc'd memory: 1028.921 MiB
#elements: 256
Variable: scalar
L2 error: 1.71161806e-04
Linf error: 2.52798214e-04
∑∂S/∂U ⋅ Uₜ : -1.97752408e-09
────────────────────────────────────────────────────────────────────────────────────────────────────
────────────────────────────────────────────────────────────────────────────────────
Trixi.jl Time Allocations
─────────────────────── ────────────────────────
Tot / % measured: 27.8s / 42.4% 245MiB / 2.1%
Section ncalls time %tot avg alloc %tot avg
────────────────────────────────────────────────────────────────────────────────────
rhs! 148k 11.8s 100.0% 79.8μs 4.78KiB 0.1% 0.03B
volume integral 148k 7.44s 63.2% 50.4μs 0.00B 0.0% 0.00B
interface flux 148k 1.21s 10.3% 8.23μs 0.00B 0.0% 0.00B
surface integral 148k 1.01s 8.6% 6.84μs 0.00B 0.0% 0.00B
Jacobian 148k 941ms 8.0% 6.38μs 0.00B 0.0% 0.00B
prolong2interfaces 148k 918ms 7.8% 6.22μs 0.00B 0.0% 0.00B
~rhs!~ 148k 118ms 1.0% 802ns 4.78KiB 0.1% 0.03B
reset ∂u/∂t 148k 114ms 1.0% 774ns 0.00B 0.0% 0.00B
prolong2boundaries 148k 5.82ms 0.0% 39.4ns 0.00B 0.0% 0.00B
prolong2mortars 148k 5.02ms 0.0% 34.0ns 0.00B 0.0% 0.00B
mortar flux 148k 4.28ms 0.0% 29.0ns 0.00B 0.0% 0.00B
boundary flux 148k 3.38ms 0.0% 22.9ns 0.00B 0.0% 0.00B
source terms 148k 2.86ms 0.0% 19.4ns 0.00B 0.0% 0.00B
I/O 3 1.89ms 0.0% 629μs 203KiB 3.9% 67.8KiB
save solution 2 1.88ms 0.0% 938μs 200KiB 3.9% 100KiB
~I/O~ 3 6.24μs 0.0% 2.08μs 1.70KiB 0.0% 581B
get element variables 2 1.38μs 0.0% 688ns 1.41KiB 0.0% 720B
get node variables 2 939ns 0.0% 470ns 0.00B 0.0% 0.00B
save mesh 2 160ns 0.0% 80.0ns 0.00B 0.0% 0.00B
analyze solution 2 1.58ms 0.0% 792μs 4.87MiB 95.9% 2.44MiB
calculate dt 19 81.1μs 0.0% 4.27μs 2.67KiB 0.1% 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 Advective: ………………………………………… Returns{Int64}(10) │
│ CFL Diffusive: ………………………………………… 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.