Introducing PlesioGeostroPy: Quickstart

PlesioGeostroPy is a python realization of the plesio-gestrophy model. Plesio-Geostrophy model (see Jackson and Maffei, 2020), or PG model, is a 2-D reduced model for 3-D MHD in rapidly rotating frame, based on near-geostrophic flows in an axisymmetric cavity. The model uses a plesio-geostrophic ansatz, where the velocity field is completely described by a scalar stream function, but remains solenoidal and satisfies the non-penetration boundary condition. In diffusionless case, the MHD equations are shown to exactly reduce to a set of equations involving only the stream function and the axially integrated magnetic moments (except for the boundary terms), which are all 2-D variables.

It should be noted that the current development involves only a full sphere geometry. However, the method would apply to general axisymmetric cavities without an inner core.

Intro

Why PlesioGeostroPy?

Indeed, we already have Daria’s Wolfram Mathematica notebook, which also heavily inspired the initial development of this small package. The decision to develop a python realization is motivated by the following observations:

• Despite the completeness of Daria’s notebook, many notations are cryptic and not easy to apprehend.

  • Problem: The moments have been replaced with notations A to H, whereas the components of magnetic fields are replaced with notations a to c. For instance, field \(\overline{M_{ss}}=\overline{B_s^2}\) is A, \(\overline{M_{s\phi}} = \overline{B_sB_\phi}\) is C (or Csp for the perturbation field). The equations also use a confusing numbering system. Not only do we have eq1 to eq11, but also eq73, eq74, eq76 and eq80. These equations seem to correspond to the not-yet-linearized version of the induction equations. None of these are physically meaningful notations, and I even have jotted down a mapping table just to decipher the Mathematica code.

  • Cause: non-physical namings.

  • Response: Given the former two observations, a rewrite of the code is desired.

• The original mathematica notebook lacks flexibility when it comes to changing expansions, plottings, numerical computations, etc.

  • Problem: Changing expansions would require switching on and off code blocks, or creating new notebooks. Plotting (according to Stefano) has always been confusing in Mathematica. In addition to that, many of the Mathematica functions, e.g. special functions, etc. do not vectorize well, and would not be efficient in numerics.

  • Cause: many of the problems can be attributed to the fact that Mathematica is a domain-specific language (DSL). The inflexibility is due to the fact that Mathematica is mostly a functional language for symbolic manipulation, whereas the support of object-oriented programming is very limited (except for a OOP repo from 30 years ago, the links from which have mostly expired). The numerics part is not comparable to numerical libaries (e.g. numpy, scipy), let alone numerics DSL such as MATLAB. In addition, it is also complicated to use external libraries for faster numerical computation.

  • Response: Program the model in a language that is a general-purpose language (GPL), which allows abstract design of the code and at the same time allows efficient numerical algorithms. Designed initially to be a scripting language or glue language, and developing later into a GPL with numerous external packages, Python is a good choice. The packages such as sympy and numpy allow good integration between the symbolic and the numeric sides of the picture, while matplotlib and other visualization packages provide an additional selling point. It is also much easier to either directly call or export intermediate variables to external numerical libraries.

• Mathematica is closed source, while python itself as well as many of the packages are open source.

Of course, switching from Daria’s Mathematica notebook to PlesioGeostroPy implementation comes with a price.

• Symbolic computation may be more robust in Mathematica than in Python. This is statement recycled from some Mathematica users. This might be true for solving differential equations or other equations symbolically, however this is not relevant for PG model. In simple algebras, simplifications and substitutions, Mathematica is just as “capricious” as Sympy, the behaviour of neither would coincide exactly as what a human would do during derivations. On top of that, there is no knowing what exactly went wrong in Mathematica, whereas in Sympy, the user can peek under the hood as the code is open-source.

• Arbitrary precision computation is slower in vanilla-version SymPy, however slightly. This may be a real issue, but may also be circumvented by using multi-precision evaluation libraries, but it would be much more desirable if computing within the double precision is sufficient.

Prerequisites

The package PlesioGeostroPy is written in python, and dependent on the following packages.

# Core packages
python=3.8.13
numpy=1.22.3
scipy=1.8.0
sympy=1.11.1
h5py=3.3.0

# Visualization and presentation
pandas=1.4.2
matplotlib=3.5.2
jupyterlab=3.4.0

# Documentation and testing
sphinx=5.0.2
pytest=7.4.2
perfplot=0.10.2
line_profiler=4.1.1

In addition, in order to harness the power of the pretty printing interface, it is highly recommended that the interactive programming is done using web-based Project Jupyter tools, such as JupyterLab, Jupyter notebook, or IPython consoles launched in Jupyter services. These tools enable LaTeX rendering, which is a significantly improvement compared to Unicode or ASCII printing methods.

---

Getting started with SymPy

Sympy, or “Symbolic Python”, is a Python library for symbolic mathematics. A very brief demo is covered here. For a longer tutorial, see the official website.

