Draft:MPDATA

In numerical analysis of partial differential equations, MPDATA is a family of iterative finite-difference/finite-volume methods for numerically integrating of hyperbolic differential equations modelling conservation laws of the form:

$$\partial _{t}(G\psi )+\nabla \cdot ({\mathbf {v}}\psi )=GR,$$ 1

where $\psi ({\mathbf {x}},t)$ is an advected scalar field (or advectee); ${\mathbf {v}}=G{\dot {\mathbf {x}}}$ is the advective velocity vector field (or advector), $G({\mathbf {x}})$ may either play the role of the fluid density, the Jacobian of coordinate transformation from Cartesian to curvilinear framework ${\mathbf {x}}$, or their product; $R({\mathbf {x}})$ combines all source terms.

MPDATA stands for Multidimensional Positive Definite Advection Transport Algorithm. The algorithm was formulated by Piotr K. Smolarkiewicz (and is also referred to as Smolarkiewicz's method) at the US National Center for Atmospheric Research - NCAR (at the time, Smolarkiewicz was a fellow of the NCAR Advanced Study Program and a recent graduate from Krzysztof Haman's group at the University of Warsaw; the seminal 1983 and 1984 MPDATA publications mention both institutions).

The crux of the method lies in iterative application of the upwind scheme. The first iteration employs the advective velocity ${\mathbf {v}}$, while each subsequent iteration employs a so-called antidiffusive velocity which corrects solution from prior iteration reducing the numerical diffusion. The antidiffusive velocities are formulated through modified equation analysis of the upwind scheme and feature cross-dimensional dependencies (i.e., applying MPDATA in multiple dimensions is not equivalent to application of one-dimensional MPDATA in all dimensions), the scheme is thus not dimensionally split, hence "M" in the algorithm name. Since each iteration of the scheme constitutes a forward-in-time upwind pass, the scheme inherits characteristics of the upwind scheme: CFL stability criterion, conservativeness, embarrassingly parallel domain decomposition, and sign-preservation. For non-negative fields $\psi$, sign-preservation translates to positive definiteness, hence the "PD" in the algorithm name. Application of the corrective iterations improves scheme convergence rate compared with first-order upwind. Depending on the MPDATA variant, convergence of up to third-order in time and space can be achieved.

While the original formulation of MPDATA employed structured grids, the algorithm has been subsequently also formulated for unstructured grids.

Description of the basic scheme in 1D

MPDATA is inherently multi-dimensional, and primarily used in computational fluid dynamics where the advective volocities and problem geometries are variable in time. Still, the key idea underlying the MPDATA approach can be conveyed with a basic example of solenoidal stationary flow in one dimension (i.e., ${\mathbf {v}}=[u]$ constant in time and space), without coordinate transformation ($G=1$), for the case of homogeneous advection ($R=0$) of a nonnegative scalar field ($\psi \geq 0$), with the following flux form of the advection equation:

$$\partial _{t}\psi +\partial _{x}(u\psi )=0.$$ 2

Upwind discretisation of the problem on a regular staggered grid with a time step $\Delta t$ and a grid step $\Delta x$, with $n=t/\Delta t$, $i=x/\Delta x$, and the half-integer spatial indices corresponding to grid-cell walls:

   $u_{i-3/2}$     $u_{i-1/2}$      $u_{i+1/2}$     $u_{i+3/2}$
···  |    ·    |    ·    |    ·    |  ···
         $\psi _{i-1}$       $\psi _{i}$       $\psi _{i+1}$

can be formulated with (for arbitrary $v$):

$${\frac {\psi _{i}^{n+1}-\psi _{i}^{n}}{\Delta t}}+{\frac {\overbrace {f(\psi _{i}^{n},\psi _{i+1}^{n},v_{i+1/2}^{n})} ^{\text{left-hand wall flux}}-\overbrace {f(\psi _{i-1}^{n},\psi _{i}^{n},v_{i-1/2}^{n})} ^{\text{right-hand wall flux}}}{\Delta x}}=0$$ 3

