Introduction
17 September 2022 by MiniUFO
[TOC]
1. Introduction
Researches on meteorology and oceanography usually encounter inversion problems that need to be solved numerically. One of the classical inversion problem is to solve for a streamfunction \(\psi\) given the vertical component of vorticity \(\zeta\) and proper boundary conditions:
Nowadays xarray becomes a popular data structure commonly used in Big Data Geoscience. Since the whole 4D data, as well as the coordinate information, are all combined into xarray, solving the inversion problem become quite straightforward and the only input would be just one xarray.DataArray of \(\zeta\). Inversion on the spherical earth, like some meteorological problems, could utilize the spherical harmonics like windspharm, which would be more efficient using FFT than SOR used here. However, in the case of ocean, SOR method is definitely a better choice in the presence of land/sea mask.
More importantly, this could be generalized into a numerical solver for second-order (or even fourth-order) partial differential equations using SOR method. Various popular inversion problems in geophysical fluid dynamics will be illustrated as examples.
One problem with SOR is that the speed of iteration using explicit python loops will be e-x-t-r-e-m-e-l-y … s-l-o-w! A very suitable solution here is to use numba. We may try our best to speed things up using more hardwares (possibly GPU).
2. Mathematical background
This package, xinvert, is designed to solve the following abstract equation:
where \(L\) is a second-order partial differential operator, \(\psi\) is the unknown to be inverted for, and \(F\) a prescribed forcing function. There could be also some parameters in the definition of \(L\), which should be specified before inverting \(\psi\).
For the 2D case, the general form of Eq. (1) is:
where coefficients \(A-G\) are all known variables. When the condition \(4AC-B^2>0\) is met everywhere in the domain, the above equation is an elliptic-type equation. In this case, one can invert \(\psi\) using the SOR iteration method. When \(4AC-B^2=0\) or \(4AC-B^2<0\), it is a parabolic or hyperbolic equation. In either case, SOR would fail to converge to the solution.
Sometimes the general form of Eq. (2) can be transformed into the standard form (i.e., standarization):
In this case, \(AD-BC>0\) should be met to insure its ellipticity. The elliptic condition has its own physical meaning in the problems of interest. That is, the system is in steady states that are stable to any small perturbation.
Many problems in meteorology and oceanography can be cast into the forms of either Eq. (2) or Eq. (3). However, some of them are formulated in 3D case (like the QG omega equation):
or in fourth-order case (Munk model):
So we implements four basic solvers to take into account the above four Eqs. (2-5) or cases. If a problem do not fit into one of these four types, we are going to add one solver for this type of problem. We hope NOT so because we want to keep the solvers as minimum and general as possible.
3. Physical problems at hands
Here we summarize all the inversion problems in meteorology and oceanography into the following table. The table can be extended further if one finds more problems that fit Eq. (1). Also, the top API calls for the users are also listed for references.
Names |
Equations |
A PI calls |
|---|---|---|
Horizontal streamfunction |
\[\nabla^2\psi=\frac{\partial^2 \psi}{\partial y^2}+\frac{\partial^2 \psi}{\partial x^2}=\zeta_k\]
|
|
MOC streamfunction |
\[\nabla^2\psi=\frac{\partial^2 \psi}{\partial z^2}+\frac{\partial^2 \psi}{\partial y^2}=\zeta_i\]
|
|
Walker streamfunction |
\[\nabla^2\psi=\frac{\partial^2 \psi}{\partial z^2}+\frac{\partial^2 \psi}{\partial x^2}=\zeta_j\]
|
|
Balanced mass field |
\[\nabla^2\Phi=\frac{\partial^2 \Phi}{\partial y^2}+\frac{\partial^2 \Phi}{\partial x^2}=F\]
|
|
Geostrophic balanced flow |
\[\frac{\partial}{\partial y}\left(f\frac{\partial \psi}{\partial y}\right)+\frac{\partial}{\partial x}\left(f\frac{\partial \psi}{\partial x}\right)=\nabla^2 \Phi\]
|
|
Eliassen model |
\[\frac{\partial}{\partial
p}\left(A\frac{\partial \psi}{\partial p}+B\frac{\partial \psi}{\partial y}\right)+\frac{\partial}{\partial y}\left(B\frac{\partial \psi}{\partial p}+C\frac{\partial \psi}{\partial y}\right)=F\]
|
|
Reference state for SWM |
\[\frac{\partial }{\partial y}\left(A\frac{\partial \Delta M}{\partial y}\right)+B\Delta M=F\]
|
|
PV inversion for balanced flow |
\[\frac{\pa\]
rtial}{partial theta}left(frac{2Lambda_0}{r^2}frac{partialLambda}{partial theta}right)+frac{partial}{partial r}left(frac{Gamma g}{Qr}frac{partialLambda}{partial r}right)=0 |
|
PV inversion for QG flow |
\[\frac{\partial}{\partial p}\left(\frac{f^2}{N^2}\frac{\partial \psi}{\partial p}\right)+\frac{\partial^2 \psi}{\partial y^2}=q\]
|
|
Gill-Matsuno model |
\[A\frac{\partial^2 \phi}{\partial y^2}+B\frac{\partial^2 \phi}{\partial x^2}+C\frac{\partial \phi}{\partial y}+D\frac{\partial \phi}{\partial x}+E\phi=Q\]
|
|
Stommel-Munk model |
\[A\nabla^4\psi-\frac{R}{D}\nabla^2\psi-\beta\frac{\partial \psi}{\partial x}=-\frac{\mathbf k\cdot\nabla\times\mathbf{\tau}}{\rho_0 D}\]
|
|
Fofonoff model |
\[\nabla^2\psi-c_0\psi=c_1-f\]
|
|
Bretherton-Haidvogel model |
\[\nabla^2\psi-\lambda D\psi=-\frac{f_0}{D}\eta_B\]
|
|
QG-Omega equation |
\[\frac{\partial}{\partial p}\left(f^2\frac{\partial \omega}{\partial p}\right)+\nabla\cdot\left(S\nabla\omega\right)=F\]
|
|
3D ocean flow |
\[\frac{\partial}{\partial p}\left(c_3\frac{\partial \phi}{\partial p}\right)+\nabla\cdot\left(c_1\nabla\phi-c_2\hat\nabla\phi\right)=F\]
|
|
… more, please add your own problems … |