As in using other Python libraries, to import utilities from the package one only needs to write

from sympy import *

Symbols, functions, expressions

As in most programming languages, one needs to define a variable before using it, and symbolic symbols are no exceptions in SymPy. This might be the only possible and reasonable thing to do for anyone with a background in C/C++, Java, C#, Fortran, Python, Julia or even MATLAB. The only scenarios, where I know the variables are readily interpreted as symbols even before definition are Mathematica and Maple.

To define a symbolic variable, which will have the name/expression \(\psi_1^{m}\), we can write the following code

psi_1m = Symbol(r"\hat{\psi}_1^m")
psi_1m

\[\displaystyle \hat{\psi}_1^m\]

This example already shows that you can use just about any LaTeX expressions as the symbol of the variable.

A symbol is just an independent variable, and has no known dependency on any other variables:

x, y, z = symbols("x, y, z")

diff_psi_1m_x = diff(psi_1m, x)
diff_psi_1m_x

\[\displaystyle 0\]

If a dependency on some variable is needed for some symbol, one needs to define a function:

psi = Function(r"\psi")(x, y, z)

diff_psi_x = diff(psi, (x, 2), (y, 1), (z, 0))
diff_psi_x

\[\displaystyle \frac{\partial^{3}}{\partial y\partial x^{2}} \psi{\left(x,y,z \right)}\]

Adding, substracting, multiplying or dividing any symbols or functions will generate an expression (Expr):

psi_expr = psi_1m**2*x + psi_1m*y + z
psi_expr

\[\displaystyle \left(\hat{\psi}_1^m\right)^{2} x + \hat{\psi}_1^m y + z\]

An expression can once again be added, substracted, multiplied, and differentiated etc. as well. To construct a derivative object without evaluation, use Derivative:

dpsi_expr_dx = Derivative(psi_expr, x)
dpsi_expr_dx

\[\displaystyle \frac{\partial}{\partial x} \left(\left(\hat{\psi}_1^m\right)^{2} x + \hat{\psi}_1^m y + z\right)\]

One can evaluate this by either invoking the method .doit(), or directly using diff method:

display(dpsi_expr_dx.doit())
display(diff(psi_expr, x))

\[\displaystyle \left(\hat{\psi}_1^m\right)^{2}\]

\[\displaystyle \left(\hat{\psi}_1^m\right)^{2}\]

In sympy, there is a clear hierarchy of classes. Expr is pretty much one of the most basic base classes. This means a Function is an Expr, a Symbol is an Expr, and a jacobi polynomial is also an Expr. The results of algebraic manipulations, which may be Add, Mul, are all Exprs. For readers who do not know object-oriented programming, by A is Expr, what I really mean is that A is a child class of, or inherit from, Expr.

Substitutions, simplifications

Substitution in Sympy is realized by the Expr.subs() method. A typical way to use the substitution method is to pass a dictionary (a key-value “map” in C++ languages) into the method.

This is much more human readable than the /. operator in Mathematica.

The following example shows how the \(\psi\) function is replaced with a Jacobi polynomial of order \(n\).

alpha = Symbol(r"\alpha")
beta = Symbol(r"\beta")
n = Symbol('n', integer=True)

diff_jacobi = diff_psi_x.subs({psi: jacobi(n, alpha, beta, x)})
diff_jacobi

\[\displaystyle \frac{\partial^{3}}{\partial y\partial x^{2}} P_{n}^{\left(\alpha,\beta\right)}\left(x\right)\]

Substituting a variable whose derivative is involved would result in delayed substitution:

diff_val = diff_psi_x.subs({x: 1, y: 1})
diff_val

\[\begin{split}\displaystyle \left. \frac{\partial^{3}}{\partial y\partial x^{2}} \psi{\left(x,y,z \right)} \right|_{\substack{ y=1\\ x=1 }}\end{split}\]

which can only be evaluated when the explicit form is calculated:

diff_val = diff_val.subs({psi: x**3*y**2}).doit()
diff_val

\[\displaystyle 12\]

Simplification is probably the most counterintuitive part of all symbolic computation libraries, simply due to the fact that the computer does not know for sure what the human programmer wants.

In many cases, when we are doing simplifications, we also do it by trial and error, and we don’t even know a priori what form we want to end up in.

We can make a dummy example here,

diff_simp = diff_psi_x.subs({psi: x**3*y**2 + 1/x/y}).simplify()
diff_simp

\[\displaystyle 2 \cdot \left(7 x y - \frac{x^{4} y^{3} + 1}{x^{3} y^{2}}\right)\]

Apparently, it is much reasonable to rewrite it as $\( 2\left(6xy - \frac{1}{x^3 y^2}\right) \)$, this can still be achieved using expand()

diff_simp.expand()

\[\displaystyle 12 x y - \frac{2}{x^{3} y^{2}}\]

For more on simplifying expressions such as rational expressions, special functions, etc., see Simplifiations in SymPy.

