Edit on GitHub

Multirate Runge-Kutta Timestepping

In this tutorial, we shall explore the use of Multirate Runge-Kutta methods for the solution of nonautonomous (or non time-invariant) equations. For our model problem, we shall reuse the acoustic wave test in the GCM configuration. See its code for details on the model and parameters. For the purposes of this tutorial, we will only run the experiment for a total of 3600 simulation seconds. Details on this test case can be found in Sec. 4.3 of Francis X Giraldo , James F Kelly , Emil M Constantinescu (2013).

using ClimateMachine
const clima_dir = dirname(dirname(pathof(ClimateMachine)));
include(joinpath(
    clima_dir,
    "tutorials",
    "Numerics",
    "TimeStepping",
    "tutorial_acousticwave_config.jl",
))

FT = Float64;

The typical context for Multirate splitting methods is given by problems in which the tendency $\mathcal{T}$ is assumed to have single parts that operate on different time rates (such as a slow time scale and a fast time scale). A general form is given by

$\dot{\boldsymbol{q}} = \mathcal{T}(\boldsymbol{q}) \equiv {\mathcal{T}}_{f}(\boldsymbol{q}) + {\mathcal{T}}_{s}(\boldsymbol{q}),$

where the right-hand side has been split into a "fast" component ${\mathcal{T}_{f}}$, and a "slow" component ${\mathcal{T}_{s}}$.

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.

For a given time-step size $\Delta t$, the two-rate method in Martin Schlegel , Oswald Knoth , Martin Arnold , Ralf Wolke (2009) is summarized as the following:

\[\begin{align} \boldsymbol{Q}_1 &= \boldsymbol{q}(t_n), \\ \boldsymbol{r}_{i} &= \sum_{j=1}^{i-1}\tilde{a}^O_{ij} {\mathcal{T}_{s}}(\boldsymbol{Q}_{j}), \\ \boldsymbol{w}_{i,1} &= \boldsymbol{Q}_{i-1},\\ \boldsymbol{w}_{i,k} &= \boldsymbol{w}_{i,k-1} + \Delta t \tilde{c}_i^O \sum_{j=1}^{k-1}\tilde{a}^I_{k,j} \left(\frac{1}{\tilde{c}_i^O}\boldsymbol{r}_i + {\mathcal{T}_{f}}(\boldsymbol{w}_{i,j})\right),\\ & \quad\quad i = 2, \cdots, s^O + 1 \text{ and } k = 2, \cdots, s^I + 1,\nonumber \\ \boldsymbol{Q}_i &= \boldsymbol{w}_{i,s^I + 1} \end{align}\]

where the tilde parameters denote increments per RK stage:

\[\begin{align} \tilde{a}_{ij} &= \begin{cases} a_{i,j} - a_{i-1, j} & \text{if } i < s + 1 \\ b_j - a_{s,j} & \text{if } i = s + 1 \end{cases},\\ \tilde{c}_{i} &= \begin{cases} c_{i} - c_{i-1} & \text{if } i < s + 1 \\ 1 - c_{s} & \text{if } i = s + 1 \end{cases}, \end{align}\]

where the coefficients $a$, $b$, and $c$ correspond to the Butcher tableau for a given RK method. The superscripts $O$ and $I$ denote the outer (slow) and inner (fast) components of the multirate method respectively. Thus, tilde coefficients should be associated with the RK method indicated by the superscripts. In other words, the RK methods for the slow ${\mathcal{T}_{s}}$ and fast ${\mathcal{T}_{f}}$ components have Butcher tables given by:

\[\begin{align} \begin{array}{c|c} \boldsymbol{c}_{O} &\boldsymbol{A}_{O} \\ \hline & \boldsymbol{b}_O^T \end{array}, \quad \begin{array}{c|c} \boldsymbol{c}_{I} &\boldsymbol{A}_{I} \\ \hline & \boldsymbol{b}_I^T \end{array}, \end{align}\]

where $\boldsymbol{A}_O = \lbrace a_{i,j}^O\rbrace$, $\boldsymbol{b}_O = \lbrace b_i^O \rbrace$, and $c_O = \lbrace c_i^O \rbrace$ (similarly for $\boldsymbol{A}_I$, $\boldsymbol{b}_I$, and $\boldsymbol{c}_I$). The method described here is for an explicit RK outer method with $s$ stages. More details can be found in M Schlegel , O Knoth , M Arnold , R Wolke (2012).

The acoustic wave test case used in this tutorial represents a global-scale problem with inertia-gravity waves traveling around the entire planet. It has a hydrostatically balanced initial state that is given a pressure perturbation. This initial pressure perturbation causes an acoustic wave to travel to the antipode, coalesce, and return to the initial position. The exact solution of this test case is simple in that the (linear) acoustic theory allows one to verify the analytic speed of sound based on the thermodynamics variables. The initial condition is defined as a hydrostatically balanced atmosphere with background (reference) potential temperature.

