How to run SpeedyWeather.jl
Creating a SpeedyWeather.jl simulation and running it consists conceptually of 4 steps. In contrast to many other models, these steps are bottom-up rather then top-down. There is no monolithic interface to SpeedyWeather.jl, instead all options that a user may want to adjust are chosen and live in their respective model components.
- Create a SpectralGrid which defines the grid and spectral resolution.
- Create model components and combine to a model.
- Initialize the model to create a simulation.
- Run that simulation.
In the following we will describe these steps in more detail. If you want to start with some examples first, see Examples which has easy and more advanced examples of how to set up SpeedyWeather.jl to run the simulation you want.
SpectralGrid
The life of every SpeedyWeather.jl simulation starts with a SpectralGrid object. A SpectralGrid defines the physical domain of the simulation and its discretization. This domain has to be a sphere because of the spherical harmonics, but it can have a different radius. The discretization is for spectral, grid-point space and the vertical as this determines the size of many arrays for preallocation, for which als the number format is essential to know. That's why SpectralGrid is the beginning of every SpeedyWeather.jl simulation and that is why it has to be passed on to (most) model components.
The default SpectralGrid is
using SpeedyWeather
spectral_grid = SpectralGrid()SpectralGrid:
├ Spectral: T31 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid: 48-ring OctahedralGaussianGrid{Float32}, 3168 grid points
├ Resolution: 401km (average)
├ Vertical: 8-layer SigmaCoordinates
└ Device: CPU using ArrayYou can also get the help prompt by typing ?SpectralGrid. Let's explain the details: The spectral resolution is T31, so the largest wavenumber in spectral space is 31, and all the complex spherical harmonic coefficients of a given 2D field (see Spherical Harmonic Transform) are stored in a LowerTriangularMatrix in the number format Float32. The radius of the sphere is 6371km by default, which is the radius of the Earth.
This spectral resolution is combined with an octahedral Gaussian grid of 3168 grid points. This grid has 48 latitude rings, 20 longitude points around the poles and up to 96 longitude points around the Equator. Data on that grid is also stored in Float32. The resolution is therefore on average about 400km. In the vertical 8 levels are used, using Sigma coordinates.
The resolution of a SpeedyWeather.jl simulation is adjusted using the trunc argument, this defines the spectral resolution and the grid resolution is automatically adjusted to keep the aliasing between spectral and grid-point space constant (see Matching spectral and grid resolution).
spectral_grid = SpectralGrid(trunc=85)SpectralGrid:
├ Spectral: T85 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid: 128-ring OctahedralGaussianGrid{Float32}, 18688 grid points
├ Resolution: 165km (average)
├ Vertical: 8-layer SigmaCoordinates
└ Device: CPU using ArrayTypical values are 31, 42, 63, 85, 127, 170, ... although you can technically use any integer, see Available horizontal resolutions for details. Now with T85 (which is a common notation for trunc=85) the grid is of higher resolution too. You may play with the dealiasing factor, a larger factor increases the grid resolution that is matched with a given spectral resolution. You don't choose the resolution of the grid directly, but using the Grid argument you can change its type (see Grids)
spectral_grid = SpectralGrid(trunc=85, dealiasing=3, Grid=HEALPixGrid)SpectralGrid:
├ Spectral: T85 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid: 179-ring HEALPixGrid{Float32}, 24300 grid points
├ Resolution: 145km (average)
├ Vertical: 8-layer SigmaCoordinates
└ Device: CPU using ArrayVertical coordinates and resolution
The number of vertical layers or levels (we use both terms often interchangeably) is determined through the nlayers argument. Especially for the BarotropicModel and the ShallowWaterModel you want to set this to
spectral_grid = SpectralGrid(nlayers=1)SpectralGrid:
├ Spectral: T31 LowerTriangularMatrix{Complex{Float32}}, radius = 6.371e6 m
├ Grid: 48-ring OctahedralGaussianGrid{Float32}, 3168 grid points
├ Resolution: 401km (average)
├ Vertical: 1-layer SigmaCoordinates
└ Device: CPU using ArrayFor a single vertical level the type of the vertical coordinates does not matter, but in general you can change the spacing of the sigma coordinates which have to be discretized in $[0, 1]$
vertical_coordinates = SigmaCoordinates(0:0.2:1)5-layer SigmaCoordinates
├─ 0.0000 k = 0.5
│× 0.1000 k = 1
├─ 0.2000 k = 1.5
│× 0.3000 k = 2
├─ 0.4000 k = 2.5
│× 0.5000 k = 3
├─ 0.6000 k = 3.5
│× 0.7000 k = 4
├─ 0.8000 k = 4.5
│× 0.9000 k = 5
└─ 1.0000 k = 5.5These are regularly spaced Sigma coordinates, defined through their half levels. The cell centers or called full levels are marked with an ×.
Creating model components
SpeedyWeather.jl deliberately avoids the namelists that are the main user interface to many old school models. Instead, we employ a modular approach whereby every non-default model component is created with its respective parameters to finally build the whole model from these components.
If you know which components you want to adjust you would create them right away, however, a new user might first want to know which components there are, so let's flip the logic around and assume you know you want to run a ShallowWaterModel. You can create a default ShallowWaterModel like so and inspect its components
model = ShallowWaterModel(spectral_grid)ShallowWaterModel <: ShallowWater
├ spectral_grid: SpectralGrid
├ device_setup: SpeedyWeather.DeviceSetup{CPU, DataType}
├ geometry: Geometry{Float32, OctahedralGaussianGrid}
├ planet: Earth{Float32}
├ atmosphere: EarthAtmosphere{Float32}
├ coriolis: Coriolis{Float32}
├ orography: EarthOrography{Float32, OctahedralGaussianGrid{Float32}}
├ forcing: NoForcing
├ drag: NoDrag
├ particle_advection: NoParticleAdvection
├ initial_conditions: InitialConditions{ZonalJet, ZeroInitially, ZeroInitially, ZeroInitially}
├ random_process: NoRandomProcess
├ time_stepping: Leapfrog{Float32}
├ spectral_transform: SpectralTransform{Float32, Array, Vector{Float32}, Vector{ComplexF32}, Vector{Int64}, Matrix{ComplexF32}, Array{ComplexF32, 3}, LowerTriangularMatrix{Float32}, LowerTriangularArray{Float32, 2, Matrix{Float32}}}
├ implicit: ImplicitShallowWater{Float32}
├ horizontal_diffusion: HyperDiffusion{Float32, Matrix{Float32}}
├ output: NetCDFOutput{FullGaussianGrid{Float32}, FullGaussianArray{Float32, 2, Matrix{Float32}}, AnvilInterpolator{Float32, OctahedralGaussianGrid}}
├ callbacks: Dict{Symbol, SpeedyWeather.AbstractCallback}
└ feedback: FeedbackSo by default the ShallowWaterModel uses a Leapfrog time_stepping, which you can inspect like so
model.time_steppingLeapfrog{Float32} <: SpeedyWeather.AbstractTimeStepper
├ trunc::Int64 = 31
├ nsteps::Int64 = 2
├ Δt_at_T31::Second = 1800 seconds
├ radius::Float32 = 6.371e6
├ adjust_with_output::Bool = true
├ robert_filter::Float32 = 0.1
├ williams_filter::Float32 = 0.53
├ Δt_millisec::Dates.Millisecond = 1800000 milliseconds
├ Δt_sec::Float32 = 1800.0
└ Δt::Float32 = 0.0002825302Model components often contain parameters from the SpectralGrid as they are needed to determine the size of arrays and other internal reasons. You should, in most cases, just ignore those. But the Leapfrog time stepper comes with Δt_at_T31 which is the parameter used to scale the time step automatically. This means at a spectral resolution of T31 it would use 30min steps, at T63 it would be ~half that, 15min, etc. Meaning that if you want to have a shorter or longer time step you can create a new Leapfrog time stepper. All time inputs are supposed to be given with the help of Dates (e.g. Minute(), Hour(), ...). But remember that (almost) every model component depends on a SpectralGrid as first argument.
spectral_grid = SpectralGrid(trunc=63, nlayers=1)
time_stepping = Leapfrog(spectral_grid, Δt_at_T31=Minute(15))Leapfrog{Float32} <: SpeedyWeather.AbstractTimeStepper
├ trunc::Int64 = 63
├ nsteps::Int64 = 2
├ Δt_at_T31::Second = 900 seconds
├ radius::Float32 = 6.371e6
├ adjust_with_output::Bool = true
├ robert_filter::Float32 = 0.1
├ williams_filter::Float32 = 0.53
├ Δt_millisec::Dates.Millisecond = 450000 milliseconds
├ Δt_sec::Float32 = 450.0
└ Δt::Float32 = 7.063255e-5The actual time step at the given resolution (here T63) is then Δt_sec, there's also Δt which is a scaled time step used internally, because SpeedyWeather.jl scales the equations with the radius of the Earth, but this is largely hidden (except here) from the user. With this new Leapfrog time stepper constructed we can create a model by passing on the components (they are keyword arguments so either use ; time_stepping for which the naming must match, or time_stepping=my_time_stepping with any name)
model = ShallowWaterModel(spectral_grid; time_stepping)ShallowWaterModel <: ShallowWater
├ spectral_grid: SpectralGrid
├ device_setup: SpeedyWeather.DeviceSetup{CPU, DataType}
├ geometry: Geometry{Float32, OctahedralGaussianGrid}
├ planet: Earth{Float32}
├ atmosphere: EarthAtmosphere{Float32}
├ coriolis: Coriolis{Float32}
├ orography: EarthOrography{Float32, OctahedralGaussianGrid{Float32}}
├ forcing: NoForcing
├ drag: NoDrag
├ particle_advection: NoParticleAdvection
├ initial_conditions: InitialConditions{ZonalJet, ZeroInitially, ZeroInitially, ZeroInitially}
├ random_process: NoRandomProcess
├ time_stepping: Leapfrog{Float32}
├ spectral_transform: SpectralTransform{Float32, Array, Vector{Float32}, Vector{ComplexF32}, Vector{Int64}, Matrix{ComplexF32}, Array{ComplexF32, 3}, LowerTriangularMatrix{Float32}, LowerTriangularArray{Float32, 2, Matrix{Float32}}}
├ implicit: ImplicitShallowWater{Float32}
├ horizontal_diffusion: HyperDiffusion{Float32, Matrix{Float32}}
├ output: NetCDFOutput{FullGaussianGrid{Float32}, FullGaussianArray{Float32, 2, Matrix{Float32}}, AnvilInterpolator{Float32, OctahedralGaussianGrid}}
├ callbacks: Dict{Symbol, SpeedyWeather.AbstractCallback}
└ feedback: FeedbackThis logic continues for all model components, see Examples. All model components are also subtype (i.e. <:) of some abstract supertype, this feature can be made use of to redefine not just constant parameters of existing model components but also to define new ones. This is more elaborated in Extending SpeedyWeather.
A model in SpeedyWeather.jl is understood as a collection of model components that roughly belong into one of three groups.
- Components (numerics, dynamics, or physics) that have parameters, possibly contain some data (immutable, precomputed) and some functions associated with them. For example, a forcing term may be defined through some parameters, but also require precomputed arrays, or data to be loaded from file and a function that instructs how to calculate this forcing on every time step.
- Components that are merely a collection of parameters that conceptually belong together. For example,
Earthis anAbstractPlanetthat contains all the orbital parameters important for weather and climate (rotation, gravity, etc). - Components for output purposes. For example,
OutputWritercontrols the NetCDF output, andFeedbackcontrols the information printed to the REPL and to file.
Creating a model
SpeedyWeather.jl currently includes 4 different models:
BarotropicModelfor the 2D barotropic vorticity equation,ShallowWaterModelfor the 2D shallow water equations,PrimitiveDryModelfor the 3D primitive equations without humidity, andPrimitiveWetModelfor the 3D primitive equations with humidity.
Overall, all above-mentioned models are kept quite similar, but there are still fundamental differences arising from the different equations. For example, the barotropic and shallow water models do not have any physical parameterizations. Conceptually you construct these different models with
spectral_grid = SpectralGrid(trunc=..., ...)
component1 = SomeComponent(spectral_grid, parameter1=..., ...)
component2 = SomeOtherComponent(spectral_grid, parameter2=..., ...)
model = BarotropicModel(spectral_grid; all_other_components..., ...)or model = ShallowWaterModel(spectral_grid; ...), etc.
Model initialization
In the previous section the model was created, but this is conceptually just gathering all its components together. However, many components need to be initialized. This step is used to precompute arrays, load necessary data from file or to communicate those between components. Furthermore, prognostic and diagnostic variables are allocated. It is (almost) all that needs to be done before the model can be run (exception being the output initialization). Many model components have a initialize! function associated with them that it executed here. However, from a user perspective all that needs to be done here is
simulation = initialize!(model)Simulation{ShallowWaterModel}
├ prognostic_variables::PrognosticVariables{...}
├ diagnostic_variables::DiagnosticVariables{...}
└ model::ShallowWaterModel{...}and we have initialized the ShallowWaterModel we have defined earlier. As initialize!(model) also initializes the prognostic (and diagnostic) variables, it also initializes the clock in simulation.prognostic_variables.clock. To initialize with a specific time, do
simulation = initialize!(model, time=DateTime(2020,5,1))
simulation.prognostic_variables.clock.time2020-05-01T00:00:00to set the time to 1st May, 2020 (but you can also do that manually). This time is used by components that depend on time, e.g. the solar zenith angle calculation.
After this step you can continue to tweak your model setup but note that some model components are immutable, or that your changes may not be propagated to other model components that rely on it. But you can, for example, change the output time step like so
simulation.model.output.output_dt = Second(3600)3600 secondsNow, if there's output, it will be every hour. Furthermore the initial conditions can be set with the initial_conditions model component which are then set during initialize!(::AbstractModel), but you can also change them now, before the model runs
simulation.prognostic_variables.vor[1][1, 1] = 00So with this we have set the zero mode of vorticity of the first (and only) layer in the shallow water to zero. Because the leapfrogging is a 2-step time stepping scheme we set here the first. As it is often tricky to set the initial conditions in spectral space, it is generally advised to do so through the initial_conditions model component.
Run a simulation
After creating a model, initializing it to a simulation, all that is left is to run the simulation.
run!(simulation) Surface relative vorticity [1/s]
┌────────────────────────────────────────────────────────────┐0.0001
90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
-90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
└────────────────────────────────────────────────────────────┘ -9e⁻⁵
0 ˚E 360 By default this runs for 10 days without output and returns a unicode plot of surface relative vorticity (see Shallow water with mountains for more on this). Now period and output are the only options to change, so with
run!(simulation, period=Day(5), output=true) Surface relative vorticity [1/s]
┌────────────────────────────────────────────────────────────┐ 7e⁻⁵
90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ ┌──┐
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
˚N │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
│▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ │▄▄│
-90 │▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄│ └──┘
└────────────────────────────────────────────────────────────┘ -8e⁻⁵
0 ˚E 360 You would continue this simulation (the previous run! call already integrated 10 days!) for another 5 days and storing default NetCDF output.