---

Getting started with PlesioGeostroPy

Now we have come to the fun part: the actual utilities concerning the PG model.

The PG equations, especially the induction equations, are of such a complex form that a sane human being is usually not overly fond of writing the expression down very often. PlesioGeostroPy does just that for you, by wrapping the governing equations in much human-readable structures.

In this section, we shall look at the symbolic part of this package. We will see how to inspect equations, transformations and their linearized versions.

Core variables

The core variables that are concerned in the model are constructed in pg_model.core module. These first include the coordinates \(s\), \(\phi\), \(z\); Cartesian coordinates \(x\), \(y\); as well as spherical coordinates \(r\), \(\theta\).

import os, sys
root_dir = "."

# The following 2 lines are a hack to import from parent directory
# If the notebook is run in the root directory, comment out these 2 lines
sys.path.append(os.path.dirname(os.getcwd()))
root_dir = ".."

from sympy import *
from pg_utils.pg_model import core
from pg_utils.pg_model.core import s, p, z, t

print("Cylindrical oordinates:")
display(*[core.cyl[i_comp] for i_comp in range(3)])

Cylindrical oordinates:

\[\displaystyle s\]

\[\displaystyle \phi\]

\[\displaystyle z\]

Next, we have the 3-D vector fields. These include

• magnetic field: in cylindrical coordinates: B_vec, background field B0_vec; in spherical coordinates B_sph, background field B0_sph;

• velocity field: in cylindrical coordinates: U_vec, background field U0_vec; in spherical coordinates U_sph, background field U0_sph;

we can take a look at one of these variables.

display(*[core.B_vec[i_comp] for i_comp in range(3)])

\[\displaystyle B_{s}{\left(s,\phi,z,t \right)}\]

\[\displaystyle B_{\phi}{\left(s,\phi,z,t \right)}\]

\[\displaystyle B_{z}{\left(s,\phi,z,t \right)}\]

And finally, the PG variables. These are stored in core.pgvar, core.pgvar_bg (background) and core.pgvar_ptb (perturbation). These are of the type

type(core.pgvar)

pg_utils.pg_model.base.CollectionPG

CollectionPG is a special class specifically designed for holding PG variables. This class allows accessing fields as attributes or by index, as well as many other useful features. For more detailed explanation of how this interface works, please refer to Demo_Collections.

To look at what variables the object holds, we can for instance make the following inquiry

core.pgvar._field_names

['Psi',
 'Mss',
 'Mpp',
 'Msp',
 'Msz',
 'Mpz',
 'zMss',
 'zMpp',
 'zMsp',
 'Bs_e',
 'Bp_e',
 'Bz_e',
 'dBs_dz_e',
 'dBp_dz_e',
 'Br_b',
 'Bs_p',
 'Bp_p',
 'Bz_p',
 'Bs_m',
 'Bp_m',
 'Bz_m']

We see that we have in all 21 variables. The first is the stream function; the second to the eight are magnetic moments; then we have the fields in the equatorial plane (ones with _e suffix). Finally, we have the boundary terms.

Note that in standard PG, only Br_b = “radial magnetic field at the boundary” is used; however, for many problems such as the eigenvalue problem, we can use Bs_p (s-magnetic field at \(z=+H\)) - Bz_m (z-magnetic field at \(z=-H\)) instead of Br_b. As mentioned, all these fields can be accessed just by typing in the attribute name:

display(core.pgvar.Psi, core.pgvar.Bz_m)

\[\displaystyle \Psi{\left(s,\phi,t \right)}\]

\[\displaystyle B^{-}_{z}{\left(s,\phi,t \right)}\]

In addition to this method of accessing fields, the same variable pgvar can also be indexed.

display(core.pgvar[0], core.pgvar[-1])

\[\displaystyle \Psi{\left(s,\phi,t \right)}\]

\[\displaystyle B^{-}_{z}{\left(s,\phi,t \right)}\]

And finally, as a sugar frosting, the fields can be traversed by iterating through the variable pgvar.

display(*[field_symb for field_symb in core.pgvar])

\[\displaystyle \Psi{\left(s,\phi,t \right)}\]

\[\displaystyle \overline{M_{ss}}{\left(s,\phi,t \right)}\]

\[\displaystyle \overline{M_{\phi\phi}}{\left(s,\phi,t \right)}\]

\[\displaystyle \overline{M_{s\phi}}{\left(s,\phi,t \right)}\]

\[\displaystyle \widetilde{M_{sz}}{\left(s,\phi,t \right)}\]

\[\displaystyle \widetilde{M_{\phi z}}{\left(s,\phi,t \right)}\]

\[\displaystyle \widetilde{zM_{ss}}{\left(s,\phi,t \right)}\]

\[\displaystyle \widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)}\]

\[\displaystyle \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)}\]

\[\displaystyle B_{s}^e{\left(s,\phi,t \right)}\]

