Lorenz-96 model
The classical form for the (single-layer) Lorenz-96 equations are defined as
\[\begin{align} \frac{\mathrm{d}\pmb{x}}{\mathrm{d} t} = \pmb{f}(\pmb{x}), \end{align}\]
where for each state component $i\in\{1,\cdots,n\}$,
\[\begin{align} f^i(\pmb{x}) &=-x^{i-2}x^{i-1} + x^{i-1}x^{i+1} - x^i + F \end{align}\]
such that the components of the vector $\pmb{x}$ are given by the variables $x^i$ with periodic boundary conditions, $x^0=x^n$, $x^{-1}=x^{n-1}$ and $x^{n+1}=x^{1}$. The term $F$ in the Lorenz-96 system is the forcing parameter that injects energy to the model. With the above definition for the classical Lorenz-96 equations, we define the L96-s model with additive noise (of scalar covariance) as
\[\begin{align} \frac{\mathrm{d} \pmb{x}}{\mathrm{d} t} = \pmb{f}(\pmb{x}) + s(t)\mathbf{I}_{n}\pmb{W}(t), \end{align}\]
where $\pmb{f}$ is defined as in the classical equations, $\mathbf{I}_n$ is the $n\times n$ identity matrix, $\pmb{W}(t)$ is an $n$-dimensional Wiener process and $s(t):\mathbb{R}\rightarrow \mathbb{R}$ is a measurable function of (possibly) time-varying diffusion coefficients. This model is analyzed in-depth for data assimilation twin experiments in the manuscript Grudzien et al. 2020 and further details of using the system for data assimilation benchmarks in stochastic dynamics are discussed there. The methods in the below define the model equations, the Jacobian, and the order 2.0 Taylor-Stratonovich scheme derived especially for statistically robust numerical simulation of the truth twin of the L96-s system.
Methods
DataAssimilationBenchmarks.L96.compute_α_ρ — Methodcompute_α_ρ(p::Int64)Computes auxiliary functions for the 2nd order Taylor-Stratonovich expansion. The constants α and ρ need to be computed once, only as a function of the order of truncation of the Fourier series, the argument p, for the integration method. These constants are then supplied as arguments to l96s_tay2_step! in kwargs. See l96s_tay2_step! for the interpretation and usage of these constants.
return α(p)::Float64, ρ(p)::Float64DataAssimilationBenchmarks.L96.dx_dt — Methoddx_dt(x::VecA(T), t::Float64, dx_params::ParamDict(T)) where T <: RealTime derivative for Lorenz-96 model, x is a model state of size state_dim and type VecA, t is a dummy time argument for consistency with integration methods, dx_params is of type ParamDict which is called for the forcing parameter.
Returns time derivative of the state vector
return dxDataAssimilationBenchmarks.L96.jacobian — Methodjacobian(x::VecA(T), t::Float64, dx_params::ParamDict(T)) where T <: RealComputes the Jacobian of Lorenz-96 about the state x of type VecA. The time variable t is a dummy variable for consistency with integration methods, dx_params is of type ParamDict which is called for the forcing parameter. Note that this is designed to load entries in a zeros array and return a sparse array to make a compromise between memory and computational resources.
return sparse(dxF)DataAssimilationBenchmarks.L96.l96s_tay2_step! — Methodl96s_tay2_step!(x::VecA(T), t::Float64, kwargs::StepKwargs) where T <: RealOne step of integration rule for l96 second order taylor rule The constants ρ and α are to be computed with compute_α_ρ, depending only on p, and supplied for all steps. This is the general formulation which includes, e.g., dependence on the truncation of terms in the auxilliary function C with respect to the parameter p. In general, truncation at p=1 is all that is necessary for order 2.0 convergence.
This method is derived in Grudzien, C. et al. (2020). NOTE: this Julia version still pending validation as in the manuscript
return xDataAssimilationBenchmarks.L96.mod_indx! — Methodmod_indx!(indx::Int64, dim::Int64)Auxiliary function to return state vector indices for the Lorenz-96 model, where indx is taken mod dim. Mod zero is replaced with dim for indexing in Julia state vectors.