ode_solver = ClimateMachine.MultirateSolverType(
    splitting_type = ClimateMachine.HEVISplitting(),
    slow_method = LSRK54CarpenterKennedy,
    fast_method = ARK2GiraldoKellyConstantinescu,
    implicit_solver_adjustable = true,
    timestep_ratio = 100,
)

timeend = FT(3600)
CFL = FT(5)
cfl_direction = HorizontalDirection()
run_acousticwave(ode_solver, CFL, cfl_direction, timeend);
ClimateMachine.array_type() = Array
┌ Info: Model composition
│     physics = ClimateMachine.Atmos.AtmosPhysics{Float64,Main.##423.EarthParameterSet,ClimateMachine.Atmos.HydrostaticState{Thermodynamics.TemperatureProfiles.DecayingTemperatureProfile{Float64},Float64},ClimateMachine.Atmos.TotalEnergyModel,ClimateMachine.Atmos.DryModel,ClimateMachine.Atmos.Compressible,ClimateMachine.TurbulenceClosures.ConstantDynamicViscosity{Float64,ClimateMachine.TurbulenceClosures.WithoutDivergence},ClimateMachine.TurbulenceConvection.NoTurbConv,ClimateMachine.TurbulenceClosures.NoHyperDiffusion,ClimateMachine.TurbulenceClosures.NoViscousSponge,ClimateMachine.Atmos.NoPrecipitation,ClimateMachine.Atmos.NoRadiation,ClimateMachine.Atmos.NoTracers,ClimateMachine.Atmos.NoLSForcing}(Main.##423.EarthParameterSet(), ClimateMachine.Atmos.HydrostaticState{Thermodynamics.TemperatureProfiles.DecayingTemperatureProfile{Float64},Float64}(Thermodynamics.TemperatureProfiles.DecayingTemperatureProfile{Float64}(300.0, 300.0, 8776.832091407436), 0.0, true), ClimateMachine.Atmos.TotalEnergyModel(), ClimateMachine.Atmos.DryModel(), ClimateMachine.Atmos.Compressible(), ClimateMachine.TurbulenceClosures.ConstantDynamicViscosity{Float64,ClimateMachine.TurbulenceClosures.WithoutDivergence}(0.0, ClimateMachine.TurbulenceClosures.WithoutDivergence()), ClimateMachine.TurbulenceConvection.NoTurbConv(), ClimateMachine.TurbulenceClosures.NoHyperDiffusion(), ClimateMachine.TurbulenceClosures.NoViscousSponge(), ClimateMachine.Atmos.NoPrecipitation(), ClimateMachine.Atmos.NoRadiation(), ClimateMachine.Atmos.NoTracers(), 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}},Main.##423.AcousticWaveSetup{Float64},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.##423.AcousticWaveSetup{Float64}(10000.0, 300.0, 3.0, 100.0, 1), ClimateMachine.Atmos.atmos_problem_init_state_auxiliary)
│     orientation = ClimateMachine.Orientations.SphericalOrientation()
│     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) │        () │
└──────────┴────────────────────────────┴─────────────────────────────────┴───────────┘

┌ Info: Establishing Atmos GCM configuration for GCM Driver: Acoustic wave test
│     precision               = Float64
│     horiz polynomial order  = 4
│     vert polynomial order   = 4
│     horiz cutoff order      = nothing
│     vert cutoff order       = nothing
│     # horiz elem            = 6
│     # vert elems            = 4
│     domain height           = 1.00e+04 m
│     MPI ranks               = 1
│     min(Δ_horz)             = 203711.48 m
│     min(Δ_vert)             = 431.68 m
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/driver_configs.jl:418
┌ Info: Initializing GCM Driver: Acoustic wave test
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/solver_configs.jl:185
┌ Info: Starting GCM Driver: Acoustic wave test
│     dt              = 1.80000e+03
│     timeend         =  3600.00
│     number of steps = 2
│     norm(Q)         = 9.5073559986777266e+13
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:802
┌ Info: Update
│     simtime =  1800.00 /  3600.00
│     wallclock = 00:01:21
│     efficiency (simtime / wallclock) =  22.1370
│     wallclock end (estimated) = 00:02:42
│     norm(Q) = 9.5073454032886547e+13
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Finished
│     norm(Q)            = 9.5073458796020750e+13
│     norm(Q) / norm(Q₀) = 9.9999893565827835e-01
│     norm(Q) - norm(Q₀) = -1.0119075651562500e+08
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:853

The interested reader can explore the combination of different slow and fast methods for Multirate solvers, consulting the ones available in ClimateMachine.jl, such as the Low-Storage-Runge-Kutta-methods, Strong-Stability-Preserving-RungeKutta-methods, and Additive-Runge-Kutta-methods.

References

This page was generated using Literate.jl.