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:

{{NumBlk|::| $\partial_t (G\psi) + \nabla \cdot (\bold{v}\psi) = GR,$ |{{EquationRef|1}}}}

where $\psi(\bold{x}, t)$ is an advected scalar field (or advectee); $\bold{v}=G\dot{\bold{x}}$ is the flow velocity vector field (or advector), $G(\bold{x})$ may either play the role of the fluid density, the Jacobian of coordinate transformation from Cartesian to curvilinear framework $\bold{x}$ , or their product; $R(\bold{x})$ combines all source terms{{Cite journal|title=MPDATA: A Finite-Difference Solver for Geophysical Flows|date=1998|journal=Journal of Computational Physics|volume=140|author=Smolarkiewicz, P.K. and Margolin, L.G.|issue=2 |pages=459–480|doi=10.1006/jcph.1998.5901|bibcode=1998JCoPh.140..459S }}.

MPDATA stands for Multidimensional Positive Definite Advection Transport Algorithm.

The algorithm was formulated by Piotr K. Smolarkiewicz{{Cite journal|author=Smolarkiewicz, P.K.|title=A simple positive definite advection scheme with small implicit diffusion|journal=Monthly Weather Review|volume=111|pages=479–186|date=1983|issue=3 |doi=10.1175/1520-0493(1983)111<0479:ASPDAS>2.0.CO;2|bibcode=1983MWRv..111..479S }}{{Cite journal|date=1984|journal=Journal of Computational Physics|author=Smolarkiewicz, P.K.|title=A fully multidimensional positive definite advection transport algorithm with small implicit diffusion|doi=10.1016/0021-9991(84)90121-9|volume=54|issue=2 |pages=325–362 |bibcode=1984JCoPh..54..325S }} (and is also referred to as Smolarkiewicz's method{{Cite book|author=Williams, M.M.R. and Loyalka, S.K.|title=Aerosol Science; theory and practice: with special applications to the nuclear industry|url=https://archive.org/details/aerosol-science-theory-and-practice-williams-loyalka|date=1991|publisher=Pergamon Press|isbn=0-08-037209-0}}) 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 $\bold{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{{Cite journal|title=The modified equation approach to the stability and accuracy analysis of finite-difference methods|author=Warming, R.F. and Hyett, B.J.|doi=10.1016/0021-9991(74)90011-4|volume=14|journal=Journal of Computational Physics|date=1974 |issue=2 |pages=159–179 |bibcode=1974JCoPh..14..159W }} 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{{Cite journal|title=MPDATA: Third-order accuracy for variable flows|author=Waruszewski, M. and Kühnlein, C. and Pawlowska H. and Smolarkiewicz, P.K.|doi=10.1016/j.jcp.2018.01.005|journal=Journal of Computational Physics|date=2018|volume=359|pages=361–379 |bibcode=2018JCoPh.359..361W }}.

While the original formulation of MPDATA employed structured grids, the algorithm has been subsequently also formulated for unstructured grids{{Cite journal|title=MPDATA: An edge-based unstructured-grid formulation|author=Smolarkiewicz, P.K. and Szmelter, J.|date=2005|journal=Journal of Computational Physics|volume=206|issue=2 |pages=624–649 |doi=10.1016/j.jcp.2004.12.021|bibcode=2005JCoPh.206..624S }}.

Despite being formulated for and named in reference to advection problems, as any other advection numerical scheme, MPDATA also applies to solutions of advection-diffusion as well as diffusion-only problems if Fickian diffusive terms are expressed as advective fluxes{{Cite journal|author=Smolarkiewicz, P.K. and Clark, T.|doi=10.1016/0021-9991(86)90270-6|journal=Journal of Computational Physics|volume=67|title=The multidimensional positive definite advection transport algorithm: Further development and applications|date=1986|issue=2 |pages=396–438 |bibcode=1986JCoPh..67..396S }} (approach referred to as the pseudo-velocity technique{{Citation|author=Lange, R.|title=ADPIC: a three-dimensional computer code for the study of pollutant dispersal and deposition under complex conditions|doi=10.2172/4308175|date=1973|url=https://digital.library.unt.edu/ark:/67531/metadc1022499/ }}{{Cite journal|doi=10.1175/1520-0450(1978)017<0320:ATDPIC>2.0.CO;2|author=Lange, R.|title=ADPIC–A Three-Dimensional Particle-in-Cell Model for the Dispersal of Atmospheric Pollutants and its Comparison to Regional Tracer Studies|journal=Journal of Applied Meteorology|date=1978 |volume=17|issue=3|page=320|bibcode=1978JApMe..17..320L }}).

== 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{{Cite journal|date=2006|doi=10.1002/fld.1071|volume=50|author=Smolarkiewicz, P.K.|title=Multidimensional positive definite advection transport algorithm: An overview|journal=International Journal for Numerical Methods in Fluids|issue=10 |pages=1123–1144 |bibcode=2006IJNMF..50.1123S }} (i.e., $\bold{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 \ge 0$ ), with the following flux form of the advection equation:

{{NumBlk|::| $\partial_t \psi + \partial_x (u\psi) = 0.$ |{{EquationRef|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:

{{NumBlk|::| $\frac{\psi^{n+1}_i - \psi^{n}_i}{\Delta t} +\frac{\overbrace{f(\psi_i^n, \psi^n_{i+1}, u^n_{i+1/2})}^\text{right-hand wall flux} - \overbrace{f(\psi_{i-1}^n, \psi^n_{i}, u^n_{i-1/2})}^\text{left-hand wall flux}}{\Delta x}=0$ |{{EquationRef|3}}}}

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

{{NumBlk|::| $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.$ |{{EquationRef|4}}}}

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

{{NumBlk|::| $\psi_i^{n+1} = \psi_i^{n} - \left[
\frac{C+|C|}{2}(\psi^n_i-\psi^n_{i-1}) + \frac{C-|C|}{2}(\psi^n_{i+1}-\psi^n_i)
\right]$ |{{EquationRef|5}}}}

which is conservative, sing-preserving and stable for $|C| \le 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 ({{EquationNote|5}}) with their Taylor expansions (using the Big O notation):

{{NumBlk|::| $\begin{align}
\psi_i^{n+1} &= \psi_i^{n} + \frac{\left.\partial_t\psi\right|_i^n}{1!} (\,+\Delta t\,) + \frac{\left.\partial^2_t\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^2_x\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^2_x\psi\right|_i^n}{2!} (-\Delta x)^2 + O(\Delta x^3)
\end{align}$ |{{EquationRef|6}}}}

which yields (up to second-order terms):

{{NumBlk|::| $\begin{align}
\left.\partial_t\psi\right|_i^n \Delta t + \underbrace{\partial^2_t\psi}_{u^2\partial_x^2\psi}|_i^n \frac{\Delta t^2}{2} &= \frac{C}{2} \left( \psi^n_{i-1} - \psi^n_{i+1}\right) + \frac$

{2} \left( \psi^n_{i-1} - 2 \psi_i^n + \psi^n_{i+1}\right)\\ &=- \,C\, \Delta x \, \partial_x \left.\psi\right|^n_i + \frac

{2} \,\Delta x^2 \, \partial_x^2\left.\psi\right|^n_i \end{align} |{{EquationRef|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. ({{EquationNote|2}}) into a time-derivative of eq. ({{EquationNote|2}}) giving $\partial^2_t \psi = -u \partial_x \partial_t \psi = u^2 \partial^2_x \psi$ ) yielding the following, so-called, modified equation:

{{NumBlk|::| $\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$ |{{EquationRef|8}}}}

which confirms that, as $\Delta t \rightarrow 0$ and $\Delta x \rightarrow 0$ , the upwind scheme approximates eq. ({{EquationNote|2}}), but the truncation error has a leading-order Fickian source term with a coefficient of numerical diffusion given by $k$ {{Cite journal|doi=10.2307/2003507|author=Roberts, K.V. and Weiss, N.O.|title=Convective Difference Schemes|volume=20|date=1966|journal=Mathematics of Computation|issue=94 |pages=272–299 |jstor=2003507 }}.

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 ({{EquationNote|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):

{{NumBlk|::| $\partial_t\psi + \partial_x \big(\!\!\!\!\!\!\!\!\!\!\!\!\underbrace{k \frac{\partial_x \psi}{\psi}}_\text{anti-diffusive velocity} \!\!\!\!\!\!\!\!\!\!\!\!\psi\big) = 0$ |{{EquationRef|9}}}}

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

{{NumBlk|::| $C^{m+1}_{i+ 1/2 } = \frac{\Delta t}{\Delta x} k^m_{i+ 1/2} \left.\frac{\partial_x \psi}{\psi}\right|^m_{i+ 1/2} \approx \begin{cases}
0 &\text{ if } \psi^m_{i+1}+\psi^m_{i} = 0\\
\left[|C^m_{i+1/2}|-(C^m_{i+1/2})^2\right]\frac{\psi^m_{i+1}-\psi^m_i}{\psi^m_{i+1}+\psi^m_i} &\text{ otherwise }
\end{cases}$ |{{EquationRef|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: stencil (numerical analysis) (and how it changes with iterations)

TODO: stability condition

TODO: positive definiteness, etc

TODO: conservativeness

TODO: anti-diffusive velocity is not solenoidal!

TODO: Last but not least ... bounadry conditions ...

Minimal implementation and convergence analysis in Python

...

Algorithm variants and techniques used in concert with MPDATA

The basic scheme... (note that many of the non-basic variants of the algorithm have larger stencils).

+ examples generated with PyMPDATA ?

double-pass donor cell ...{{Cite conference|author=Beason, C.W. and Margolin, L.G.|title=DPDC (double-pass donor cell): A second-order monotone scheme for advection|url=https://www.osti.gov/biblio/7049237|date=1988|osti=7049237 |conference=Nuclear explosives code developer's conference, Boulder, CO, USA, 11 Oct 1988}}

infinite-gauge variant ...{{Cite journal|title=

MPDATA: gauge transformations, limiters and monotonicity

|doi=10.1002/fld.1070|date=2006|author=Margolin, L. and Shashkov, M.|journal=International Journal for Numerical Methods in Fluids|volume=50|issue=10 |pages=1193–1206 |bibcode=2006IJNMF..50.1193M }}

non-oscillatory variant ... flux limiter... {{Cite journal|author=Grabowski, W.W. and Smolarkiewicz, P.K.|title=Monotone finite-difference approximations to the advection-condensation problem|date=1990|journal=Monthly Weather Review|volume=118|issue=10 |page=2082 |doi=10.1175/1520-0493(1990)118<2082:MFDATT>2.0.CO;2|bibcode=1990MWRv..118.2082G }}

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 ...{{Cite journal|doi=10.1016/j.jcp.2007.08.005|journal=Journalof Computational Physics|author=Smolarkiewicz, P.K. and Sharman, R. and Weil, J. and Perry, S.G. and Heist, D. and Bowker, G.|title=Building resolving large-eddy simulations and comparison with wind tunnel experiments|date=2007 |volume=227|issue=1 |pages=633–653 |bibcode=2007JCoPh.227..633S }}

adaptive mesh refinement controlled by antidiffusive velocities{{Cite journal|doi=10.1002/fld.1118|author=Szmelter, J. and Smolarkiewicz, P.K.|title=MPDATA error estimator for mesh adaptivity|journal=International Journal for Numerical Methods in Fluids|date=2006 |volume=50|issue=10 |pages=1269–1293 |bibcode=2006IJNMF..50.1269S }}

Algorithm steps ([[shallow water equations|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):

{{NumBlk|::| $\begin{cases}
\,\,\, \partial_t h \!\!\!\!\!&\,\, +\, \nabla\cdot \,\,(\bold{v} h) \!\!\!&= 0 \\
\partial_t (uh) \!\!\!\!\!&\,\, +\, \nabla\cdot (\bold{v} uh) \!\!\!&= ... \\
\partial_t (vh) \!\!\!\!\!&\,\, +\, \nabla\cdot (\bold{v} vh) \!\!\!&= ... \\
\end{cases}$ |{{EquationRef|N}}}}

where $\bold{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:

== Applications ==

Equation ({{EquationNote|1}}), which is numerically approximated with MPDATA, is used in modelling a wide range of transport phenomena across different scales and flow regimes. MPDATA has documented applications in the following domains:

numerical weather prediction (NWP) systems:
compressible dynamical core for the COSMO system{{Cite journal|title=Compressible EULAG Dynamical Core in COSMO: Convective-Scale Alpine Weather Forecasts|author=Ziemiański, M.Z. and Wójcik, D.K. and Rosa, B. and Piotrowski, Z.|journal=Monthly Weather Review|date=2021 |volume=149|issue=10 |pages=3563–3583 |doi=10.1175/MWR-D-20-0317.1|bibcode=2021MWRv..149.3563Z }}
nonhydrostatic finite-volume dynamical core for the ECMWF Integrated Forecast System (IFS){{Cite journal|title=FVM 1.0: a nonhydrostatic finite-volume dynamical core for the IFS|author=Kühnlein, C. and Deconinck, W. and Klein, R. and Malardel, S. and Piotrowski, Z.P. and Smolarkiewicz, P.K. and Szmelter, J. and Wedi, N.P.|doi=10.5194/gmd-12-651-2019|journal=Geoscientific Model Development|date=2019 |volume=12|issue=2 |pages=651–676 |doi-access=free |bibcode=2019GMD....12..651K }}

ocean modelling:
passive tracer advection in terrain-following coordinates in the Regional Ocean Modeling System (ROMS){{Cite journal|title=Developments in terrain-following ocean models: intercomparisons of numerical aspects|author=Ezer, T. and Arango, H. and Shchepetkin, A.F.|doi=10.1016/s1463-5003(02)00003-3|journal=Ocean Modelling|volume=4|date=2002|issue=3–4 |pages=249–267 |bibcode=2002OcMod...4..249E }}

large eddy simulations:
moist atmospheric flows (including representation of cloud and precipitation processes){{Cite journal|author=Pedersen, J.G. and Malinowski, S.P. and Grabowski, W.W.|title=Resolution and domain-size sensitivity in implicit large-eddy simulation of the stratocumulus-topped boundary layer|doi=10.1002/2015MS000572|journal=Journal of Advances in Modeling Earth Systems|volume=8|year=2016|issue=2 |pages=885–903 |bibcode=2016JAMES...8..885P }}
multi-phase particle-laden flows (e.g., contrail formation{{Cite journal|doi=10.1017/aer.2020.3|author=Unterstrasser, S. and Stephan, A.|title=Far field wake vortex evolution of two aircraft formation flight and implications on young contrails|volume=124|date=2020|journal=The Aeronautical Journal|issue=1275 |pages=667–702 }}, aerosol-cloud-precipitation interactions{{Cite journal|doi=10.5194/gmd-17-759-2024|journal=Geoscientific Model Development|volume=17|author=Zmijewski, P. and Dziekan, P. and Pawlowska, H.|title=Modeling collision–coalescence in particle microphysics: numerical convergence of mean and variance of precipitation in cloud simulations using the University of Warsaw Lagrangian Cloud Model (UWLCM) 2.1|date=2024|issue=2 |pages=759–780 |doi-access=free |bibcode=2024GMD....17..759Z }})
turbulent magnetohydrodynamics{{Cite journal|journal=Journal of Computational Physics|title=EULAG, a computational model for multiscale flows: An MHD extension|author=Smolarkiewicz, P.K. and Charbonneau, P.|date=2013|volume=236|pages=608–623 |doi=10.1016/j.jcp.2012.11.008|bibcode=2013JCoPh.236..608S }} (including Hall magnetohydrodynamics{{Cite journal|author=Bora, K. and Bhattacharyya, R. and Smolarkiewicz, P.LK.|title=Evolution of Three-dimensional Coherent Structures in Hall Magnetohydrodynamics|journal=The Astrophysical Journal|volume=906|date=2021|issue=2 |page=102 |doi=10.3847/1538-4357/abc8f7|doi-access=free |arxiv=2011.04223 |bibcode=2021ApJ...906..102B }})
simulations of wind-driven sand dune dynamics{{Cite journal|doi=10.1002/fld.1138|date=2006|title=Numerical simulation of sand dune evolution in severe winds|author=Ortiz, P. and Smolarkiewicz, P.K.|journal=International Journal for Numerical Methods in Fluids|volume=50|issue=10 |pages=1229–1246 |bibcode=2006IJNMF..50.1229O }}
simulations of flows past complex structures (e.g. the Pentagon building and the Palace of Culture and Science in Warsaw{{Cite journal|doi=10.1007/s10652-016-9470-3|author=Korycki, M. and Łobocki, L. and Wyszogrodzki, A.|title=Numerical simulation of stratified flow around a tall building of a complex shape|journal=Environmental Fluid Mechanics|volume=16|date=2016|issue=6 |pages=1143–1171 |bibcode=2016EFM....16.1143K }})

direct numerical simulations:
multi-phase flows (fog formation){{Cite journal|title=Direct Numerical Simulation of the Moist Stably Stratified Surface Layer: Turbulence and Fog Formation|author=Michael MacDonald, M. and Kurowski, M.J. and Teixeira, J.|date=2020|journal=Boundary-Layer Meteorology|volume=175|issue=3 |pages=343–368 |doi=10.1007/s10546-020-00511-2|arxiv=2012.04128 |bibcode=2020BoLMe.175..343M }}

solvers for shallow water equations resolving:
double-gyre circulation{{Cite journal|author=Kuang, D. and Lee, L.|date=2015|title=A Conservative Formulation and a Numerical Algorithm for the Double-Gyre Nonlinear Shallow-Water Model|journal=Numerical Mathematics: Theory, Methods and Applications|volume=8|issue=4 |pages=634–650 |doi=10.4208/nmtma.2015.m1408}}
wetting–drying areas{{Cite journal|doi=10.1016/j.jcp.2015.02.003|author=Jarecka, D. and Jaruga, A. and Smolarkiewicz, P.K.|journal=Journal of Computational Physics|volume=289|date=2015|title=A spreading drop of shallow water|pages=53–61 |bibcode=2015JCoPh.289...53J }}

aerosol size-spectrum dynamics:
advection problem in particle size dimension{{Cite journal|author=Unterstrasser, S. and Hoffmann, F. and Lerch, M.|title=Collisional growth in a particle-based cloud microphysical model: insights from column model simulations using LCM1D (v1.0)|journal=Geoscientific Model Development|volume=13|date=2020|issue=11 |pages=5119–5145 |doi=10.5194/gmd-13-5119-2020|doi-access=free |bibcode=2020GMD....13.5119U }}
two-dimensional advection in size-spectral and spatial dimensions{{Cite journal|author=Olesik, M.A. and Banaśkiewicz, J. and Bartman, P. and Baumgartner, M. and Unterstrasser, S. and Arabas, S.|title=On numerical broadening of particle-size spectra: a condensational growth study using PyMPDATA 1.0|doi=10.5194/gmd-15-3879-2022|journal=Geoscientific Model Development|volume=15|date=2022|issue=9 |pages=3879–3899 |doi-access=free |bibcode=2022GMD....15.3879O }}

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{{Cite journal|date=2020|title=Derivative pricing as a transport problem: MPDATA solutions to Black–Scholes-type equations|author=Arabas, S. and Farhat, A.|journal=Journal of Computational and Applied Mathematics|volume=373|doi=10.1016/j.cam.2019.05.023}}

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{{Cite journal|title=libmpdata++ 1.0: a library of parallel MPDATA solvers for systems of generalised transport equations|author=Jaruga, A. and Arabas, S. and Jarecka, D. and Pawlowska, H. and Smolarkiewicz, P. and Waruszewski, M.|doi=10.5194/gmd-8-1005-2015|journal=Geoscientific Model Development|volume=8|date=2015|issue=4 |pages=1005–1032 |doi-access=free |arxiv=1407.1309 |bibcode=2015GMD.....8.1005J }}. Featuring distributed-memory parallelism using Message Passing Interface (MPI){{Cite journal|title=University of Warsaw Lagrangian Cloud Model (UWLCM) 2.0: adaptation of a mixed Eulerian–Lagrangian numerical model for heterogeneous computing clusters|author=Dziekan, P. and Zmijewski, P.|journal=Geoscientific Model Development|volume=15|date=2022|issue=11 |pages=4489–4501 |doi=10.5194/gmd-15-4489-2022|doi-access=free |bibcode=2022GMD....15.4489D }}. 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{{Cite journal|author=Dziekan, P. and Waruszewski, M. and Pawlowska, H.|title=University of Warsaw Lagrangian Cloud Model (UWLCM) 1.0: a modern large-eddy simulation tool for warm cloud modeling with Lagrangian microphysics|journal=Geosci. Model Dev.|date=2019 |volume=12|issue=6 |pages=2587–2606 |doi=10.5194/gmd-12-2587-2019|doi-access=free |bibcode=2019GMD....12.2587D }} 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{{Cite journal|title=PyMPDATA v1: Numba-accelerated implementation of MPDATA with examples in Python, Julia and Matlab|author=Bartman, P. and Banaśkiewicz, J. and Drenda, S. and Manna, M. and Olesik, M.A. and Rozwoda, P. and Sadowski, M. and Arabas, S.|doi=10.21105/joss.03896|journal=Journal of Open Source Software|volume=7|date=2022|issue=77 |page=3896 |bibcode=2022JOSS....7.3896B }}, supporting distributed-memory parallelism using MPI via Numba-MPI{{Cite journal|author=Derlatka, K. and Manna, M. and Bulenok, O. and Zwicker, D. and Arabas, S.|title=Numba-MPI v1.0: Enabling MPI communication within Numba/LLVM JIT-compiled Python code|journal=SoftwareX|volume=28|date=2024|doi=10.1016/j.softx.2024.101897|arxiv=2407.13712 |bibcode=2024SoftX..2801897D }}. 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{{Citation|title=PyMPDATA Zenodo releases|date=2024 |doi=10.5281/zenodo.6329302 |last1=Arabas |first1=Sylwester |last2=Banaśkiewicz |first2=Jakub |last3=Bartman |first3=Piotr |last4=Derlatka |first4=Kacper |last5=Drenda |first5=Szymon |last6=Manna |first6=Maciej |last7=Olesik |first7=Michael |last8=Magnuszewski |first8=Paweł |last9=Manna |first9=Maciej |last10=Rozwoda |first10=Paweł |last11=Sadowski |first11=Michał }}. 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{{Citation|url=https://open-atmos.github.io/PyMPDATA|title=PyMPDATA documentation}}.

MPDATA implementations integrated in other software:
PISM (Parallel Ice Sheet Model, GPL v3 license) includes implementation of the non-oscillatory variant of two-dimensional MPDATA in C++
AtmosFOAM (a derivative of OpenFOAM; OpenFOAM is distributed under the GPL v3 license) features implementations of MPDATA in both C++ and Python/Numba{{Citation|url=https://github.com/search?q=org%3AAtmosFOAM%20mpdata&type=code|title=GitHub code search for MPDATA in AtmosFOAM repositories}}.
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{{Citation|title=bE_SDs v0.1|date=2017|doi=10.5281/zenodo.1050260 |author=Grabowski, W.W.}} as a part of babyEULAG code{{Cite journal|author=Grabowski, W.W. and Dziekan, P. and Pawlowska, H.|title=Lagrangian condensation microphysics with Twomey CCN activation|journal=Geosci. Model Dev.|date=2018 |volume=11|issue=1 |pages=103–120 |doi=10.5194/gmd-11-103-2018|doi-access=free |bibcode=2018GMD....11..103G }}. 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{{Citation|title=apc-llc/mpdata|url=https://github.com/apc-llc/mpdata|author=Mikushin, D.|publisher= Applied Parallel Computing LLC}} 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{{Citation|doi=10.5281/zenodo.8178549|author=Piotrowski, Z.|title=piotrows/advectionHPCtester: Advection HPC tester: MPDATA eDSL (1.0.0)|date=2023 }}. Released under the GNU LGPL.

References

{{draft categories|1=Category:Computational fluid dynamics Category:Numerical differential equations Category:Iterative methods Category:Numerical analysis Category:Conservation equations}}