with the flux function defined using positive and negative parts of $v_{i\pm 1/2}$ as:

$$f(\psi _{l},\psi _{r},u)=\!\!\overbrace {\frac {u+|u|}{2}} ^{\text{positive part}}\!\psi _{l}+\!\!\overbrace {\frac {u-|u|}{2}} ^{\text{negative part}}\!\psi _{r}.$$ 4

Introducing the non-dimensional Courant number $C=v{\Delta t}/{\Delta x}$, the resultant explicit-in-time scheme (referred to as "upwind", "upstream" or "donor-cell"), for a constant $C$ reads:

$$\psi _{i}^{n+1}=\psi _{i}^{n}-\left[{\frac {C+|C|}{2}}(\psi _{i}^{n}-\psi _{i-1}^{n})+{\frac {C-|C|}{2}}(\psi _{i+1}^{n}-\psi _{i}^{n})\right]$$ 5

which is conservative, sing-preserving and stable for $|C|\leq 1$. However, the scheme is only first-order accurate in space ($\Delta x$) and time ($\Delta t$), and incurs numerical diffusion which can be quantified with the modified-equation analysis by substituting the discretised values other than $\psi _{i}^{n}$ in (5) with their Taylor expansions (using the Big O notation):

$${\begin{aligned}\psi _{i}^{n+1}&=\psi _{i}^{n}+{\frac {\left.\partial _{t}\psi \right|_{i}^{n}}{1!}}(\,+\Delta t\,)+{\frac {\left.\partial _{t}^{2}\psi \right|_{i}^{n}}{2!}}(\,+\Delta t\,)^{2}+O(\Delta t^{3})\\\psi _{i+1}^{n}&=\psi _{i}^{n}+{\frac {\left.\partial _{x}\psi \right|_{i}^{n}}{1!}}(+\Delta x)+{\frac {\left.\partial _{x}^{2}\psi \right|_{i}^{n}}{2!}}(+\Delta x)^{2}+O(\Delta x^{3})\\\psi _{i-1}^{n}&=\psi _{i}^{n}+{\frac {\left.\partial _{x}\psi \right|_{i}^{n}}{1!}}(-\Delta x)+{\frac {\left.\partial _{x}^{2}\psi \right|_{i}^{n}}{2!}}(-\Delta x)^{2}+O(\Delta x^{3})\end{aligned}}$$ 6

which yields (up to second-order terms):

$${\begin{aligned}\left.\partial _{t}\psi \right|_{i}^{n}\Delta t+\underbrace {\partial _{t}^{2}\psi } _{u^{2}\partial _{x}^{2}\psi }|_{i}^{n}{\frac {\Delta t^{2}}{2}}&={\frac {C}{2}}\left(\psi _{i-1}^{n}-\psi _{i+1}^{n}\right)+{\frac {|C|}{2}}\left(\psi _{i-1}^{n}-2\psi _{i}^{n}+\psi _{i+1}^{n}\right)\\&=-\,C\,\Delta x\,\partial _{x}\left.\psi \right|_{i}^{n}+{\frac {|C|}{2}}\,\Delta x^{2}\,\partial _{x}^{2}\left.\psi \right|_{i}^{n}\end{aligned}}$$ 7

in which the highlighted second-order time derivative can be replaced with second-order spatial derivative using the Cauchy-Kovalevskaya procedure (i.e., by substituting eq. (2) into a time-derivative of eq. (2) giving $\partial _{t}^{2}\psi =-u\partial _{x}\partial _{t}\psi =u^{2}\partial _{x}^{2}\psi$) yielding the following, so-called, modified equation:

$$\partial _{t}\psi |_{i}^{n}+u\partial _{x}\psi |_{i}^{n}=\underbrace {\left(|u|{\frac {\Delta x}{2}}-u^{2}{\frac {\Delta t}{2}}\right)} _{k}\,\partial _{x}^{2}\psi |_{i}^{n}$$ 8

which confirms that, as $\Delta t\rightarrow 0$ and $\Delta x\rightarrow 0$, the upwind scheme approximates eq. (2), but the truncation error has a leading-order Fickian source term with a coefficient of numerical diffusion given by $k$.