\[\displaystyle B_{\phi}^e{\left(s,\phi,t \right)}\]

\[\displaystyle B_{z}^e{\left(s,\phi,t \right)}\]

\[\displaystyle B_{s, z}^e{\left(s,\phi,t \right)}\]

\[\displaystyle B_{\phi, z}^e{\left(s,\phi,t \right)}\]

\[\displaystyle B_{r1}{\left(\theta,\phi,t \right)}\]

\[\displaystyle B^{+}_{s}{\left(s,\phi,t \right)}\]

\[\displaystyle B^{+}_{\phi}{\left(s,\phi,t \right)}\]

\[\displaystyle B^{+}_{z}{\left(s,\phi,t \right)}\]

\[\displaystyle B^{-}_{s}{\left(s,\phi,t \right)}\]

\[\displaystyle B^{-}_{\phi}{\left(s,\phi,t \right)}\]

\[\displaystyle B^{-}_{z}{\left(s,\phi,t \right)}\]

The same thing holds for pgvar_bg and pgvar_ptb. The background fields are denoted with an additional \(0\) superscript, and the perturbed fields are denoted with lower-case letters.

display(core.pgvar_bg.Mss, core.pgvar_ptb.Mss)

\[\displaystyle \overline{M_{ss}}^0{\left(s,\phi \right)}\]

\[\displaystyle \overline{m_{ss}}{\left(s,\phi,t \right)}\]

Last but not least, the core module also stores the PG ansatz of the velocity, which can be invoked whenever necessary:

display(*[core.U_pg[i_comp] for i_comp in range(3)])

\[\displaystyle \frac{\frac{\partial}{\partial \phi} \Psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}}\]

\[\displaystyle - \frac{\frac{\partial}{\partial s} \Psi{\left(s,\phi,t \right)}}{H{\left(s \right)}}\]

\[\displaystyle \frac{z \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \Psi{\left(s,\phi,t \right)}}{s H^{2}{\left(s \right)}}\]

Equations and linearizations

With all the variables defined, we can now introduce the equations. The 15 PG equations have been typed in manually in the equations.py module. Since these equations are constructed when the module is imported, the first import of equations module takes some time (~10s or so).

These are similarly stored in a CollectionPG object called equations.eqs_pg. Let us inspect these equations:

from pg_utils.pg_model import equations as eqs

display(*[eq for eq in eqs.eqs_pg])

\[\displaystyle \left(- \frac{\frac{d}{d s} H{\left(s \right)}}{2 H^{2}{\left(s \right)}} + \frac{1}{s H{\left(s \right)}}\right) \frac{\partial^{3}}{\partial t\partial \phi^{2}} \Psi{\left(s,\phi,t \right)} + \frac{\partial}{\partial s} \frac{s \frac{\partial^{2}}{\partial t\partial s} \Psi{\left(s,\phi,t \right)}}{H{\left(s \right)}} = - \frac{s \left(\frac{s \frac{\partial}{\partial s} \overline{f_\phi}{\left(s,\phi,t \right)} + \overline{f_\phi}{\left(s,\phi,t \right)}}{s} - \frac{\frac{\partial}{\partial \phi} \overline{f_s}{\left(s,\phi,t \right)}}{s}\right)}{2 H{\left(s \right)}} + \left(\frac{s f_{\phi}^e{\left(s,\phi,t \right)}}{H{\left(s \right)}} + \frac{\frac{\partial}{\partial \phi} \widetilde{f_z}{\left(s,\phi,t \right)}}{2 H{\left(s \right)}}\right) \frac{d}{d s} H{\left(s \right)} - \frac{2 \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \Psi{\left(s,\phi,t \right)}}{H^{2}{\left(s \right)}}\]

\[\displaystyle \frac{\partial}{\partial t} \overline{M_{ss}}{\left(s,\phi,t \right)} = - \left(U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \frac{\overline{M_{ss}}{\left(s,\phi,t \right)}}{H{\left(s \right)}} + \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \frac{\overline{M_{ss}}{\left(s,\phi,t \right)}}{H{\left(s \right)}}}{s}\right) H{\left(s \right)} + 2 \overline{M_{ss}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + \frac{2 \overline{M_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} \overline{M_{\phi\phi}}{\left(s,\phi,t \right)} = 2 s \overline{M_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{\phi}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} H{\left(s \right)} \overline{M_{\phi\phi}}{\left(s,\phi,t \right)} + \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} H{\left(s \right)} \overline{M_{\phi\phi}}{\left(s,\phi,t \right)}}{s}}{H{\left(s \right)}} - 2 \overline{M_{\phi\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)}\]

