Multirate Infinitesimal Step (MIS) Timestepping
In this tutorial, we shall explore the use of Multirate Infinitesimal Step (MIS) methods for the solution of nonautonomous (or non time-invariant) equations. For our model problem, we shall reuse the rising thermal bubble tutorial. See its tutorial page for details on the model and parameters. For the purposes of this tutorial, we will only run the experiment for a total of 500 simulation seconds.
using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
clima_dir,
"tutorials",
"Numerics",
"TimeStepping",
"tutorial_risingbubble_config.jl",
))
FT = Float64;Referencing the formulation introduced in the previous Multirate RK methods tutorial, we can describe Multirate Infinitesimal Step (MIS) methods by
\[\begin{align} v_i (0) &= q^n + \sum_{j=1}^{i-1} \alpha_{ij} (Q^{(j)} - q^n) \\ \frac{dv_i}{d\tau} &= \sum_{j=1}^{i-1} \frac{\gamma_{ij}}{d_i \Delta t} (Q^{(j)} - q^n) + \sum_{j=1}^i \frac{\beta_{ij}}{d_i} \mathcal{T}_S (Q^{(j)}, t + \Delta t c_i) + \mathcal{T}_F(v_i, t^n + \Delta t \tilde c_i + \frac{c_i - \tilde c_i}{d_i} \tau), \quad \tau \in [0, \Delta t d_i] \\ Q^{(i)} &= v_i(\Delta t d_i), \end{align}\]
where we have used the the stage values $Q^{(i)} = v_i(\tau_i)$ as the solution to the inner ODE problem, ${\mathcal{T}_{s}}$ for the slow component, and `{\mathcal{T}_{f}} for the fast one, as in the Multirate RK methods tutorial.
Referencing the canonical form introduced in Time integration, both ${\mathcal{T}_{f}}$ and ${\mathcal{T}_{s}}$ could be discretized either explicitly or implicitly, hence, they could belong to either $\mathcal{F}(t, \boldsymbol{q})$ or $\mathcal{G}(t, \boldsymbol{q})$ term.
The method is defined in terms of the lower-triangular matrices $\alpha$, $\beta$ and $\gamma$, with $d_i = \sum_j \beta_{ij}$, $c_i = (I - \alpha - \gamma)^{-1} d$ and $\tilde c = \alpha c$. More details can be found in J. Wensch , Knoth O. , A. Galant (2009) and Oswald Knoth , Joerg Wensch (2014).
ode_solver = ClimateMachine.MISSolverType(;
mis_method = MIS2,
fast_method = LSRK144NiegemannDiehlBusch,
nsubsteps = (40,),
)
timeend = FT(500)
CFL = FT(20)
run_simulation(ode_solver, CFL, timeend);ClimateMachine.array_type() = Array
┌ Info: Model composition
│ physics = ClimateMachine.Atmos.AtmosPhysics{Float64,Main.##418.EarthParameterSet,ClimateMachine.Atmos.HydrostaticState{Thermodynamics.TemperatureProfiles.DryAdiabaticProfile{Float64},Float64},ClimateMachine.Atmos.TotalEnergyModel,ClimateMachine.Atmos.DryModel,ClimateMachine.Atmos.Compressible,ClimateMachine.TurbulenceClosures.SmagorinskyLilly{Float64},ClimateMachine.TurbulenceConvection.NoTurbConv,ClimateMachine.TurbulenceClosures.NoHyperDiffusion,ClimateMachine.TurbulenceClosures.NoViscousSponge,ClimateMachine.Atmos.NoPrecipitation,ClimateMachine.Atmos.NoRadiation,ClimateMachine.Atmos.NTracers{4,Float64},ClimateMachine.Atmos.NoLSForcing}(Main.##418.EarthParameterSet(), ClimateMachine.Atmos.HydrostaticState{Thermodynamics.TemperatureProfiles.DryAdiabaticProfile{Float64},Float64}(Thermodynamics.TemperatureProfiles.DryAdiabaticProfile{Float64}(300.0, 0.0), 0.0, true), ClimateMachine.Atmos.TotalEnergyModel(), ClimateMachine.Atmos.DryModel(), ClimateMachine.Atmos.Compressible(), ClimateMachine.TurbulenceClosures.SmagorinskyLilly{Float64}(0.21), ClimateMachine.TurbulenceConvection.NoTurbConv(), ClimateMachine.TurbulenceClosures.NoHyperDiffusion(), ClimateMachine.TurbulenceClosures.NoViscousSponge(), ClimateMachine.Atmos.NoPrecipitation(), ClimateMachine.Atmos.NoRadiation(), ClimateMachine.Atmos.NTracers{4,Float64}([1.0, 2.0, 3.0, 4.0]), ClimateMachine.Atmos.NoLSForcing())
│ problem = ClimateMachine.Atmos.AtmosProblem{Tuple{ClimateMachine.Atmos.AtmosBC{ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip},ClimateMachine.Atmos.Insulating,ClimateMachine.Atmos.Impermeable,ClimateMachine.Atmos.OutflowPrecipitation,ClimateMachine.Atmos.ImpermeableTracer,ClimateMachine.TurbulenceConvection.NoTurbConvBC},ClimateMachine.Atmos.AtmosBC{ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip},ClimateMachine.Atmos.Insulating,ClimateMachine.Atmos.Impermeable,ClimateMachine.Atmos.OutflowPrecipitation,ClimateMachine.Atmos.ImpermeableTracer,ClimateMachine.TurbulenceConvection.NoTurbConvBC}},typeof(Main.##418.init_risingbubble!),typeof(ClimateMachine.Atmos.atmos_problem_init_state_auxiliary)}((ClimateMachine.Atmos.AtmosBC{ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip},ClimateMachine.Atmos.Insulating,ClimateMachine.Atmos.Impermeable,ClimateMachine.Atmos.OutflowPrecipitation,ClimateMachine.Atmos.ImpermeableTracer,ClimateMachine.TurbulenceConvection.NoTurbConvBC}(ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip}(ClimateMachine.Atmos.FreeSlip()), ClimateMachine.Atmos.Insulating(), ClimateMachine.Atmos.Impermeable(), ClimateMachine.Atmos.OutflowPrecipitation(), ClimateMachine.Atmos.ImpermeableTracer(), ClimateMachine.TurbulenceConvection.NoTurbConvBC()), ClimateMachine.Atmos.AtmosBC{ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip},ClimateMachine.Atmos.Insulating,ClimateMachine.Atmos.Impermeable,ClimateMachine.Atmos.OutflowPrecipitation,ClimateMachine.Atmos.ImpermeableTracer,ClimateMachine.TurbulenceConvection.NoTurbConvBC}(ClimateMachine.Atmos.Impenetrable{ClimateMachine.Atmos.FreeSlip}(ClimateMachine.Atmos.FreeSlip()), ClimateMachine.Atmos.Insulating(), ClimateMachine.Atmos.Impermeable(), ClimateMachine.Atmos.OutflowPrecipitation(), ClimateMachine.Atmos.ImpermeableTracer(), ClimateMachine.TurbulenceConvection.NoTurbConvBC())), Main.##418.init_risingbubble!, ClimateMachine.Atmos.atmos_problem_init_state_auxiliary)
│ orientation = ClimateMachine.Orientations.FlatOrientation()
│ source = DispatchedTuples.DispatchedTuple{Tuple{Tuple{ClimateMachine.Atmos.Momentum,ClimateMachine.Atmos.Gravity}},DispatchedTuples.NoDefaults} with 1 entries:
│ ClimateMachine.Atmos.Momentum() => ClimateMachine.Atmos.Gravity()
│ default => ()
│
│ data_config = nothing
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/driver_configs.jl:188
PDE: ∂_t Y_i + (∇•F_1(Y))_i + (∇•F_2(Y,G)))_i = (S(Y,G))_i
┌────────────┬────────────────────────────┬─────────────────────────────────┬───────────┐
│ Equation │ Flux{FirstOrder} │ Flux{SecondOrder} │ Source │
│ (Y_i) │ (F_1) │ (F_2) │ (S) │
├────────────┼────────────────────────────┼─────────────────────────────────┼───────────┤
│ Mass │ (Advect) │ () │ () │
│ Momentum │ (Advect, PressureGradient) │ (ViscousStress) │ (Gravity) │
│ Energy │ (Advect, Pressure) │ (ViscousFlux, DiffEnthalpyFlux) │ () │
│ Tracers{4} │ (Advect) │ (Diffusion) │ () │
└────────────┴────────────────────────────┴─────────────────────────────────┴───────────┘
┌ Info: Establishing Atmos LES configuration for DryRisingBubble
│ precision = Float64
│ horiz polynomial order = 4
│ vert polynomial order = 4
│ horiz cutoff order = nothing
│ vert cutoff order = nothing
│ domain_min = 0.00 m, 0.00 m, 0.00 m
│ domain_max = 10000.00 m, 500.00 m, 10000.00 m
│ resolution = 125x125x125
│ MPI ranks = 1
│ min(Δ_horz) = 86.34 m
│ min(Δ_vert) = 86.34 m
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/driver_configs.jl:291
┌ Info: Initializing DryRisingBubble
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/solver_configs.jl:185
solver_config.ode_solver_type = MISSolverType{ClimateMachine.SlowFastSplitting}(ClimateMachine.SlowFastSplitting(), ClimateMachine.Atmos.AtmosAcousticGravityLinearModel, ClimateMachine.ODESolvers.MIS2, ClimateMachine.ODESolvers.LSRK144NiegemannDiehlBusch, (40,), false)
┌ Info: Starting DryRisingBubble
│ dt = 4.95050e+00
│ timeend = 500.00
│ number of steps = 101
│ norm(Q) = 5.3834517091450920e+09
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:802
┌ Info: Update
│ simtime = 4.95 / 500.00
│ wallclock = 00:01:38
│ efficiency (simtime / wallclock) = 0.0502
│ wallclock end (estimated) = 02:46:07
│ norm(Q) = 5.3834240236291618e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 49.50 / 500.00
│ wallclock = 00:02:42
│ efficiency (simtime / wallclock) = 0.3042
│ wallclock end (estimated) = 00:27:23
│ norm(Q) = 5.3834650176991339e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 99.01 / 500.00
│ wallclock = 00:03:46
│ efficiency (simtime / wallclock) = 0.4364
│ wallclock end (estimated) = 00:19:05
│ norm(Q) = 5.3833240611431246e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 148.51 / 500.00
│ wallclock = 00:04:50
│ efficiency (simtime / wallclock) = 0.5116
│ wallclock end (estimated) = 00:16:17
│ norm(Q) = 5.3834020452896700e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 198.02 / 500.00
│ wallclock = 00:05:52
│ efficiency (simtime / wallclock) = 0.5614
│ wallclock end (estimated) = 00:14:50
│ norm(Q) = 5.3834917108780994e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 247.52 / 500.00
│ wallclock = 00:06:54
│ efficiency (simtime / wallclock) = 0.5974
│ wallclock end (estimated) = 00:13:56
│ norm(Q) = 5.3835253391216793e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 301.98 / 500.00
│ wallclock = 00:07:59
│ efficiency (simtime / wallclock) = 0.6294
│ wallclock end (estimated) = 00:13:14
│ norm(Q) = 5.3835701374433317e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 356.44 / 500.00
│ wallclock = 00:09:05
│ efficiency (simtime / wallclock) = 0.6538
│ wallclock end (estimated) = 00:12:44
│ norm(Q) = 5.3834984434495621e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 410.89 / 500.00
│ wallclock = 00:10:10
│ efficiency (simtime / wallclock) = 0.6728
│ wallclock end (estimated) = 00:12:23
│ norm(Q) = 5.3834351201074266e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│ simtime = 465.35 / 500.00
│ wallclock = 00:11:15
│ efficiency (simtime / wallclock) = 0.6885
│ wallclock end (estimated) = 00:12:06
│ norm(Q) = 5.3834971668947382e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Finished
│ norm(Q) = 5.3835509731386242e+09
│ norm(Q) / norm(Q₀) = 1.0000184387264706e+00
│ norm(Q) - norm(Q₀) = 9.9263993532180786e+04
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:853
┌ Info: Euclidean distance
│ norm(Q - Qe) = 4.1623723276456207e+07
│ norm(Q - Qe) / norm(Qe) = 7.7317909633606010e-03
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:869
The reader can compare the Courant number (denoted by CFL in the code snippet) used in this example, with the adopted in the single-rate explicit timestepping tutorial page in which we use the same scheme as the fast method employed in this case, and notice that with this MIS method we are able to take a much larger Courant number.
References
This page was generated using Literate.jl.