Single-rate Explicit Timestepping

In this tutorial, we shall explore the use of explicit Runge-Kutta 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 100 simulation seconds.

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

FT = Float64;

After discretizing the spatial terms in the equation, the semi-discretization of the governing equations have the form:

$\begin{aligned} \frac{\mathrm{d} \boldsymbol{q}}{ \mathrm{d} t} &= M^{-1}\left(M S + D^{T} M (F^{adv} + F^{visc}) + \sum_{f=1}^{N_f} L^T M_f(\widehat{F}^{adv} + \widehat{F}^{visc}) \right) \equiv \mathcal{T}(\boldsymbol{q}). \end{aligned}$

Referencing the canonical form introduced in Time integration we have that in any explicit formulation $\mathcal{F}(t, \boldsymbol{q}) \equiv 0$ and, in this particular forumlation, $\mathcal{T}(t, \boldsymbol{q}) \equiv \mathcal{G}(t, \boldsymbol{q})$.

The time step restriction for an explicit method must satisfy the stable Courant number for the specific time-integrator and must be selected from the following constraints

$\Delta t_{\mathrm{explicit}} = min \left( \frac{C \Delta x_i}{u_i + a}, \frac{C \Delta x_i^2}{\nu} \right)$

where $C$ is the stable Courant number, $u_i$ denotes the velocity components, $a$ the speed of sound, $\Delta x_i$ the grid spacing (non-uniform in case of spectral element methods) along the direction $(x_1,x_2,x_3)$, and $\nu$ the kinematic viscosity. The first term on the right is the time step condition due to the non-dissipative components, while the second term to the dissipation. For explicit time-integrators, we have to find the minimum time step that satisfies this condition along all three spatial directions.

Runge-Kutta methods

A single step of an $s$-stage Runge-Kutta (RK) method for solving the resulting ODE problem presented above and can be expressed as the following:

\[\begin{align} \boldsymbol{q}^{n+1} = \boldsymbol{q}^n + \Delta t \sum_{i=1}^{s} b_i \mathcal{T}(\boldsymbol{Q}^i), \end{align}\]

where $\boldsymbol{\mathcal{T}}(\boldsymbol{Q}^i)$ is the evaluation of the right-hand side tendency at the stage value $\boldsymbol{Q}^i$, defined at each stage of the RK method:

\[\begin{align} \boldsymbol{Q}^i = \boldsymbol{q}^{n} + \Delta t \sum_{j=1}^{s} a_{i,j} \mathcal{T}(\boldsymbol{Q}^j). \end{align}\]

The first stage is initialized using the field at the previous time step: $\boldsymbol{Q}^{1} \leftarrow \boldsymbol{q}^n$.

In the above expressions, we define $\boldsymbol{A} = \lbrace a_{i,j} \rbrace \in \mathbb{R}^{s\times s}$, $\boldsymbol{b} = \lbrace b_i \rbrace \in \mathbb{R}^s$, and $\boldsymbol{c} = \lbrace c_i \rbrace \in \mathbb{R}^s$ as the characteristic coefficients of a given RK method. This means we can associate any RK method with its so-called Butcher tableau:

\[\begin{align} \begin{array}{c|c} \boldsymbol{c} &\boldsymbol{A}\\ \hline & \boldsymbol{b}^T \end{array} = \begin{array}{c|c c c c} c_1 & a_{1,1} & a_{1,2} & \cdots & a_{1,s}\\ c_2 & a_{2,1} & a_{2,2} & \cdots & a_{2,s}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ c_s & a_{s,1} & a_{s,2} & \cdots & a_{s,s}\\ \hline & b_1 & b_2 & \cdots & b_s \end{array}. \end{align}\]

The vector $\boldsymbol{c}$ is often called the consistency vector, and is typically subjected to the row-sum condition:

$c_i = \sum_{j=1}^{s} a_{i,j}, \quad \forall i = 1, \cdots, s.$

This simplifies the order conditions for higher-order RK methods. For more information on general RK methods, we refer the interested reader to Ch. 5.2 of Kendall Atkinson , Weimin Han , David E Stewart (2011).

Low-storage Runge-Kutta (LSRK) methods

ClimateMachine.jl contains the following low-storage methods:

Forward Euler LowStorageRungeKutta2N,
A 5-stage 4th-order Runge-Kutta method of Carpenter and Kennedy LSRK54CarpenterKennedy
A 14-stage 4th-order Runge-Kutta method developed by Niegemann, Diehl, and Busch LSRK144NiegemannDiehlBusch.

To start, let's try using the 5-stage method: LSRK54CarpenterKennedy.