The concept of MPDATA lies in reversing the effect of numerical diffusion by using the pseudo-velocity technique that allows to express diffusion terms in transport problems using an equation of the same form as (2). To this end, Smolarkiewicz introduced the anti-diffusive pseudo velocity (note the opposite sign stemming from, de facto, integrating backwards in time to reverse the effect of diffusion):

$$\partial _{t}\psi +\partial _{x}{\big (}\!\!\!\!\!\!\!\!\!\!\!\!\underbrace {k{\frac {\partial _{x}\psi }{\psi }}} _{\text{anti-diffusive velocity}}\!\!\!\!\!\!\!\!\!\!\!\!\psi {\big )}=0$$ 9

and proposed discretisation allowing to perform a corrective step using upwind integration with an anti-diffusive Courant number field:

$$C_{i+1/2}^{m+1}={\frac {\Delta t}{\Delta x}}k_{i+1/2}^{m}\left.{\frac {\partial _{x}\psi }{\psi }}\right|_{i+1/2}^{m}\approx {\begin{cases}0&{\text{ if }}\psi _{i+1}^{m}+\psi _{i}^{m}=0\\\left[|C_{i+1/2}^{m}|-(C_{i+1/2}^{m})^{2}\right]{\frac {\psi _{i+1}^{m}-\psi _{i}^{m}}{\psi _{i+1}^{m}+\psi _{i}^{m}}}&{\text{ otherwise }}\end{cases}}$$ 10

where $m>0$ numbers algorithm iterations. The first pass of the scheme is an ordinary upwind integration using $C^{m=1}=C$ and yielding $\psi ^{m=1}$. In the second pass, the values of the anti-diffusive Courant number $C^{m=2}$ are computed based on $C^{m=1}$ and $\psi ^{m=1}$, and used to perform an upwind anti-diffusive pass yielding $\psi ^{m=2}$. Note that even for constant-in-$x$ physical $C^{m=1}$, the anti-diffusive$C^{m=1}$ varies in $x$ (and thus corresponds to a divergent flow). Subsequent iterations are correcting the integration of the anti-diffusive corrections from previous passes. If $M$ iterations are used, $\psi ^{m=M}$ is used as the result of integration of the advective term over one $\Delta t$.

TODO: A bounded

TODO: (note: multiple dimensions)

TODO: iterations

TODO: stability condition

TODO: positive definiteness, etc

TODO: conservativeness

TODO: anti-diffusive velocity is not solenoidal!

Minimal implementation and convergence analysis in Python

...

Algorithm variants and techniques used in concert with MPDATA

The basic scheme... + examples generated with PyMPDATA ?

• double-pass donor cell ...

• infinite-gauge variant ...

• non-oscillatory variant ... flux limiter...

• divergent-flow correction ...

• higher-order terms ...

• application to advection-diffusion problems

• application to inhomogeneous problems

• application to problems featuring pressure-gradient source terms (e.g., Boussinesq approximation (buoyancy) of the Navier-Stokes equations)

• immersed boundary method ...

• adaptive mesh refinement controlled by antidiffusive velocities

Algorithm steps (shallow-water system example)

The algorithm steps listed below constitute a solver to a minimal fluid dynamics problem - an inviscid shallow water equations system (in two spatial dimensions, assuming flat bathymetry):

$${\begin{cases}\,\,\,\partial _{t}h\!\!\!\!\!&\,\,+\,\nabla \cdot \,\,({\mathbf {v}}h)\!\!\!&=0\\\partial _{t}(uh)\!\!\!\!\!&\,\,+\,\nabla \cdot ({\mathbf {v}}uh)\!\!\!&=...\\\partial _{t}(vh)\!\!\!\!\!&\,\,+\,\nabla \cdot ({\mathbf {v}}vh)\!\!\!&=...\\\end{cases}}$$ N

where ${\mathbf {v}}=[u,v]$ is the advector (one vector field discretised at cell walls); fluid height $h$, momentum components $uh$ and $vh$ are the advectees (three scalar fields discretised at cell centers).

