Review waiting, please be patient. This may take 3 months or more, since drafts are reviewed in no specific order. There are 4,419 pending submissions waiting for review. If the submission is accepted, then this page will be moved into the article space. If the submission is declined, then the reason will be posted here. In the meantime, you can continue to improve this submission by editing normally. Where to get help If you need help editing or submitting your draft, please ask us a question at the AfC Help Desk or get live help from experienced editors. These venues are only for help with editing and the submission process, not to get reviews. If you need feedback on your draft, or if the review is taking a lot of time, you can try asking for help on the talk page of a relevant WikiProject. Some WikiProjects are more active than others so a speedy reply is not guaranteed. How to improve a draft Wikipedia:Contributing to Wikipedia – a basic overview on how to edit Wikipedia. Help:Wikitext – how to use the markup Help:Referencing for beginners – how to include references Wikipedia:Article development – how to develop your article Wikipedia:Writing better articles – how to improve your article Wikipedia:Verifiability – make sure your article includes reliable third-party sources You can also browse Wikipedia:Featured articles and Wikipedia:Good articles to find examples of Wikipedia's best writing on topics similar to your proposed article. Improving your odds of a speedy review To improve your odds of a faster review, tag your draft with relevant WikiProject tags using the button below. This will let reviewers know a new draft has been submitted in their area of interest. For instance, if you wrote about a female astronomer, you would want to add the Biography, Astronomy, and Women scientists tags. Add tags to your draft Editor resources Easy tools: Citation bot (help) | Advanced: Fix bare URLs Reviewer tools Instructions · What links here · MPDATA (talk: + · bio) · (log) · Copyvios report · reFill · Citation Bot · (Search: Google, Wikipedia) · Submitted 18 days ago by Slayoo (talk: D · +) · Last edited 18 days ago by Slayoo

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

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

1

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

MPDATA stands for Multidimensional Positive Definite Advection Transport Algorithm. The algorithm was formulated by Piotr K. Smolarkiewicz^[2][3] (and is also referred to as Smolarkiewicz's method^[4]) 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 ${\displaystyle {\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 can be derived through modified equation analysis^[5] 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 ${\displaystyle \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^[6].

While the original formulation of MPDATA employed structured grids, the algorithm has been subsequently also formulated for unstructured grids^[7]. 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^[8] (approach referred to as the pseudo-velocity technique^[9][10]).

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^[11] (i.e., ${\displaystyle {\mathbf {v}}=[u]}$ constant in time and space), without coordinate transformation (${\displaystyle G=1}$), for the case of homogeneous advection (${\displaystyle R=0}$) of a nonnegative scalar field (${\displaystyle \psi \geq 0}$), with the following flux form of the advection equation:

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

2

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

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

can be formulated with:

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

3

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

$${\displaystyle 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 ${\displaystyle C=u{\Delta t}/{\Delta x}}$, the resultant explicit-in-time scheme (referred to as "upwind", "upstream" or "donor-cell"), for a constant ${\displaystyle C}$ reads:

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

$${\displaystyle {\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):

$${\displaystyle {\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 ${\displaystyle \partial _{t}^{2}\psi =-u\partial _{x}\partial _{t}\psi =u^{2}\partial _{x}^{2}\psi }$) yielding the following, so-called, modified equation:

$${\displaystyle \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 ${\displaystyle \Delta t\rightarrow 0}$ and ${\displaystyle \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 ${\displaystyle k}$^[12].

The concept of MPDATA lies in reversing the effect of numerical diffusion by using the pseudo-velocity technique^[9][10] 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):

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

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

Note that despite a constant flow field, the anti-diffusive velocity varies in space and time and is not solenoidal.

The values of ${\displaystyle \psi _{i-1}^{n}}$ at the left edge of the domain, and ${\displaystyle \psi _{i+1}^{n}}$ at the right edge of the domain are obtained by specifying boundary conditions.

The basic formulation of the scheme has been extended yielding a multitude of MPDATA variants, including (in chronological order of literature references):

• Divergent-flow correction applicable for non-solenoidal flows (including the antidiffusive velocities)^[3],

• Double-Pass Donor Cell (DPDC) in which an analytical approximation of the sum of all corrective iterations is used to define a single corrective pass^[13],

• Non-oscillatory variant featuring a flux limiter to ensure monotone solutions^[14],

• Infinite-gauge variant for transport of variable-sign fields (e.g., momentum)^[15],

• third-order accurate anti-diffusive velocities for a wide range of flow types^[6].

Note that some of the non-basic variants of the algorithm have larger stencils) extent than the original formulation.

Some of the pertinent numerical techniques used together with MPDATA include:

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

• immersed boundary method^[16]

• adaptive mesh refinement controlled by antidiffusive velocities^[17].

Equation (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^[18]
  • nonhydrostatic finite-volume dynamical core for the ECMWF Integrated Forecast System (IFS)^[19]

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

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

• direct numerical simulations:
  • multi-phase flows (fog formation)^[28]

• solvers for shallow water equations resolving:
  • double-gyre circulation^[29]
  • wetting–drying areas^[30]

• aerosol size-spectrum dynamics:
  • advection problem in particle size dimension^[31]
  • two-dimensional advection in size-spectral and spatial dimensions^[32]

• 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^[33]

• 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^[34]. Featuring distributed-memory parallelism using Message Passing Interface (MPI)^[35]. 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^[36] 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^[37], supporting distributed-memory parallelism using MPI via Numba-MPI^[38]. 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^[39]. 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^[40].

• 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^[41].
  • 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^[42] as a part of babyEULAG code^[43]. 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^[44] 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^[45]. Released under the GNU LGPL.