As is the case for all explicit methods, we are limited by the fastest propogating waves described by our governing equations. In our case, these are the acoustic waves (with approximate wave speed given by the speed of sound of 343 m/s). For the rising bubble example used here, we use 4th order polynomials in a discontinuous Galerkin approximation, with a domain resolution of 125 meters in each spatial direction. This gives an effective minimanl nodal distance (distance between LGL nodes) of 86 meters over the entire mesh. Using the equation for the explcit time-step above, we can determine the $\Delta t$ by specifying the desired Courant number $C$ (denoted CFL in the code below). In our case, a heuristically determined value of 0.4 is used.

timeend = FT(100)
ode_solver =
    ClimateMachine.ExplicitSolverType(solver_method = LSRK54CarpenterKennedy)
CFL = FT(0.4)
run_simulation(ode_solver, CFL, timeend);

ClimateMachine.array_type() = Array
┌ Info: Model composition
│     physics = ClimateMachine.Atmos.AtmosPhysics{Float64,Main.##393.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.##393.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.##393.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.##393.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 = ExplicitSolverType(ClimateMachine.ODESolvers.LSRK54CarpenterKennedy)
┌ Info: Starting DryRisingBubble
│     dt              = 9.94036e-02
│     timeend         =   100.00
│     number of steps = 1006
│     norm(Q)         = 5.3834517091450920e+09
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:802
┌ Info: Update
│     simtime =     0.80 /   100.00
│     wallclock = 00:01:00
│     efficiency (simtime / wallclock) =   0.0132
│     wallclock end (estimated) = 02:06:33
│     norm(Q) = 5.3834501606460686e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    18.59 /   100.00
│     wallclock = 00:02:00
│     efficiency (simtime / wallclock) =   0.1547
│     wallclock end (estimated) = 00:10:46
│     norm(Q) = 5.3834165047290468e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    35.69 /   100.00
│     wallclock = 00:03:00
│     efficiency (simtime / wallclock) =   0.1977
│     wallclock end (estimated) = 00:08:25
│     norm(Q) = 5.3833080082318859e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    54.37 /   100.00
│     wallclock = 00:04:00
│     efficiency (simtime / wallclock) =   0.2258
│     wallclock end (estimated) = 00:07:22
│     norm(Q) = 5.3834528277298279e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    73.66 /   100.00
│     wallclock = 00:05:00
│     efficiency (simtime / wallclock) =   0.2448
│     wallclock end (estimated) = 00:06:48
│     norm(Q) = 5.3834548619051781e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    94.93 /   100.00
│     wallclock = 00:06:00
│     efficiency (simtime / wallclock) =   0.2630
│     wallclock end (estimated) = 00:06:20
│     norm(Q) = 5.3833067714414816e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Finished
│     norm(Q)            = 5.3833346937422333e+09
│     norm(Q) / norm(Q₀) = 9.9997826387062039e-01
│     norm(Q) - norm(Q₀) = -1.1701540285873413e+05
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:853
┌ Info: Euclidean distance
│     norm(Q - Qe)            = 7.6157180010283887e+06
│     norm(Q - Qe) / norm(Qe) = 1.4146533511373853e-03
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:869

What if we wish to take a larger timestep size? We could try to increase the target Courant number, say $C = 1.7$, and re-run the simulation. But this would break! In fact, in this case the numerical scheme would fall outside of the stability region and the simulation would crash. This occurs when the time step exceeds the maximal stable time-step size of the method. For the 5-stage method, one can typically get away with using time-step sizes corresponding to a Courant number of $C \approx 0.4$ but typically not much larger. In contrast, we can use an LSRK method with a larger stability region. Let's illustrate this by using the 14-stage method with a $C = 1.7$ instead.

ode_solver = ClimateMachine.ExplicitSolverType(
    solver_method = LSRK144NiegemannDiehlBusch,
)
CFL = FT(1.7)
run_simulation(ode_solver, CFL, timeend);

ClimateMachine.array_type() = Array
┌ Info: Model composition
│     physics = ClimateMachine.Atmos.AtmosPhysics{Float64,Main.##393.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.##393.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.##393.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.##393.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 = ExplicitSolverType(ClimateMachine.ODESolvers.LSRK144NiegemannDiehlBusch)
┌ Info: Starting DryRisingBubble
│     dt              = 4.21941e-01
│     timeend         =   100.00
│     number of steps = 237
│     norm(Q)         = 5.3834517091450920e+09
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:802
┌ Info: Update
│     simtime =    31.65 /   100.00
│     wallclock = 00:01:00
│     efficiency (simtime / wallclock) =   0.5201
│     wallclock end (estimated) = 00:03:12
│     norm(Q) = 5.3833059070775442e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    62.03 /   100.00
│     wallclock = 00:02:00
│     efficiency (simtime / wallclock) =   0.5132
│     wallclock end (estimated) = 00:03:14
│     norm(Q) = 5.3834364205979176e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    97.47 /   100.00
│     wallclock = 00:03:01
│     efficiency (simtime / wallclock) =   0.5375
│     wallclock end (estimated) = 00:03:06
│     norm(Q) = 5.3833158225763550e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Finished
│     norm(Q)            = 5.3833337636290865e+09
│     norm(Q) / norm(Q₀) = 9.9997809109798363e-01
│     norm(Q) - norm(Q₀) = -1.1794551600551605e+05
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:853
┌ Info: Euclidean distance
│     norm(Q - Qe)            = 7.6257985218708757e+06
│     norm(Q - Qe) / norm(Qe) = 1.4165258525337223e-03
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:869

And it successfully completes. Currently, the 14-stage LSRK method LSRK144NiegemannDiehlBusch contains the largest stability region of the low-storage methods available in ClimateMachine.jl.

Strong Stability Preserving Runge–Kutta (SSPRK) methods

Just as with the LSRK methods, the SSPRK methods are self-starting, with $\boldsymbol{Q}^{1} = \boldsymbol{q}^n$, and stage-values are of the form

\[\begin{align} \boldsymbol{Q}^{i+1} = a_{i,1} \boldsymbol{q}^n + a_{i,2} \boldsymbol{Q}^{i} + \Delta t b_i\mathcal{T}(\boldsymbol{Q}^{i}) \end{align}\]

with the value at the next step being the $(N+1)$-th stage value $\boldsymbol{q}^{n+1} = \boldsymbol{Q}^{(N+1)}$. This allows the updates to be performed with only three copies of the state vector (storing $\boldsymbol{q}^n$, $\boldsymbol{Q}^{i}$ and $\mathcal{T}(\boldsymbol{Q}^{i})$). We illustrate here the use of a SSPRK33ShuOsher method.

ode_solver = ClimateMachine.ExplicitSolverType(solver_method = SSPRK33ShuOsher)
CFL = FT(0.2)
run_simulation(ode_solver, CFL, timeend);

ClimateMachine.array_type() = Array
┌ Info: Model composition
│     physics = ClimateMachine.Atmos.AtmosPhysics{Float64,Main.##393.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.##393.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.##393.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.##393.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 = ExplicitSolverType(ClimateMachine.ODESolvers.SSPRK33ShuOsher)
┌ Info: Starting DryRisingBubble
│     dt              = 4.97265e-02
│     timeend         =   100.00
│     number of steps = 2011
│     norm(Q)         = 5.3834517091450920e+09
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:802
┌ Info: Update
│     simtime =    19.09 /   100.00
│     wallclock = 00:01:00
│     efficiency (simtime / wallclock) =   0.3170
│     wallclock end (estimated) = 00:05:15
│     norm(Q) = 5.3834106509736691e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    38.19 /   100.00
│     wallclock = 00:02:00
│     efficiency (simtime / wallclock) =   0.3178
│     wallclock end (estimated) = 00:05:14
│     norm(Q) = 5.3833309166308546e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    57.53 /   100.00
│     wallclock = 00:03:00
│     efficiency (simtime / wallclock) =   0.3193
│     wallclock end (estimated) = 00:05:13
│     norm(Q) = 5.3834443216400681e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    76.88 /   100.00
│     wallclock = 00:04:00
│     efficiency (simtime / wallclock) =   0.3199
│     wallclock end (estimated) = 00:05:12
│     norm(Q) = 5.3834604979871988e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Update
│     simtime =    96.32 /   100.00
│     wallclock = 00:05:00
│     efficiency (simtime / wallclock) =   0.3207
│     wallclock end (estimated) = 00:05:11
│     norm(Q) = 5.3833105569188194e+09
└ @ ClimateMachine.Callbacks /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Callbacks/Callbacks.jl:75
┌ Info: Finished
│     norm(Q)            = 5.3833337634148645e+09
│     norm(Q) / norm(Q₀) = 9.9997809105819091e-01
│     norm(Q) - norm(Q₀) = -1.1794573022747040e+05
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:853
┌ Info: Euclidean distance
│     norm(Q - Qe)            = 7.6257870188512690e+06
│     norm(Q - Qe) / norm(Qe) = 1.4165237157968797e-03
└ @ ClimateMachine /central/scratch/climaci/climatemachine-docs/1430/climatemachine-docs/src/Driver/Driver.jl:869

References

This page was generated using Literate.jl.