Since the scalar fields representing momentum components $(uh)$ and $(vh)$ are not of constant sign, it is recommended to use the infinite-gauge variant of MPDATA. This variant implies non-monotonic solutions, hence in practice it is only practical together with the non-oscillatory option used as well. Steps pertaining to both options are featured below.

An MPDATA-based solution consists of the following steps:

1. ...
2. ...

Applications

MPDATA has documented applications in the following domains:

• numerical weather prediction (NWP) systems:
  • compressible dynamical core for the COSMO system
  • nonhydrostatic finite-volume dynamical core for the ECMWF Integrated Forecast System (IFS)

• ocean modelling:
  • passive tracer advection in terrain-following coordinates in the Regional Ocean Modeling System (ROMS)

• large eddy simulations:
  • moist atmospheric flows (including representation of cloud and precipitation processes)
  • multi-phase particle-laden flows (e.g., contrail formation, aerosol-cloud-precipitation interactions)
  • turbulent magnetohydrodynamics (including Hall magnetohydrodynamics)
  • simulations of wind-driven sand dune dynamics
  • simulations of flows past complex structures (e.g. the Pentagon building and the Palace of Culture and Science in Warsaw)

• direct numerical simulations:
  • multi-phase flows (fog formation)

• solvers for shallow water equations resolving:
  • double-gyre circulation
  • wetting–drying areas

• aerosol size-spectrum dynamics:
  • advection problem in particle size dimension
  • two-dimensional advection in size-spectral and spatial dimensions

• mathematical finance:
  • solutions to the Black-Scholes-type equations (a stochastic differential equation transformed into a partial differential equation using Itô's lemma) cast as an advective transport problem for pricing of European and American options

Open-source implementations

• reusable MPDATA libraries and packages:
  • libmpdata++: structured-grid MPDATA implemented in C++ using Blitz++ for array handling and array-valued expressions, and OpenMP for multi-threading. API and usage examples documented in a journal article. Featuring distributed-memory parallelism using Message Passing Interface (MPI). Supports integration in 1-, 2- and 3-dimensional domains. Implements arbitrary number of corrective iterations, the infinite-gauge, non-oscillatory and fully-third-order variants. Features implementations of a shallow-water equation system solver, and a conjugate-residual Poisson equation solver for handling pressure terms in momentum conservation equations. The library has been developed as a dynamical core for the UWLCM large eddy simulation system. Released under the GNU GPL v3 license.
  • PyMPDATA structured-grid MPDATA implemented in Python using NumPy and Numba for array handling, just-in-time compilation and multi-threading, supporting distributed-memory parallelism using MPI via Numba-MPI. Supports integration in 1-, 2- and 3-dimensional domains, arbitrary number of corrective iterations, implements infinite-gauge, non-oscillatory and double-pass-donor-cell algorithm variants. Pure-Python wheels available on PyPI, release tarballs archived at Zenodo. Released under the GNU GPL v3 license together with a set of examples in a form of Jupyter notebooks. Web-based code-coupled documentation in English.

• MPDATA implementations integrated in other software:
  • AtmosFOAM (a derivative of OpenFOAM; OpenFOAM is distributed under the GPL v3 license) features implementations of MPDATA in both C++ and Python/Numba.
  • mpdat_2d: a FORTRAN 77 subroutine implementing 2-dimensional structured-grid MPDATA including the non-oscillatory option (likely a derivative of the original MPDATA implementation created at NCAR) has been publicly released as a part of babyEULAG code. No documentation or README, unspecified license and authorship.
  • smolar: structured-grid MPDATA implemented as a C library supporting multi-threading and multi-process parallelism, as well as offloading of anti-diffusive velocity computation to GPU using CUDA. Code publicly available on GitHub without release versioning, unspecified license, documentation in Russian language credited to Russian Academy of Sciences.
  • advectionHPCtester: structured-grid 3-dimensional implementation of the basic scheme in FORTRAN. Code archived on Zenodo. Released under the GNU LGPL.