\[\displaystyle \frac{\partial}{\partial t} \overline{M_{s\phi}}{\left(s,\phi,t \right)} = s \overline{M_{ss}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{\phi}{\left(s,\phi,z,t \right)}}{s} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \overline{M_{s\phi}}{\left(s,\phi,t \right)} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \overline{M_{s\phi}}{\left(s,\phi,t \right)}}{s} + \frac{\overline{M_{\phi\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} \widetilde{M_{sz}}{\left(s,\phi,t \right)} = \left(\frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + 2 \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)}\right) \widetilde{M_{sz}}{\left(s,\phi,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \widetilde{M_{sz}}{\left(s,\phi,t \right)} + \widetilde{zM_{ss}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{s}{\left(s,\phi,z,t \right)} \frac{d}{d s} H{\left(s \right)}}{H{\left(s \right)}} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \widetilde{M_{sz}}{\left(s,\phi,t \right)}}{s} + \frac{\widetilde{M_{\phi z}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s} + \frac{\widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s H{\left(s \right)}}\]

\[\displaystyle \frac{\partial}{\partial t} \widetilde{M_{\phi z}}{\left(s,\phi,t \right)} = s \widetilde{M_{sz}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{\phi}{\left(s,\phi,z,t \right)}}{s} + \left(- \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)}\right) \widetilde{M_{\phi z}}{\left(s,\phi,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \widetilde{M_{\phi z}}{\left(s,\phi,t \right)} + \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{s}{\left(s,\phi,z,t \right)} \frac{d}{d s} H{\left(s \right)}}{H{\left(s \right)}} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \widetilde{M_{\phi z}}{\left(s,\phi,t \right)}}{s} + \frac{\widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s H{\left(s \right)}}\]

