API¶
Base Model Class¶
This is the base class from which all other models inherit. All of these initialization arguments are available to all of the other model types. This class is not called directly.
- class pyqg.Model(nz=1, nx=64, ny=None, L=1000000.0, W=None, dt=7200.0, twrite=1000.0, tmax=1576800000.0, tavestart=315360000.0, taveint=86400.0, useAB2=False, rek=5.787e-07, filterfac=23.6, f=None, g=9.81, q_parameterization=None, uv_parameterization=None, parameterization=None, diagnostics_list='all', ntd=1, log_level=1, logfile=None)¶
A generic pseudo-spectral inversion model.
- Attributes:
- nx, nyint
Number of real space grid points in the x, y directions (cython)
- nk, nlint
Number of spectral space grid points in the k, l directions (cython)
- nzint
Number of vertical levels (cython)
- kk, llreal array
Zonal and meridional wavenumbers (nk) (cython)
- areal array
inversion matrix (nk, nk, nl, nk) (cython)
- qreal array
Potential vorticity in real space (nz, ny, nx) (cython)
- qhcomplex array
Potential vorticity in spectral space (nk, nl, nk) (cython)
- phcomplex array
Streamfunction in spectral space (nk, nl, nk) (cython)
- u, vreal array
Zonal and meridional velocity anomalies in real space (nz, ny, nx) (cython)
- Ubgreal array
Background zonal velocity (nk) (cython)
- Qyreal array
Background potential vorticity gradient (nk) (cython)
- ufull, vfullreal arrays
Zonal and meridional full velocities in real space (nz, ny, nx) (cython)
- uh, vhcomplex arrays
Velocity anomaly components in spectral space (nk, nl, nk) (cython)
- rekfloat
Linear drag in lower layer (cython)
- tfloat
Model time (cython)
- tcint
Model timestep (cython)
- dtfloat
Numerical timestep (cython)
- L, Wfloat
Domain length in x and y directions
- filterfacfloat
Amplitdue of the spectral spherical filter
- twriteint
Interval for cfl writeout (units: number of timesteps)
- tmaxfloat
Total time of integration (units: model time)
- tavestartfloat
Start time for averaging (units: model time)
- tsnapstartfloat
Start time for snapshot writeout (units: model time)
- taveintfloat
Time interval for accumulation of diagnostic averages. (units: model time)
- tsnapintfloat
Time interval for snapshots (units: model time)
- ntdint
Number of threads to use. Should not exceed the number of cores on your machine.
- pmodesreal array
Vertical pressure modes (unitless)
- radiireal array
Deformation radii (units: model length)
- q_parameterizationfunction or pyqg.Parameterization
Optional
Parameterization
object or function which takes the model as input and returns anumpy
array of shape(nz, ny, nx)
to be added to \(\partial_t q\) before stepping forward. This can be used to implement subgrid forcing parameterizations.- uv_parameterizationfunction or pyqg.Parameterization
Optional
Parameterization
object or function which takes the model as input and returns a tuple of twonumpy
arrays, each of shape(nz, ny, nx)
, to be added to the zonal and meridional velocity derivatives (respectively) at each timestep (by adding their curl to \(\partial_t q\)). This can also be used to implemented subgrid forcing parameterizations, but expressed in terms of velocity rather than potential vorticity.
Note
All of the test cases use
nx==ny
. Expect bugs if you choose these parameters to be different.Note
All time intervals will be rounded to nearest dt interval.
- Parameters:
- nxint
Number of grid points in the x direction.
- nyint
Number of grid points in the y direction (default: nx).
- Lnumber
Domain length in x direction. Units: meters.
- W
Domain width in y direction. Units: meters (default: L).
- reknumber
linear drag in lower layer. Units: seconds ^{-1}.
- filterfacnumber
amplitdue of the spectral spherical filter (originally 18.4, later changed to 23.6).
- dtnumber
Numerical timstep. Units: seconds.
- twriteint
Interval for cfl writeout. Units: number of timesteps.
- tmaxnumber
Total time of integration. Units: seconds.
- tavestartnumber
Start time for averaging. Units: seconds.
- tsnapstartnumber
Start time for snapshot writeout. Units: seconds.
- taveintnumber
Time interval for accumulation of diagnostic averages. Units: seconds. (For performance purposes, averaging does not have to occur every timestep)
- tsnapintnumber
Time interval for snapshots. Units: seconds.
- ntdint
Number of threads to use. Should not exceed the number of cores on your machine.
- q_parameterizationfunction or pyqg.Parameterization
Optional
Parameterization
object or function which takes the model as input and returns anumpy
array of shape(nz, ny, nx)
to be added to \(\partial_t q\) before stepping forward. This can be used to implement subgrid forcing parameterizations.- uv_parameterizationfunction or pyqg.Parameterization
Optional
Parameterization
object or function which takes the model as input and returns a tuple of twonumpy
arrays, each of shape(nz, ny, nx)
, to be added to the zonal and meridional velocity derivatives (respectively) at each timestep (by adding their curl to \(\partial_t q\)). This can also be used to implemented subgrid forcing parameterizations, but expressed in terms of velocity rather than potential vorticity.- parameterizationpyqg.Parameterization
An explicit
Parameterization
object representing either a velocity or potential vorticity parameterization, whose type will be inferred.
- describe_diagnostics()¶
Print a human-readable summary of the available diagnostics.
- modal_projection(p, forward=True)¶
Performs a field p into modal amplitudes pn using the basis [pmodes]. The inverse transform calculates p from pn
- property parameterization¶
Return the model’s parameterization if present (either in terms of PV or velocity, warning if there are both).
- Returns:
- parameterizationpyqg.Parameterization or function
- run()¶
Run the model forward without stopping until the end.
- run_with_snapshots(tsnapstart=0.0, tsnapint=432000.0)¶
Run the model forward, yielding to user code at specified intervals.
- Parameters:
- tsnapstartint
The timestep at which to begin yielding.
- tstapintint
The interval at which to yield.
- spec_var(ph)¶
compute variance of p from Fourier coefficients ph
- stability_analysis(bottom_friction=False)¶
- Performs the baroclinic linear instability analysis given
given the base state velocity :math: (U, V) and the stretching matrix :math: S:
\[A \Phi = \omega B \Phi,\]where
\[A = B (U k + V l) + I (k Q_y - l Q_x) + 1j \delta_{N N} r_{ek} I \kappa^2\]where \(\delta_{N N} = [0,0,\dots,0,1] ,\)
and
\[B = S - I \kappa^2 .\]The eigenstructure is
\[\Phi\]and the eigenvalue is
\[`\omega`\]The growth rate is Im\(\{\omega\}\).
- Parameters:
- bottom_friction: optional inclusion linear bottom drag
in the linear stability calculation (default is False, as if :math: r_{ek} = 0)
- Returns:
- omega: complex array
The eigenvalues with largest complex part (units: inverse model time)
- phi: complex array
The eigenvectors associated associated with omega (unitless)
- to_dataset()¶
Convert outputs from model to an xarray dataset
- Returns:
- dsxarray.Dataset
- vertical_modes()¶
Calculate standard vertical modes. Simply the eigenvectors of the stretching matrix S
Specific Model Types¶
These are the actual models which are run.
- class pyqg.QGModel(beta=1.5e-11, rd=15000.0, delta=0.25, H1=500, U1=0.025, U2=0.0, **kwargs)¶
Two layer quasigeostrophic model.
This model is meant to representflows driven by baroclinic instabilty of a base-state shear \(U_1-U_2\). The upper and lower layer potential vorticity anomalies \(q_1\) and \(q_2\) are
\[\begin{split}q_1 &= \nabla^2\psi_1 + F_1(\psi_2 - \psi_1) \\ q_2 &= \nabla^2\psi_2 + F_2(\psi_1 - \psi_2)\end{split}\]with
\[\begin{split}F_1 &\equiv \frac{k_d^2}{1 + \delta^2} \\ F_2 &\equiv \delta F_1 \ .\end{split}\]The layer depth ratio is given by \(\delta = H_1 / H_2\). The total depth is \(H = H_1 + H_2\).
The background potential vorticity gradients are
\[\begin{split}\beta_1 &= \beta + F_1(U_1 - U_2) \\ \beta_2 &= \beta - F_2( U_1 - U_2) \ .\end{split}\]The evolution equations for \(q_1\) and \(q_2\) are
\[\begin{split}\partial_t\,{q_1} + J(\psi_1\,, q_1) + \beta_1\, {\psi_1}_x &= \text{ssd} \\ \partial_t\,{q_2} + J(\psi_2\,, q_2)+ \beta_2\, {\psi_2}_x &= -r_{ek}\nabla^2 \psi_2 + \text{ssd}\,.\end{split}\]where ssd represents small-scale dissipation and \(r_{ek}\) is the Ekman friction parameter.
- Parameters:
- betanumber
Gradient of coriolis parameter. Units: meters ^{-1} seconds ^{-1}
- reknumber
Linear drag in lower layer. Units: seconds ^{-1}
- rdnumber
Deformation radius. Units: meters.
- deltanumber
Layer thickness ratio (H1/H2)
- U1number
Upper layer flow. Units: meters seconds ^{-1}
- U2number
Lower layer flow. Units: meters seconds ^{-1}
- set_U1U2(U1, U2)¶
Set background zonal flow.
- Parameters:
- U1number
Upper layer flow. Units: meters seconds ^{-1}
- U2number
Lower layer flow. Units: meters seconds ^{-1}
- set_q1q2(q1, q2, check=False)¶
Set upper and lower layer PV anomalies.
- Parameters:
- q1array-like
Upper layer PV anomaly in spatial coordinates.
- q1array-like
Lower layer PV anomaly in spatial coordinates.
- class pyqg.LayeredModel(beta=1.5e-11, nz=3, rd=15000.0, f=0.0001236812857687059, H=None, U=None, V=None, rho=None, delta=None, **kwargs)¶
Layered quasigeostrophic model.
This model is meant to represent flows driven by baroclinic instabilty of a base-state shear. The potential vorticity anomalies qi are related to the streamfunction psii through
\[ \begin{align}\begin{aligned}{q_i} = \nabla^2\psi_i + \frac{f_0^2}{H_i} \left(\frac{\psi_{i-1}- \psi_i}{g'_{i-1}}- \frac{\psi_{i}-\psi_{i+1}}{g'_{i}}\right)\,, \qquad i = 2,\textsf{N}-1\,,\\{q_1} = \nabla^2\psi_1 + \frac{f_0^2}{H_1} \left(\frac{\psi_{2}- \psi_1}{g'_{1}}\right)\,, \qquad i =1\,,\\{q_\textsf{N}} = \nabla^2\psi_\textsf{N} + \frac{f_0^2}{H_\textsf{N}} \left(\frac{\psi_{\textsf{N}-1}- \psi_\textsf{N}}{g'_{\textsf{N}}}\right) + \frac{f_0}{H_\textsf{N}}h_b\,, \qquad i =\textsf{N}\,,\end{aligned}\end{align} \]where the reduced gravity, or buoyancy jump, is
\[g'_i \equiv g \frac{\rho_{i+1}-\rho_i}{\rho_i}\,.\]The evolution equations are
\[\,{q_{i}}_t + \mathsf{J}\left(\psi_i\,, q_i\right) + \textsf{Q}_y {\psi_i}_x - \textsf{Q}_x {\psi_i}_y = \text{ssd} - r_{ek} \delta_{i\textsf{N}} \nabla^2 \psi_i\,, \qquad i = 1,\textsf{N}\,,\]where the mean potential vorticy gradients are
\[\textsf{Q}_x = \textsf{S}\textsf{V}\,,\]and
\[\textsf{Q}_y = \beta\,\textsf{I} - \textsf{S}\textsf{U}\,\,,\]where S is the stretching matrix, I is the identity matrix, and the background velocity is
\(\vec{\textsf{V}}(z) = \left(\textsf{U},\textsf{V}\right)\).
- Parameters:
- nzinteger number
Number of layers (> 1)
- betanumber
Gradient of coriolis parameter. Units: meters ^{-1} seconds ^{-1}
- rdnumber
Deformation radius. Units: meters. Only necessary for the two-layer (nz=2) case.
- deltanumber
Layer thickness ratio (H1/H2). Only necessary for the two-layer (nz=2) case. Unitless.
- Ulist of size nz
Base state zonal velocity. Units: meters seconds ^{-1}
- Varray of size nz
Base state meridional velocity. Units: meters seconds ^{-1}
- Harray of size nz
Layer thickness. Units: meters
- rho: array of size nz.
Layer density. Units: kilograms meters ^{-3}
- class pyqg.BTModel(beta=0.0, rd=0.0, H=1.0, U=0.0, **kwargs)¶
Single-layer (barotropic) quasigeostrophic model. This class can represent both pure two-dimensional flow and also single reduced-gravity layers with deformation radius
rd
.The equivalent-barotropic quasigeostrophic evolution equations is
\[\partial_t q + J(\psi, q ) + \beta \psi_x = \text{ssd}\]The potential vorticity anomaly is
\[q = \nabla^2 \psi - \kappa_d^2 \psi\]- Parameters:
- betanumber, optional
Gradient of coriolis parameter. Units: meters ^{-1} seconds ^{-1}
- rdnumber, optional
Deformation radius. Units: meters.
- Unumber, optional
Upper layer flow. Units: meters seconds ^{-1}.
- set_U(U)¶
Set background zonal flow.
- Parameters:
- Unumber
Upper layer flow. Units: meters seconds ^{-1}.
- class pyqg.SQGModel(beta=0.0, Nb=1.0, f_0=1.0, H=1.0, U=0.0, **kwargs)¶
Surface quasigeostrophic model.
- Parameters:
- betanumber
Gradient of coriolis parameter. Units: meters ^{-1} seconds ^{-1}
- Nbnumber
Buoyancy frequency. Units: seconds ^{-1}.
- Unumber
Background zonal flow. Units: meters seconds ^{-1}.
- set_U(U)¶
Set background zonal flow
Lagrangian Particles¶
- class pyqg.LagrangianParticleArray2D(x0, y0, periodic_in_x=False, periodic_in_y=False, xmin=- inf, xmax=inf, ymin=- inf, ymax=inf, particle_dtype='f8')¶
A class for keeping track of a set of lagrangian particles in two-dimensional space. Tries to be fast.
- Parameters:
- x0, y0array-like
Two arrays (same size) representing the particle initial positions.
- periodic_in_xbool
Whether the domain wraps in the x direction.
- periodic_in_ybool
Whether the domain ‘wraps’ in the y direction.
- xmin, xmaxnumbers
Maximum and minimum values of x coordinate
- ymin, ymaxnumbers
Maximum and minimum values of y coordinate
- particle_dtypedtype
Data type to use for particles
- step_forward_with_function(uv0fun, uv1fun, dt)¶
Advance particles using a function to determine u and v.
- Parameters:
- uv0funfunction
Called like
uv0fun(x,y)
. Should return the velocity field u, v at time t.- uv1fun(x,y)function
Called like
uv1fun(x,y)
. Should return the velocity field u, v at time t + dt.- dtnumber
Timestep.
- class pyqg.GriddedLagrangianParticleArray2D(x0, y0, Nx, Ny, grid_type='A', **kwargs)¶
Lagrangian particles with velocities given on a regular cartesian grid.
- Parameters:
- x0, y0array-like
Two arrays (same size) representing the particle initial positions.
- Nx, Ny: int
Number of grid points in the x and y directions
- grid_type: {‘A’}
Arakawa grid type specifying velocity positions.
- interpolate_gridded_scalar(x, y, c, order=1, pad=1, offset=0)¶
Interpolate gridded scalar C to points x,y.
- Parameters:
- x, yarray-like
Points at which to interpolate
- carray-like
The scalar, assumed to be defined on the grid.
- orderint
Order of interpolation
- padint
Number of pad cells added
- offsetint
???
- Returns:
- ciarray-like
The interpolated scalar
- step_forward_with_gridded_uv(U0, V0, U1, V1, dt, order=1)¶
Advance particles using a gridded velocity field. Because of the Runga-Kutta timestepping, we need two velocity fields at different times.
- Parameters:
- U0, V0array-like
Gridded velocity fields at time t - dt.
- U1, V1array-like
Gridded velocity fields at time t.
- dtnumber
Timestep.
- orderint
Order of interpolation.
Diagnostic Tools¶
Utility functions for pyqg model data.
- pyqg.diagnostic_tools.calc_ispec(model, _var_dens, averaging=True, truncate=True, nd_wavenumber=False, nfactor=1)¶
Compute isotropic spectrum phr from 2D spectrum of variable signal2d such that signal2d.var() = phr.sum() * (kr[1] - kr[0]).
- Parameters:
- modelpyqg.Model instance
The model object from which var_dens originates
- var_denssquared modulus of fourier coefficients like this:
np.abs(signal2d_fft)**2/m.M**2
- averaging: If True, spectral density is estimated with averaging over circles,
otherwise summation is used and Parseval identity holds
- truncate: If True, maximum wavenumber corresponds to inner circle in Fourier space,
otherwise - outer circle
- nd_wavenumber: If True, wavenumber is nondimensional:
minimum wavenumber is 1 and corresponds to domain length/width, otherwise - wavenumber is dimensional [m^-1]
- nfactor: width of the bin in sqrt(dk^2+dl^2) units
- Returns:
- krarray
isotropic wavenumber
- phrarray
isotropic spectrum
- pyqg.diagnostic_tools.diagnostic_differences(m1, m2, reduction='rmse', instantaneous=False)¶
Compute a dictionary of differences in the diagnostics of two models at possibly different resolutions (e.g. for quantifying the effects of parameterizations). Applies normalization/isotropization to certain diagnostics before comparing them and skips others. Also computes differences for each vertical layer separately.
- Parameters:
- m1pyqg.Model instance
The first model to compare
- m2pyqg.Model instance
The second model to compare
- reductionstring or function
A function that takes two arrays of diagnostics and computes a distance metric. Defaults to the root mean squared difference (‘rmse’).
- instantaneousboolean
If true, compute difference metrics for the instantaneous values of a diagnostic, rather than its time average. Defaults to false.
- Returns:
- diffsdict
A dictionary of diagnostic name => distance. If the diagnostic is defined over multiple layers, separate keys are included with an appended z index.
- pyqg.diagnostic_tools.diagnostic_similarities(model, target, baseline, **kw)¶
Like diagnostic_differences, but returning a dictionary of similarity scores between negative infinity and 1 which quantify how much closer the diagnostics of a given model are to a target with respect to a baseline. Scores approach 1 when the distance between the model and the target is small compared to the baseline and are negative when that distance is greater.
- Parameters:
- modelpyqg.Model instance
The model for which we want to compute similiarity scores (e.g. a parameterized low resolution model)
- targetpyqg.Model instance
The target model (e.g. a high resolution model)
- baselinepyqg.Model instance
The baseline against which we check for improvement or degradation (e.g. an unparameterized low resolution model)
- Returns:
- simsdict
A dictionary of diagnostic name => similarity. If the diagnostic is defined over multiple layers, separate keys are included with an appended z index.
- pyqg.diagnostic_tools.spec_sum(ph2)¶
Compute total spectral sum of the real spectral quantity``ph^2``.
- Parameters:
- modelpyqg.Model instance
The model object from which ph originates
- ph2real array
The field on which to compute the sum
- Returns:
- var_densfloat
The sum of ph2
- pyqg.diagnostic_tools.spec_var(model, ph)¶
Compute variance of
p
from Fourier coefficientsph
.- Parameters:
- modelpyqg.Model instance
The model object from which ph originates
- phcomplex array
The field on which to compute the variance
- Returns:
- var_densfloat
The variance of ph
Parameterizations¶
- class pyqg.Parameterization¶
A generic class representing a subgrid parameterization. Inherit from this class, \(UVParameterization\), or \(QParameterization\) to define a new parameterization.
- abstract __call__(m)¶
Call the parameterization given a pyqg.Model. Override this function in the subclass when defining a new parameterization.
- Parameters:
- mModel
The model for which we are evaluating the parameterization.
- Returns:
- forcingreal array or tuple
The forcing associated with the model. If the model has been initialized with this parameterization as its
q_parameterization
, this should be an array of shape(nz, ny, nx)
. Foruv_parameterization
, this should be a tuple of two such arrays or a single array of shape(2, nz, ny, nx)
.
- abstract property parameterization_type¶
Whether the parameterization applies to velocity (in which case this property should return
"uv_parameterization"
) or potential vorticity (in which case this property should return"q_parameterization"
). If you inherit fromUVParameterization
orQParameterization
, this will be defined automatically.- Returns:
- parameterization_typestring
Either
"uv_parameterization"
or"q_parameterization"
, depending on how the output should be interpreted.
- __add__(other)¶
Add two parameterizations (returning a new object).
- Parameters:
- otherParameterization
The parameterization to add to this one.
- Returns:
- sumParameterization
The sum of the two parameterizations.
- __mul__(constant)¶
Multiply a parameterization by a constant (returning a new object).
- Parameters:
- constantnumber
Multiplicative factor for scaling the parameterization.
- Returns:
- productParameterization
The parameterization times the constant.
- class pyqg.parameterizations.UVParameterization¶
A generic class representing a subgrid parameterization in terms of velocity. Inherit from this to define a new velocity parameterization.
- class pyqg.parameterizations.QParameterization¶
A generic class representing a subgrid parameterization in terms of potential vorticity. Inherit from this to define a new potential vorticity parameterization.
- class pyqg.parameterizations.Smagorinsky(constant=0.1)¶
Velocity parameterization from Smagorinsky 1963.
This parameterization assumes that due to subgrid stress, there is an effective eddy viscosity
\[\nu = (C_S \Delta)^2 \sqrt{2(S_{x,x}^2 + S_{y,y}^2 + 2S_{x,y}^2)}\]which leads to updated velocity tendencies \(\Pi_{i}, i \in \{1,2\}\) corresponding to \(x\) and \(y\) respectively (equation is the same in each layer):
\[\Pi_{i} = 2 \partial_i(\nu S_{i,i}) + \partial_{2-i}(\nu S_{i,2-i})\]where \(C_S\) is a tunable Smagorinsky constant, \(\Delta\) is the grid spacing, and
\[S_{i,j} = \frac{1}{2}(\partial_i \mathbf{u}_j + \partial_j \mathbf{u}_i)\]- Parameters:
- constantnumber
Smagorinsky constant \(C_S\). Defaults to 0.1.
- class pyqg.parameterizations.BackscatterBiharmonic(smag_constant=0.08, back_constant=0.99, eps=1e-32)¶
PV parameterization based on Jansen and Held 2014 and Jansen et al. 2015 (adapted by Pavel Perezhogin). Assumes that a configurable fraction of Smagorinsky dissipation is scattered back to larger scales in an energetically consistent way.
- Parameters:
- smag_constantnumber
Smagorinsky constant \(C_S\) for the dissipative model. Defaults to 0.08.
- back_constantnumber
Backscatter constant \(C_B\) describing the fraction of Smagorinsky-dissipated energy which should be scattered back to larger scales. Defaults to 0.99. Normally should be less than 1, but larger values may still be stable, e.g. due to additional dissipation in the model from numerical filtering.
- epsnumber
Small constant to add to the denominator of the backscatter formula to prevent division by zero errors. Defaults to 1e-32.
- class pyqg.parameterizations.ZannaBolton2020(constant=- 46761284)¶
Velocity parameterization derived from equation discovery by Zanna and Bolton 2020 (Eq. 6).
- Parameters:
- constantnumber
Scaling constant \(\kappa_{BC}\). Units: meters ^{-2}. Defaults to \(\approx -4.68 \times 10^7\), a value obtained by empirically minimizing squared error with respect to the subgrid forcing that results from applying the filtering method of Guan et al. 2022 to a two-layer QGModel with default parameters.
- class pyqg.parameterizations.RingForcing(k_in_forc=0, k_out_forc=0, mag_noise_forc=0, layers='all')¶
Stochastically force PV in spectral space on a ring associated with a band of given wavenumbers from Uchida et al. (2023).
This parametrization introduces a vertically uniform stochastic forcing decorrelated in time by inverse Fourier transforming white noise in the wavenumber domain
\[\hat{w}(t,k^y,k^x) = a(t,k^y,k^x) + ib(t,k^y,k^x)\]in the wavenumber band between \(k_{in}<\sqrt{{k^x}^2 + {k^y}^2}<k_{out}\) and is zero outside of this band. \(k^x, k^y\) are the zonal and meridional wavenumbers and \(a, b\) are Gaussian random variables.
After taking the inverse Fourier transform, the horizontal mean is removed and then divided by the horizontal mean of the absolute values to have the amplitudes on the order of unity.
The magnitude of the noise can be adjusted by the input variable mag_noise_forc.
- Parameters:
- k_in_forcnumber
Inner wave number of the ring Defaults to 0.0.
- k_out_forcnumber
Outer wave number of the ring Defaults to 0.0.
- mag_noise_forcnumber
Amplitude of the forcing