Using an implicit-explicit (IMEX) Runge-Kutta solver based on Ariadne with Trixi.jl

using Trixi
using Theseus
using CairoMakie

Notes: You 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_LOOPVECTORIZATION

First call to load callbacks

trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion.jl"), sol = nothing);

stepsize_callback = StepsizeCallback(cfl = 1.0)

callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, stepsize_callback)

trixi_include(@__MODULE__, joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_diffusion.jl"), sol = nothing, callbacks = callbacks);
###############################################################################
[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.
[ Info: You just called `trixi_include`. Julia may now compile the code, please be patient.

run the simulation

sol = solve(
    ode,
    Theseus.SSP3433();
    dt = 0.01, # solve needs some value here but it will be overwritten by the stepsize_callback
    ode_default_options()..., callback = callbacks,
    krylov_algo = :gmres,
);

████████╗██████╗ ██╗██╗  ██╗██╗
╚══██╔══╝██╔══██╗██║╚██╗██╔╝██║
   ██║   ██████╔╝██║ ╚███╔╝ ██║
   ██║   ██╔══██╗██║ ██╔██╗ ██║
   ██║   ██║  ██║██║██╔╝ ██╗██║
   ╚═╝   ╚═╝  ╚═╝╚═╝╚═╝  ╚═╝╚═╝

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ SemidiscretizationHyperbolicParabolic                                                            │
│ ═════════════════════════════════════                                                            │
│ #spatial dimensions: ………………………… 2                                                                │
│ mesh: ………………………………………………………………… TreeMesh{2, Trixi.SerialTree{2, Float64}} with length 341        │
│ hyperbolic equations: ……………………… LinearScalarAdvectionEquation2D                                  │
│ parabolic equations: ………………………… LaplaceDiffusion2D                                               │
│ initial condition: ……………………………… initial_condition_diffusive_convergence_test                     │
│ source terms: …………………………………………… nothing                                                          │
│ source terms parabolic: ………………… nothing                                                          │
│ solver: …………………………………………………………… DG                                                               │
│ parabolic solver: ………………………………… ViscousFormulationBassiRebay1                                    │
│ 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                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ AliveCallback                                                                                    │
│ ═════════════                                                                                    │
│ interval: ……………………………………………………… 10                                                               │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ StepsizeCallback                                                                                 │
│ ════════════════                                                                                 │
│ CFL Advective: ………………………………………… Returns{Float64}(1.0)                                            │
│ CFL Diffusive: ………………………………………… Returns{Float64}(0.0)                                            │
│ Interval: ……………………………………………………… 1                                                                │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Time integration                                                                                 │
│ ════════════════                                                                                 │
│ Start time: ………………………………………………… 0.0                                                              │
│ Final time: ………………………………………………… 1.5                                                              │
│ time integrator: …………………………………… SSP3433                                                          │
│ adaptive: ……………………………………………………… false                                                            │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘
┌──────────────────────────────────────────────────────────────────────────────────────────────────┐
│ Environment information                                                                          │
│ ═══════════════════════                                                                          │
│ #threads: ……………………………………………………… 1                                                                │
│ threading backend: ……………………………… static                                                           │
│ LoopVectorization: ……………………………… disabled                                                         │
└──────────────────────────────────────────────────────────────────────────────────────────────────┘

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                  0                run time:       4.57000000e-07 s
 Δt:             1.00000000e-02                └── GC time:    0.00000000e+00 s (0.000%)
 sim. time:      0.00000000e+00 (0.000%)       time/DOF/rhs!:         NaN s
                                               PID:                   Inf s
 #DOFs per field:          4096                alloc'd memory:       1007.002 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       3.94120825e-06
 Linf error:     1.23670425e-05
 ∑∂S/∂U ⋅ Uₜ :  -1.23370054e-01
────────────────────────────────────────────────────────────────────────────────────────────────────

┌ 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:956
│   [7] preprocess_args
│     @ ./broadcast.jl:955
│   [8] preprocess
│     @ ./broadcast.jl:952
│   [9] preprocess_args
│     @ ./broadcast.jl:956
│  [10] preprocess_args
│     @ ./broadcast.jl:955
│  [11] preprocess
│     @ ./broadcast.jl:952
│  [12] preprocess_args
│     @ ./broadcast.jl:955
│  [13] preprocess
│     @ ./broadcast.jl:952
│  [14] override_bc_copyto!
│     @ ~/.julia/packages/Enzyme/C5rSn/src/compiler/interpreter.jl:818
│  [15] copyto!
│     @ ./broadcast.jl:925
│  [16] materialize!
│     @ ./broadcast.jl:883
│  [17] materialize!
│     @ ./broadcast.jl:880
│  [18] RKIMEX
│     @ ~/work/Ariadne.jl/Ariadne.jl/libs/Theseus/src/imex/imex.jl:33
└ @ Enzyme.Compiler ~/.julia/packages/Enzyme/C5rSn/src/rules/llvmrules.jl:812
┌ 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:956
│   [7] preprocess_args
│     @ ./broadcast.jl:955
│   [8] preprocess
│     @ ./broadcast.jl:952
│   [9] preprocess_args
│     @ ./broadcast.jl:956
│  [10] preprocess_args
│     @ ./broadcast.jl:955
│  [11] preprocess
│     @ ./broadcast.jl:952
│  [12] override_bc_copyto!
│     @ ~/.julia/packages/Enzyme/C5rSn/src/compiler/interpreter.jl:818
│  [13] copyto!
│     @ ./broadcast.jl:925
│  [14] materialize!
│     @ ./broadcast.jl:883
│  [15] materialize!
│     @ ./broadcast.jl:880
│  [16] RKIMEX
│     @ ~/work/Ariadne.jl/Ariadne.jl/libs/Theseus/src/imex/imex.jl:33
└ @ Enzyme.Compiler ~/.julia/packages/Enzyme/C5rSn/src/rules/llvmrules.jl:812
#timesteps:     10 │ Δt: 1.2500e-02 │ sim. time: 1.2500e-01 (8.333%)   │ run time: 4.9393e+01 s
#timesteps:     20 │ Δt: 1.2500e-02 │ sim. time: 2.5000e-01 (16.667%)  │ run time: 4.9544e+01 s
#timesteps:     30 │ Δt: 1.2500e-02 │ sim. time: 3.7500e-01 (25.000%)  │ run time: 4.9679e+01 s
#timesteps:     40 │ Δt: 1.2500e-02 │ sim. time: 5.0000e-01 (33.333%)  │ run time: 4.9824e+01 s
#timesteps:     50 │ Δt: 1.2500e-02 │ sim. time: 6.2500e-01 (41.667%)  │ run time: 4.9962e+01 s
#timesteps:     60 │ Δt: 1.2500e-02 │ sim. time: 7.5000e-01 (50.000%)  │ run time: 5.0099e+01 s
#timesteps:     70 │ Δt: 1.2500e-02 │ sim. time: 8.7500e-01 (58.333%)  │ run time: 5.0244e+01 s
#timesteps:     80 │ Δt: 1.2500e-02 │ sim. time: 1.0000e+00 (66.667%)  │ run time: 5.0396e+01 s
#timesteps:     90 │ Δt: 1.2500e-02 │ sim. time: 1.1250e+00 (75.000%)  │ run time: 5.0557e+01 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                100                run time:       5.07217941e+01 s
 Δt:             1.25000000e-02                └── GC time:    6.62086671e-01 s (1.305%)
 sim. time:      1.25000000e+00 (83.333%)      time/DOF/rhs!:  6.97139091e-08 s
                                               PID:            3.09575368e-05 s
 #DOFs per field:          4096                alloc'd memory:       1089.346 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       4.05704067e-05
 Linf error:     7.68749323e-05
 ∑∂S/∂U ⋅ Uₜ :  -1.04542969e-02
────────────────────────────────────────────────────────────────────────────────────────────────────

#timesteps:    110 │ Δt: 1.2500e-02 │ sim. time: 1.3750e+00 (91.667%)  │ run time: 5.0889e+01 s

────────────────────────────────────────────────────────────────────────────────────────────────────
 Simulation running 'LinearScalarAdvectionEquation2D' with DGSEM(polydeg=3)
────────────────────────────────────────────────────────────────────────────────────────────────────
 #timesteps:                120                run time:       5.10751563e+01 s
 Δt:             1.25000000e-02                └── GC time:    6.78699125e-01 s (1.329%)
 sim. time:      1.50000000e+00 (100.000%)     time/DOF/rhs!:  6.99734998e-08 s
                                               PID:            1.07525479e-06 s
 #DOFs per field:          4096                alloc'd memory:       1109.377 MiB
 #elements:                 256

 Variable:       scalar
 L2 error:       3.82554699e-05
 Linf error:     1.65324860e-04
 ∑∂S/∂U ⋅ Uₜ :  -6.38135743e-03
────────────────────────────────────────────────────────────────────────────────────────────────────

────────────────────────────────────────────────────────────────────────────────────────────────────
Trixi.jl simulation finished.  Final time: 1.5  Time steps: 120 (accepted), 120 (total)
────────────────────────────────────────────────────────────────────────────────────────────────────

───────────────────────────────────────────────────────────────────────────────────────
              Trixi.jl                        Time                    Allocations
                                     ───────────────────────   ────────────────────────
          Tot / % measured:               51.1s /   2.3%           4.24GiB /   0.2%

Section                      ncalls     time    %tot     avg     alloc    %tot      avg
───────────────────────────────────────────────────────────────────────────────────────
parabolic rhs!                5.48k    1.11s   96.3%   203μs   9.00KiB    0.1%    1.68B
  calculate gradient          5.48k    497ms   43.1%  90.7μs   3.33KiB    0.0%    0.62B
    volume integral           5.48k    300ms   26.0%  54.7μs     0.00B    0.0%    0.00B
    Jacobian                  5.48k   56.0ms    4.9%  10.2μs     0.00B    0.0%    0.00B
    interface flux            5.48k   40.9ms    3.5%  7.45μs     0.00B    0.0%    0.00B
    surface integral          5.48k   39.8ms    3.5%  7.27μs     0.00B    0.0%    0.00B
    prolong2interfaces        5.48k   35.8ms    3.1%  6.53μs     0.00B    0.0%    0.00B
    reset gradients           5.48k   16.8ms    1.5%  3.07μs     0.00B    0.0%    0.00B
    ~calculate gradient~      5.48k   6.15ms    0.5%  1.12μs   3.33KiB    0.0%    0.62B
    prolong2mortars           5.48k    698μs    0.1%   127ns     0.00B    0.0%    0.00B
    mortar flux               5.48k    473μs    0.0%  86.4ns     0.00B    0.0%    0.00B
    prolong2boundaries        5.48k    385μs    0.0%  70.3ns     0.00B    0.0%    0.00B
    boundary flux             5.48k    134μs    0.0%  24.5ns     0.00B    0.0%    0.00B
  volume integral             5.48k    313ms   27.1%  57.0μs     0.00B    0.0%    0.00B
  calculate viscous fluxes    5.48k    116ms   10.1%  21.2μs     0.00B    0.0%    0.00B
  interface flux              5.48k   36.8ms    3.2%  6.71μs     0.00B    0.0%    0.00B
  surface integral            5.48k   34.3ms    3.0%  6.25μs     0.00B    0.0%    0.00B
  transform variables         5.48k   33.6ms    2.9%  6.13μs     0.00B    0.0%    0.00B
  prolong2interfaces          5.48k   32.3ms    2.8%  5.90μs     0.00B    0.0%    0.00B
  Jacobian                    5.48k   27.4ms    2.4%  5.01μs     0.00B    0.0%    0.00B
  reset ∂u/∂t                 5.48k   9.14ms    0.8%  1.67μs     0.00B    0.0%    0.00B
  ~parabolic rhs!~            5.48k   8.96ms    0.8%  1.63μs   5.67KiB    0.1%    1.06B
  prolong2mortars             5.48k    618μs    0.1%   113ns     0.00B    0.0%    0.00B
  mortar flux                 5.48k    451μs    0.0%  82.2ns     0.00B    0.0%    0.00B
  prolong2boundaries          5.48k    434μs    0.0%  79.1ns     0.00B    0.0%    0.00B
  source terms parabolic      5.48k    121μs    0.0%  22.1ns     0.00B    0.0%    0.00B
  boundary flux               5.48k    118μs    0.0%  21.5ns     0.00B    0.0%    0.00B
rhs!                            480   39.7ms    3.4%  82.7μs   4.78KiB    0.1%    10.2B
  volume integral               480   24.7ms    2.1%  51.5μs     0.00B    0.0%    0.00B
  interface flux                480   3.80ms    0.3%  7.92μs     0.00B    0.0%    0.00B
  surface integral              480   3.46ms    0.3%  7.20μs     0.00B    0.0%    0.00B
  Jacobian                      480   3.26ms    0.3%  6.80μs     0.00B    0.0%    0.00B
  prolong2interfaces            480   3.06ms    0.3%  6.37μs     0.00B    0.0%    0.00B
  ~rhs!~                        480    723μs    0.1%  1.51μs   4.78KiB    0.1%    10.2B
  reset ∂u/∂t                   480    423μs    0.0%   882ns     0.00B    0.0%    0.00B
  prolong2boundaries            480   69.4μs    0.0%   145ns     0.00B    0.0%    0.00B
  prolong2mortars               480   65.1μs    0.0%   136ns     0.00B    0.0%    0.00B
  mortar flux                   480   51.7μs    0.0%   108ns     0.00B    0.0%    0.00B
  boundary flux                 480   14.6μs    0.0%  30.3ns     0.00B    0.0%    0.00B
  source terms                  480   11.4μs    0.0%  23.8ns     0.00B    0.0%    0.00B
analyze solution                  3   3.13ms    0.3%  1.04ms   7.24MiB   99.6%  2.41MiB
calculate dt                    121    246μs    0.0%  2.04μs   17.0KiB    0.2%     144B
───────────────────────────────────────────────────────────────────────────────────────

This page was generated using Literate.jl.