\[\displaystyle \frac{\partial}{\partial t} \widetilde{zM_{ss}}{\left(s,\phi,t \right)} = \left(2 \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + 2 \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)}\right) \widetilde{zM_{ss}}{\left(s,\phi,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \widetilde{zM_{ss}}{\left(s,\phi,t \right)} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \widetilde{zM_{ss}}{\left(s,\phi,t \right)}}{s} + \frac{2 \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} \widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)} = 2 s \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{\phi}{\left(s,\phi,z,t \right)}}{s} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)} - 2 \widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} = s \widetilde{zM_{ss}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \frac{U_{\phi}{\left(s,\phi,z,t \right)}}{s} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} + \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} \widetilde{zM_{s\phi}}{\left(s,\phi,t \right)}}{s} + \frac{\widetilde{zM_{\phi\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{s}^e{\left(s,\phi,t \right)} = B_{s}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{s}^e{\left(s,\phi,t \right)} + \frac{B_{\phi}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{s}^e{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{\phi}^e{\left(s,\phi,t \right)} = B_{s}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{\phi}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{\phi}^e{\left(s,\phi,t \right)} + \frac{B_{\phi}^e{\left(s,\phi,t \right)} U_{s}{\left(s,\phi,z,t \right)} - B_{s}^e{\left(s,\phi,t \right)} U_{\phi}{\left(s,\phi,z,t \right)}}{s} + \frac{B_{\phi}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{\phi}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{\phi}^e{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{z}^e{\left(s,\phi,t \right)} = B_{z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{z}^e{\left(s,\phi,t \right)} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{z}^e{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{s, z}^e{\left(s,\phi,t \right)} = B_{s, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} - B_{s, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{s, z}^e{\left(s,\phi,t \right)} + \frac{B_{\phi, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{s, z}^e{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{\phi, z}^e{\left(s,\phi,t \right)} = - B_{\phi, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} + B_{s, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{\phi}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{\phi, z}^e{\left(s,\phi,t \right)} + \frac{B_{\phi, z}^e{\left(s,\phi,t \right)} U_{s}{\left(s,\phi,z,t \right)} - B_{s, z}^e{\left(s,\phi,t \right)} U_{\phi}{\left(s,\phi,z,t \right)}}{s} + \frac{B_{\phi, z}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{\phi}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{\phi, z}^e{\left(s,\phi,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B_{r1}{\left(\theta,\phi,t \right)} = - \frac{\frac{\partial}{\partial \phi} B_{r1}{\left(\theta,\phi,t \right)} U_{\phi}{\left(r,\theta,\phi,t \right)}}{r \sin{\left(\theta \right)}} - \frac{\frac{\partial}{\partial \theta} B_{r1}{\left(\theta,\phi,t \right)} U_{\theta}{\left(r,\theta,\phi,t \right)} \sin{\left(\theta \right)}}{r \sin{\left(\theta \right)}}\]

\[\displaystyle \frac{\partial}{\partial t} B^{+}_{s}{\left(s,\phi,t \right)} = B^{+}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + B^{+}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{s}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{s}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{s}{\left(s,\phi,z,t \right)} + \frac{B^{+}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B^{+}_{\phi}{\left(s,\phi,t \right)} = B^{+}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{\phi}{\left(s,\phi,z,t \right)} + B^{+}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{\phi}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{\phi}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{\phi}{\left(s,\phi,z,t \right)} + \frac{B^{+}_{\phi}{\left(s,\phi,t \right)} U_{s}{\left(s,\phi,z,t \right)} - B^{+}_{s}{\left(s,\phi,t \right)} U_{\phi}{\left(s,\phi,z,t \right)}}{s} + \frac{B^{+}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{\phi}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{\phi}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B^{+}_{z}{\left(s,\phi,t \right)} = B^{+}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{z}{\left(s,\phi,z,t \right)} + B^{+}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{z}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{z}{\left(s,\phi,z,t \right)} + \frac{B^{+}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{z}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{z}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B^{-}_{s}{\left(s,\phi,t \right)} = B^{-}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{s}{\left(s,\phi,z,t \right)} + B^{-}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{s}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{s}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{s}{\left(s,\phi,z,t \right)} + \frac{B^{-}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{s}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{s}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B^{-}_{\phi}{\left(s,\phi,t \right)} = B^{-}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{\phi}{\left(s,\phi,z,t \right)} + B^{-}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{\phi}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{\phi}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{\phi}{\left(s,\phi,z,t \right)} + \frac{B^{-}_{\phi}{\left(s,\phi,t \right)} U_{s}{\left(s,\phi,z,t \right)} - B^{-}_{s}{\left(s,\phi,t \right)} U_{\phi}{\left(s,\phi,z,t \right)}}{s} + \frac{B^{-}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{\phi}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{\phi}{\left(s,\phi,z,t \right)}}{s}\]

\[\displaystyle \frac{\partial}{\partial t} B^{-}_{z}{\left(s,\phi,t \right)} = B^{-}_{s}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U_{z}{\left(s,\phi,z,t \right)} + B^{-}_{z}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U_{z}{\left(s,\phi,z,t \right)} - U_{s}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial s} B_{z}{\left(s,\phi,z,t \right)} - U_{z}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial z} B_{z}{\left(s,\phi,z,t \right)} + \frac{B^{-}_{\phi}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U_{z}{\left(s,\phi,z,t \right)}}{s} - \frac{U_{\phi}{\left(s,\phi,z,t \right)} \frac{\partial}{\partial \phi} B_{z}{\left(s,\phi,z,t \right)}}{s}\]

Of course, as other CollectionPG objects as core.pgvar, the equations in equations.eqs_pg can also be accessed either by attribute (e.g. equations.eqs_pg.Mss) or by index (e.g. equations.eqs_pg[1]).

These equations are then linearized automatically in the same module. The procedure is as follows: each fields are expanded as the background state + \(\epsilon\cdot\) a perturbation state. Then, all the terms that are linear in \(\epsilon\) are collected to yields the linearized form. This method is borrowed from Daria’s Mathematica notebook.

The linearized equations are then stored in equations.eqs_lin. For instance, the linearized evolution equation for \(\widetilde{zM_{s\phi}}\) takes the form

eqs.eqs_pg_lin.zMsp

\[\displaystyle \frac{\partial}{\partial t} \widetilde{zm_{s\phi}}{\left(s,\phi,t \right)} = - U^{0}_{s}{\left(s,\phi,z \right)} \frac{\partial}{\partial s} \widetilde{zm_{s\phi}}{\left(s,\phi,t \right)} + \widetilde{zm_{s\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial z} U^{0}_{z}{\left(s,\phi,z \right)} + \widetilde{zm_{ss}}{\left(s,\phi,t \right)} \frac{\partial}{\partial s} U^{0}_{\phi}{\left(s,\phi,z \right)} - \frac{\widetilde{zM_{ss}}^0{\left(s,\phi \right)} \frac{\partial^{2}}{\partial s^{2}} \psi{\left(s,\phi,t \right)}}{H{\left(s \right)}} + \frac{\widetilde{zM_{ss}}^0{\left(s,\phi \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial s} \psi{\left(s,\phi,t \right)}}{H^{2}{\left(s \right)}} - \frac{U^{0}_{\phi}{\left(s,\phi,z \right)} \widetilde{zm_{ss}}{\left(s,\phi,t \right)}}{s} - \frac{U^{0}_{\phi}{\left(s,\phi,z \right)} \frac{\partial}{\partial \phi} \widetilde{zm_{s\phi}}{\left(s,\phi,t \right)}}{s} + \frac{\widetilde{zm_{\phi\phi}}{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} U^{0}_{s}{\left(s,\phi,z \right)}}{s} + \frac{\widetilde{zM_{ss}}^0{\left(s,\phi \right)} \frac{\partial}{\partial s} \psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}} - \frac{\frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)} \frac{\partial}{\partial s} \widetilde{zM_{s\phi}}^0{\left(s,\phi \right)}}{s H{\left(s \right)}} + \frac{\frac{\partial}{\partial s} \psi{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} \widetilde{zM_{s\phi}}^0{\left(s,\phi \right)}}{s H{\left(s \right)}} + \frac{\widetilde{zM_{s\phi}}^0{\left(s,\phi \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{s H^{2}{\left(s \right)}} + \frac{\widetilde{zM_{\phi\phi}}^0{\left(s,\phi \right)} \frac{\partial^{2}}{\partial \phi^{2}} \psi{\left(s,\phi,t \right)}}{s^{2} H{\left(s \right)}}\]

• The vorticity equation has been validated against Daria’s notebook quad_malkus_reg_diff.nb, and has passed the test.

• Induction terms of \(\overline{m_{ss}}\), \(\overline{m_{\phi\phi}}\), \(\overline{m_{s\phi}}\), \(\widetilde{m_{sz}}\), \(\widetilde{m_{\phi z}}\), \(\widetilde{zm_{ss}}\), \(z\widetilde{m_{\phi\phi}}\), \(z\widetilde{m_{s\phi}}\), as well as \(b_{es}\) and \(b_{e\phi}\) (i.e. the linearized induction term), have been validated against quad_malkus_reg_diff.nb and have passed the test.

eqs.eqs_pg_lin.Bs_p.subs({core.U0_vec.s: 0, core.U0_vec.p: 0, core.U0_vec.z: 0}).doit()

\[\displaystyle \frac{\partial}{\partial t} b^{+}_{s}{\left(s,\phi,t \right)} = - \frac{z \frac{\partial}{\partial z} B^{0}_{s}{\left(s,\phi,z \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{s H^{2}{\left(s \right)}} + \frac{B_s^{0+}{\left(s,\phi \right)} \frac{\partial^{2}}{\partial s\partial \phi} \psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}} - \frac{B_s^{0+}{\left(s,\phi \right)} \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{s H^{2}{\left(s \right)}} + \frac{\frac{\partial}{\partial \phi} B^{0}_{s}{\left(s,\phi,z \right)} \frac{\partial}{\partial s} \psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}} - \frac{\frac{\partial}{\partial s} B^{0}_{s}{\left(s,\phi,z \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}} + \frac{B_\phi^{0+}{\left(s,\phi \right)} \frac{\partial^{2}}{\partial \phi^{2}} \psi{\left(s,\phi,t \right)}}{s^{2} H{\left(s \right)}} - \frac{B_s^{0+}{\left(s,\phi \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{s^{2} H{\left(s \right)}}\]

• Induction terms of \(b_s^\pm\), \(b_\phi^\pm\), \(b_z^\pm\) has been validated against the expressions in the document PG_Assim.pdf, which are derived by hand. The validation is successful.

• The validation is performed under zero background velocity.

We note that in the vorticity equation, the placeholders for the forces are used:

eqs.eqs_pg_lin.Psi

\[\displaystyle \frac{s \frac{\partial^{3}}{\partial t\partial s^{2}} \psi{\left(s,\phi,t \right)}}{H{\left(s \right)}} - \frac{s \frac{d}{d s} H{\left(s \right)} \frac{\partial^{2}}{\partial t\partial s} \psi{\left(s,\phi,t \right)}}{H^{2}{\left(s \right)}} + \frac{\frac{\partial^{2}}{\partial t\partial s} \psi{\left(s,\phi,t \right)}}{H{\left(s \right)}} - \frac{\frac{d}{d s} H{\left(s \right)} \frac{\partial^{3}}{\partial t\partial \phi^{2}} \psi{\left(s,\phi,t \right)}}{2 H^{2}{\left(s \right)}} + \frac{\frac{\partial^{3}}{\partial t\partial \phi^{2}} \psi{\left(s,\phi,t \right)}}{s H{\left(s \right)}} = \frac{s f_{\phi}^e{\left(s,\phi,t \right)} \frac{d}{d s} H{\left(s \right)}}{H{\left(s \right)}} - \frac{s \frac{\partial}{\partial s} \overline{f_\phi}{\left(s,\phi,t \right)}}{2 H{\left(s \right)}} - \frac{\overline{f_\phi}{\left(s,\phi,t \right)}}{2 H{\left(s \right)}} + \frac{\frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \widetilde{f_z}{\left(s,\phi,t \right)}}{2 H{\left(s \right)}} + \frac{\frac{\partial}{\partial \phi} \overline{f_s}{\left(s,\phi,t \right)}}{2 H{\left(s \right)}} - \frac{2 \frac{d}{d s} H{\left(s \right)} \frac{\partial}{\partial \phi} \psi{\left(s,\phi,t \right)}}{H^{2}{\left(s \right)}}\]

where \(\overline{f_s}\), \(\overline{f_\phi}\), \(\widetilde{f_z}\) and \(f_{e\phi}\) are the vertically evenly integrated \(s\)-component, \(\phi\)-component, vertically oddly integrated \(z\)-component and the equatorial \(\phi\)-component of arbitrary body force \(\mathbf{f}\), respectively.

This allows some flexibility regarding which forces are included in the system. Currently, only the forms of the Lorentz force are compiled in the package, but adding viscous diffusion as well as buoyancy force will be straightforward.

Precomputed equations

As you probably have noticed, loading the equations from the pg_utils.pg_model.equations module is slooowwww!

Loading the equations now may take 20-30s (tested on my laptop, AMD R5-4500U).

The reason is that it takes a while to compile the equations, as they are not all written out in explicit forms in the code. In linearization especially, the program needs to simplify the expressions to a form where the perturbation \(\epsilon\) is taken out of all brackets. Another problem is that so far, compilation and derivation of equations is done sequentially, not in parallel.

While the compilation of equations is not expected to be the bottleneck of the program in the end, it would still be desirable if it can be loaded faster. A method is to compute the equations and store them in a string that can be quickly loaded. This has already been done and the outputs have been saved under out/symbolic/. One can use the following code to load a set of equations:

from pg_utils.pg_model import base

with open(os.path.join(root_dir, "out/symbolic/eqs_pg.json"), 'r') as fread:
    pgeqs = base.CollectionPG.load_json(fread, parser=parse_expr)

which cuts back the cost of loading from 20s to 2s. The following four sets are given in the folder:

• eqs_pg.json: original PG equations

• eqs_pg_lin.json: linearized PG equations

• eqs_cg.json: conjugate variable equations

• eqs_cg_lin.json: linearized conjugate variable equations

This is now the recommended way to load equations. However, it also means that if the equations in pg_utils.pg_model.equations are changed, they would need to be re-compiled and re-saved.

Forcing and linearization

As mentioned, the forcing in the vorticity equation is only given by placeholder functions. It remains to be defined what forces are involved in the system. This is the task for the forcing.py module.

Currently, only the forms of the Lorentz force are compiled in the package. The explicit forms for Lorentz force are stored in forcing.Ls_sym_expr, forcing.Lp_sym_expr, forcing.Lz_asym_expr and forcing.Le_p_expr.

from pg_utils.pg_model import forcing

Eq(forcing.Le_p, forcing.Le_p_expr)

\[\displaystyle L_{\phi}^e{\left(s,\phi,t \right)} = B_{\phi, z}^e{\left(s,\phi,t \right)} B_{z}^e{\left(s,\phi,t \right)} + B_{s}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial s} B_{\phi}^e{\left(s,\phi,t \right)} + \frac{B_{\phi}^e{\left(s,\phi,t \right)} B_{s}^e{\left(s,\phi,t \right)}}{s} + \frac{B_{\phi}^e{\left(s,\phi,t \right)} \frac{\partial}{\partial \phi} B_{\phi}^e{\left(s,\phi,t \right)}}{s}\]

The linearization of the forces are the same as linearization of equations and is done automatically. The linearized forcing terms are stored in forcing.Ls_sym_lin, forcing.Lp_sym_lin, forcing.Lz_asym_lin and forcing.Le_p_lin expressions.

Eq(forcing.Ls_sym, forcing.Ls_sym_lin)

\[\displaystyle \overline{L_s}{\left(s,\phi,t \right)} = \frac{2 s B_s^{0+}{\left(s,\phi \right)} b^{+}_{s}{\left(s,\phi,t \right)}}{H{\left(s \right)}} + \frac{2 s B_s^{0-}{\left(s,\phi \right)} b^{-}_{s}{\left(s,\phi,t \right)}}{H{\left(s \right)}} + B_s^{0+}{\left(s,\phi \right)} b^{+}_{z}{\left(s,\phi,t \right)} - B_s^{0-}{\left(s,\phi \right)} b^{-}_{z}{\left(s,\phi,t \right)} + B_z^{0+}{\left(s,\phi \right)} b^{+}_{s}{\left(s,\phi,t \right)} - B_z^{0-}{\left(s,\phi \right)} b^{-}_{s}{\left(s,\phi,t \right)} + \frac{\partial}{\partial s} \overline{m_{ss}}{\left(s,\phi,t \right)} - \frac{\overline{m_{\phi\phi}}{\left(s,\phi,t \right)}}{s} + \frac{\overline{m_{ss}}{\left(s,\phi,t \right)}}{s} + \frac{\frac{\partial}{\partial \phi} \overline{m_{s\phi}}{\left(s,\phi,t \right)}}{s}\]

Linearized expressions for Lorentz force are validated against Daria’s notebook quad_malkus_reg_diff.nb.

• \(\overline{L_\phi}\), \(\widetilde{L_z}\) and \(L_\phi(z=0)\) passed the validation;

• The \(s\)-component of the equatorial Lorentz force in quad_malkus_reg_diff has a \(\frac{s^2}{H}\) coefficient, while the current derivation (and others) indicate this coefficient should be \(\frac{s}{H}\)

For a more detailed explanation on the ingredients, see Variables, equations and forcing.

To be continued

Here are more demos to learn about PlesioGeostroPy:

• The collection interface

• The variables, equations and forcing

• Regularity conditions of tensor components in polar coordinates