PyMPDATA


NamePyMPDATA JSON
Version 1.1.6 PyPI version JSON
download
home_pageNone
SummaryNumba-accelerated Pythonic implementation of MPDATA with examples in Python, Julia and Matlab
upload_time2024-09-21 15:03:38
maintainerNone
docs_urlNone
authorhttps://github.com/open-atmos/PyMPDATA/graphs/contributors
requires_pythonNone
licenseGPL-3.0
keywords atmospheric-modelling numba numerical-integration advection pde-solver advection-diffusion
VCS
bugtrack_url
requirements No requirements were recorded.
Travis-CI No Travis.
coveralls test coverage No coveralls.
            # <img src="https://raw.githubusercontent.com/open-atmos/PyMPDATA/main/.github/pympdata_logo.png" width=100 height=113 alt="pympdata logo">

# PyMPDATA

[![Python 3](https://img.shields.io/static/v1?label=Python&logo=Python&color=3776AB&message=3)](https://www.python.org/)
[![LLVM](https://img.shields.io/static/v1?label=LLVM&logo=LLVM&color=gold&message=Numba)](https://www.numba.org)
[![Linux OK](https://img.shields.io/static/v1?label=Linux&logo=Linux&color=yellow&message=%E2%9C%93)](https://en.wikipedia.org/wiki/Linux)
[![macOS OK](https://img.shields.io/static/v1?label=macOS&logo=Apple&color=silver&message=%E2%9C%93)](https://en.wikipedia.org/wiki/macOS)
[![Windows OK](https://img.shields.io/static/v1?label=Windows&logo=Windows&color=white&message=%E2%9C%93)](https://en.wikipedia.org/wiki/Windows)
[![Jupyter](https://img.shields.io/static/v1?label=Jupyter&logo=Jupyter&color=f37626&message=%E2%9C%93)](https://jupyter.org/)
[![Maintenance](https://img.shields.io/badge/Maintained%3F-yes-green.svg)](https://github.com/open-atmos/PyMPDATA/graphs/commit-activity)
[![OpenHub](https://www.openhub.net/p/atmos-cloud-sim-uj-PyMPDATA/widgets/project_thin_badge?format=gif)](https://www.openhub.net/p/atmos-cloud-sim-uj-PyMPDATA)   
[![status](https://joss.theoj.org/papers/10e7361e43785dbb1b3d659c5b01757a/status.svg)](https://joss.theoj.org/papers/10e7361e43785dbb1b3d659c5b01757a)
[![DOI](https://zenodo.org/badge/225610671.svg)](https://zenodo.org/badge/latestdoi/225610671)     
[![EU Funding](https://img.shields.io/static/v1?label=EU%20Funding%20by&color=103069&message=FNP&logoWidth=25&logo=image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAC4AAAAeCAYAAABTwyyaAAAEzklEQVRYw9WYS2yUVRiGn3P5ZzozpZ3aUsrNgoKlKBINmkhpCCwwxIAhsDCpBBIWhmCMMYTEhSJ4i9EgBnSBEm81MRrFBhNXEuUSMCopiRWLQqEGLNgr085M5//POS46NNYFzHQ6qGc1i5nzP/P973m/9ztCrf7A8T9csiibCocUbvTzfxLcAcaM3cY3imXz25lT3Y34G7gQYAKV3+bFAHcATlBTPogJNADG92iY28FHW97kyPbnuW/W7xgzAhukQ9xe04PJeOT0HkQRwK0TlEeGWb/kOO9v3kdD3a8YK9GhDMfa6mg9fxunOm/lWPtcpDI4K7n/jnN8+uQbrFrUSiwU/DtSEUB/MsKKBT+zYslJqiYNgVE4JwhHkzy86wlWvrKVWDSZ/YFjZlU39yw4y/rGoyQGowWB67zl4QQue+jssMdXrQvZ/00jyeHwqCgDKwnsiJjSvkYAxsG5K9WsenYbJdqAtAjhCIxCSZt/4fK1w5A2WCvxrUAKCHwNVoA2aGmvq11jJQQapEXrgMBKqmJJugejKGWLIxXrBPFoigfv/omd675gRkU/xgqUDlAhH3UDaAAlLSqUQekAYyVTyhLs3tDMsvntlIYzOFcEcOcEGd9jx9oDbGs6QO0t/Tijxi9S4bhzxiWaVh5m94Zm0n7oui4ybo0raUlcncQnxx+g+WgDF/vLoYDmoqSl/dJUnt7XRCoTZjij0Z6Pc2LiNS4EBBkNvoeOJXN+yPWWSZeANOhwJq/98nKVwNdoL8B5AROxBKBL0gjh8DMhdCh3eJnrA0yqhLpplwmyup6IajvAOIGfKGVx3VmCRGnOMpe5QAdG0bT8CAeeep0d6z6nqjSJnQiZWEllLMWrmz6k+fE9rGk8MVqYgsGv5ZH2i1Opr+9kajzB5d74hKQ+KS3d/WVMLhtgdu1lriRiOR/4nDVunaR24x7qp3UV5Cb/fJvC83nv26W81LIK58SYNFmwq4hsGx/5BwKlzYRma2NUthgOJSew4i7ru9nJYCQF5tApb2yvjiDQKJV/IfJKh0o6qssSLKv/jcAoRKHQQzE2Lj2OMV5OkWFc4MZIpsev8uXWXRx6ZicbGk8QZLxxgwe+x/rlR3h3816+f2E7lbEU+ZDn3vKVpePCdFovzCISHqbl5EIoQOteKMPB1rto65zNyfOz+KOrGl06lHPQyi/WOohH0/T0l1MZH6A3GUEKl7Pmr2la6wBrBWWRDP2DUcqjKVKBGom9RZmABAykwnglafpSJSPQvsfiOR0EQ7ExVmazA8cY6N4K1iw6RdAXRwi4mgrheT5Dvs4LeuS81a15Ll/3dQisFVSVpnj7sf1sX/sZvhAc+6UOrQyBVUQ8gx/orFmDsZqtaw/y1qZ9zKjp5vDpenyjcNe+cLNmTiUdf/bEOddVQ0VpgsOn54ET+EYxvWKALSu+5tGG76it7MNaiZKGQ23zCIcMfUMxBnrjN3fmHHvCAlp+vJcXWx6itqoXpAEnUNLx8iMfo5Xh1i17R3PJYCpC2cZ3qK3sQ8WGEDDuXlAQuFKGHzpmopXhTNfk0bmxs7uC1w6uJul79AxFkMIiBJy5UoUWjrZLU5DCFdTARDHuDqVw+OkSwI0MCEW4gtNF2BPrBCo8fKNbtILWX9aUDqFqHnn7AAAAAElFTkSuQmCC)](https://www.fnp.org.pl/en/)
[![PL Funding](https://img.shields.io/static/v1?label=PL%20Funding%20by&color=d21132&message=NCN&logoWidth=25&logo=image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABQAAAANCAYAAACpUE5eAAAABmJLR0QA/wD/AP+gvaeTAAAAKUlEQVQ4jWP8////fwYqAiZqGjZqIHUAy4dJS6lqIOMdEZvRZDPcDQQAb3cIaY1Sbi4AAAAASUVORK5CYII=)](https://www.ncn.gov.pl/?language=en)
[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)

[![Github Actions Build Status](https://github.com/open-atmos/PyMPDATA/actions/workflows/tests+pypi.yml/badge.svg?branch=main)](https://github.com/open-atmos/PyMPDATA/actions)
[![Appveyor Build status](http://ci.appveyor.com/api/projects/status/github/open-atmos/PyMPDATA?branch=main&svg=true)](https://ci.appveyor.com/project/slayoo/pympdata/branch/main)
[![Coverage Status](https://codecov.io/gh/open-atmos/PyMPDATA/branch/main/graph/badge.svg)](https://app.codecov.io/gh/open-atmos/PyMPDATA)

[![PyPI version](https://badge.fury.io/py/PyMPDATA.svg)](https://pypi.org/project/PyMPDATA)
[![API docs](https://shields.mitmproxy.org/badge/docs-pdoc.dev-brightgreen.svg)](https://open-atmos.github.io/PyMPDATA/index.html)

PyMPDATA is a high-performance Numba-accelerated Pythonic implementation of the MPDATA 
  algorithm of Smolarkiewicz et al. used in geophysical fluid dynamics and beyond.
MPDATA numerically solves generalised transport equations -
  partial differential equations used to model conservation/balance laws, scalar-transport problems,
  convection-diffusion phenomena.
As of the current version, PyMPDATA supports homogeneous transport
  in 1D, 2D and 3D using structured meshes, optionally
  generalised by employment of a Jacobian of coordinate transformation. 
PyMPDATA includes implementation of a set of MPDATA variants including
  the non-oscillatory option, infinite-gauge, divergent-flow, double-pass donor cell (DPDC) and 
  third-order-terms options.
It also features support for integration of Fickian-terms in advection-diffusion
  problems using the pseudo-transport velocity approach.
In 2D and 3D simulations, domain-decomposition is used for multi-threaded parallelism.

PyMPDATA is engineered purely in Python targeting both performance and usability,
    the latter encompassing research users', developers' and maintainers' perspectives.
From researcher's perspective, PyMPDATA offers hassle-free installation on multitude
  of platforms including Linux, OSX and Windows, and eliminates compilation stage
  from the perspective of the user.
From developers' and maintainers' perspective, PyMPDATA offers a suite of unit tests, 
  multi-platform continuous integration setup,
  seamless integration with Python development aids including debuggers and profilers.

PyMPDATA design features
  a custom-built multi-dimensional Arakawa-C grid layer allowing
  to concisely represent multi-dimensional stencil operations on both
  scalar and vector fields.
The grid layer is built on top of NumPy's ndarrays (using "C" ordering)
  using the Numba's ``@njit`` functionality for high-performance array traversals.
It enables one to code once for multiple dimensions, and automatically
  handles (and hides from the user) any halo-filling logic related with boundary conditions.
Numba ``prange()`` functionality is used for implementing multi-threading 
  (it offers analogous functionality to OpenMP parallel loop execution directives).
The Numba's deviation from Python semantics rendering [closure variables
  as compile-time constants](https://numba.pydata.org/numba-doc/dev/reference/pysemantics.html#global-and-closure-variables)
  is extensively exploited within ``PyMPDATA``
  code base enabling the just-in-time compilation to benefit from 
  information on domain extents, algorithm variant used and problem
  characteristics (e.g., coordinate transformation used, or lack thereof).

A separate project called [``PyMPDATA-MPI``](https://github.com/open-atmos/PyMPDATA-MPI) 
  depicts how [``numba-mpi``](https://pypi.org/project/numba-mpi) can be used
  to enable distributed memory parallelism in PyMPDATA.

The [``PyMPDATA-examples``](https://pypi.org/project/PyMPDATA-examples/) 
  package covers a set of examples presented in the form of Jupyer notebooks
  offering single-click deployment in the cloud using [mybinder.org](https://mybinder.org)
  or using [colab.research.google.com](https://colab.research.google.com/).
The examples reproduce results from several published
  works on MPDATA and its applications, and provide a validation of the implementation
  and its performance.
 
## Dependencies and installation

To install PyMPDATA, one may use: ``pip install PyMPDATA`` (or 
``pip install git+https://github.com/open-atmos/PyMPDATA.git`` to get updates beyond the latest release).
PyMPDATA depends on ``NumPy`` and ``Numba``.

Running the tests shipped with the package requires additional packages listed in the 
[test-time-requirements.txt](https://github.com/open-atmos/PyMPDATA/blob/main/test-time-requirements.txt) file
(which include ``PyMPDATA-examples``, see below).

## Examples (Jupyter notebooks reproducing results from literature):

PyMPDATA examples are bundled with PyMPDATA and located in the `examples` subfolder.
They constitute a separate [``PyMPDATA_examples``](https://pypi.org/p/PyMPDATA-examples) Python package which is also available at PyPI.
The examples have additional dependencies listed in [``PyMPDATA_examples`` package ``setup.py``](https://github.com/open-atmos/PyMPDATA/blob/main/examples/setup.py) file.
Running the examples requires the ``PyMPDATA_examples`` package to be installed.
Since the examples package includes Jupyter notebooks (and their execution requires write access), the suggested install and launch steps are:
```
git clone https://github.com/open-atmos/PyMPDATA-examples.git
cd PyMPDATA-examples
pip install -e .
jupyter-notebook
```
Alternatively, one can also install the examples package from pypi.org by using ``pip install PyMPDATA-examples``.
  
## Package structure and API:

In short, PyMPDATA numerically solves the following equation:

$$ \partial_t (G \psi) + \nabla \cdot (Gu \psi) + \mu \Delta (G \psi) = 0 $$

where scalar field $\psi$ is referred to as the advectee,
vector field u is referred to as advector, and the G factor corresponds to optional coordinate transformation.
The inclusion of the Fickian diffusion term is optional and is realised through modification of the
advective velocity field with MPDATA handling both the advection and diffusion (for discussion
see, e.g. [Smolarkiewicz and Margolin 1998](https://doi.org/10.1006/jcph.1998.5901), sec. 3.5, par. 4).

The key classes constituting the PyMPDATA interface are summarised below with code
snippets exemplifying usage of PyMPDATA from Python, Julia and Matlab.

A [pdoc-generated](https://pdoc.dev/) documentation of PyMPDATA public API is maintained at: [https://open-atmos.github.io/PyMPDATA](  https://open-atmos.github.io/PyMPDATA/index.html) 

#### Options class

The [``Options``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/options.html) class
groups both algorithm variant options as well as some implementation-related
flags that need to be set at the first place. All are set at the time
of instantiation using the following keyword arguments of the constructor 
(all having default values indicated below):
- ``n_iters: int = 2``: number of iterations (2 means upwind + one corrective iteration)
- ``infinite_gauge: bool = False``: flag enabling the infinite-gauge option (does not maintain sign of the advected field, thus in practice implies switching flux corrected transport on)
- ``divergent_flow: bool = False``: flag enabling divergent-flow terms when calculating antidiffusive velocity
- ``nonoscillatory: bool = False``: flag enabling the non-oscillatory or monotone variant (a.k.a flux-corrected transport option, FCT)
- ``third_order_terms: bool = False``: flag enabling third-order terms
- ``epsilon: float = 1e-15``: value added to potentially zero-valued denominators 
- ``non_zero_mu_coeff: bool = False``: flag indicating if code for handling the Fickian term is to be optimised out
- ``DPDC: bool = False``: flag enabling double-pass donor cell option (recursive pseudovelocities)
- ``dimensionally_split: bool = False``: flag disabling cross-dimensional terms in antidiffusive velocity
- ``dtype: np.floating = np.float64``: floating point precision

For a discussion of the above options, see e.g., [Smolarkiewicz & Margolin 1998](https://doi.org/10.1006/jcph.1998.5901),
[Jaruga, Arabas et al. 2015](https://doi.org/10.5194/gmd-8-1005-2015) and [Olesik, Arabas et al. 2020](https://arxiv.org/abs/2011.14726)
(the last with examples using PyMPDATA).

In most use cases of PyMPDATA, the first thing to do is to instantiate the 
[``Options``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/options.html) class 
with arguments suiting the problem at hand, e.g.:
<details>
<summary>Julia code (click to expand)</summary>

```Julia
using Pkg
Pkg.add("PyCall")
using PyCall
Options = pyimport("PyMPDATA").Options
options = Options(n_iters=2)
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
Options = py.importlib.import_module('PyMPDATA').Options;
options = Options(pyargs('n_iters', 2));
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
from PyMPDATA import Options
options = Options(n_iters=2)
```
</details>

#### Arakawa-C grid layer and boundary conditions

In PyMPDATA, the solution domain is assumed to extend from the
first cell's boundary to the last cell's boundary (thus the
first scalar field value is at $\[\Delta x/2, \Delta y/2\]$.
The [``ScalarField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/scalar_field.html)
and [``VectorField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/vector_field.html) classes implement the
[Arakawa-C staggered grid](https://en.wikipedia.org/wiki/Arakawa_grids#Arakawa_C-grid) logic
in which:
- scalar fields are discretised onto cell centres (one value per cell),
- vector field components are discretised onto cell walls.

The schematic of the employed grid/domain layout in two dimensions is given below
(with the Python code snippet generating the figure as a part of CI workflow):
<details>
<summary>Python code (click to expand)</summary>

```Python
import numpy as np
from matplotlib import pyplot

dx, dy = .2, .3
grid = (10, 5)

pyplot.scatter(*np.mgrid[
        dx / 2 : grid[0] * dx : dx, 
        dy / 2 : grid[1] * dy : dy
    ], color='red', 
    label='scalar-field values at cell centres'
)
pyplot.quiver(*np.mgrid[
        0 : (grid[0]+1) * dx : dx, 
        dy / 2 : grid[1] * dy : dy
    ], 1, 0, pivot='mid', color='green', width=.005,
    label='vector-field x-component values at cell walls'
)
pyplot.quiver(*np.mgrid[
        dx / 2 : grid[0] * dx : dx,
        0: (grid[1] + 1) * dy : dy
    ], 0, 1, pivot='mid', color='blue', width=.005,
    label='vector-field y-component values at cell walls'
)
pyplot.xticks(np.linspace(0, grid[0]*dx, grid[0]+1))
pyplot.yticks(np.linspace(0, grid[1]*dy, grid[1]+1))
pyplot.title(f'staggered grid layout (grid={grid}, dx={dx}, dy={dy})')
pyplot.xlabel('x')
pyplot.ylabel('y')
pyplot.legend(bbox_to_anchor=(.1, -.1), loc='upper left', ncol=1)
pyplot.grid()
pyplot.savefig('readme_grid.png')
```
</details>

![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_grid.png)

The ``__init__`` methods of
[``ScalarField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/scalar_field.html)
and [``VectorField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/vector_field.html)
have the following signatures:
- [``ScalarField(data: np.ndarray, halo: int, boundary_conditions)``](https://github.com/open-atmos/PyMPDATA/blob/main/PyMPDATA/scalar_field.py)
- [``VectorField(data: Tuple[np.ndarray, ...], halo: int, boundary_conditions)``](https://github.com/open-atmos/PyMPDATA/blob/main/PyMPDATA/vector_field.py)
The ``data`` parameters are expected to be Numpy arrays or tuples of Numpy arrays, respectively.
The ``halo`` parameter is the extent of ghost-cell region that will surround the
data and will be used to implement boundary conditions. 
Its value (in practice 1 or 2) is
dependent on maximal stencil extent for the MPDATA variant used and
can be easily obtained using the ``Options.n_halo`` property.

As an example, the code below shows how to instantiate a scalar
and a vector field given a 2D constant-velocity problem,
using a grid of 24x24 points, Courant numbers of -0.5 and -0.25
in "x" and "y" directions, respectively, with periodic boundary 
conditions and with an initial Gaussian signal in the scalar field
(settings as in Fig.&nbsp;5 in [Arabas et al. 2014](https://doi.org/10.3233/SPR-140379)):
<details>
<summary>Julia code (click to expand)</summary>

```Julia
ScalarField = pyimport("PyMPDATA").ScalarField
VectorField = pyimport("PyMPDATA").VectorField
Periodic = pyimport("PyMPDATA.boundary_conditions").Periodic

nx, ny = 24, 24
Cx, Cy = -.5, -.25
idx = CartesianIndices((nx, ny))
halo = options.n_halo
advectee = ScalarField(
    data=exp.(
        -(getindex.(idx, 1) .- .5 .- nx/2).^2 / (2*(nx/10)^2) 
        -(getindex.(idx, 2) .- .5 .- ny/2).^2 / (2*(ny/10)^2)
    ),  
    halo=halo, 
    boundary_conditions=(Periodic(), Periodic())
)
advector = VectorField(
    data=(fill(Cx, (nx+1, ny)), fill(Cy, (nx, ny+1))),
    halo=halo,
    boundary_conditions=(Periodic(), Periodic())    
)
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
ScalarField = py.importlib.import_module('PyMPDATA').ScalarField;
VectorField = py.importlib.import_module('PyMPDATA').VectorField;
Periodic = py.importlib.import_module('PyMPDATA.boundary_conditions').Periodic;

nx = int32(24);
ny = int32(24);
  
Cx = -.5;
Cy = -.25;

[xi, yi] = meshgrid(double(0:1:nx-1), double(0:1:ny-1));

halo = options.n_halo;
advectee = ScalarField(pyargs(...
    'data', py.numpy.array(exp( ...
        -(xi+.5-double(nx)/2).^2 / (2*(double(nx)/10)^2) ...
        -(yi+.5-double(ny)/2).^2 / (2*(double(ny)/10)^2) ...
    )), ... 
    'halo', halo, ...
    'boundary_conditions', py.tuple({Periodic(), Periodic()}) ...
));
advector = VectorField(pyargs(...
    'data', py.tuple({ ...
        Cx * py.numpy.ones(int32([nx+1 ny])), ... 
        Cy * py.numpy.ones(int32([nx ny+1])) ...
     }), ...
    'halo', halo, ...
    'boundary_conditions', py.tuple({Periodic(), Periodic()}) ...
));
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
from PyMPDATA import ScalarField
from PyMPDATA import VectorField
from PyMPDATA.boundary_conditions import Periodic
import numpy as np

nx, ny = 24, 24
Cx, Cy = -.5, -.25
halo = options.n_halo

xi, yi = np.indices((nx, ny), dtype=float)
advectee = ScalarField(
  data=np.exp(
    -(xi+.5-nx/2)**2 / (2*(nx/10)**2)
    -(yi+.5-ny/2)**2 / (2*(ny/10)**2)
  ),
  halo=halo,
  boundary_conditions=(Periodic(), Periodic())
)
advector = VectorField(
  data=(np.full((nx + 1, ny), Cx), np.full((nx, ny + 1), Cy)),
  halo=halo,
  boundary_conditions=(Periodic(), Periodic())
)
```
</details>

Note that the shapes of arrays representing components 
of the velocity field are different than the shape of
the scalar field array due to employment of the staggered grid.

Besides the exemplified [``Periodic``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/periodic.html) class representing 
periodic boundary conditions, PyMPDATA supports 
[``Extrapolated``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/extrapolated.html), 
[``Constant``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/constant.html) and
[``Polar``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/polar.html) 
boundary conditions.

#### Stepper

The logic of the MPDATA iterative solver is represented
in PyMPDATA by the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) class.

When instantiating the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html), the user has a choice 
of either supplying just the  number of dimensions or specialising the stepper for a given grid:
<details>
<summary>Julia code (click to expand)</summary>

```Julia
Stepper = pyimport("PyMPDATA").Stepper

stepper = Stepper(options=options, n_dims=2)
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
Stepper = py.importlib.import_module('PyMPDATA').Stepper;

stepper = Stepper(pyargs(...
  'options', options, ...
  'n_dims', int32(2) ...
));
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
from PyMPDATA import Stepper

stepper = Stepper(options=options, n_dims=2)
```
</details>
or
<details>
<summary>Julia code (click to expand)</summary>

```Julia
stepper = Stepper(options=options, grid=(nx, ny))
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
stepper = Stepper(pyargs(...
  'options', options, ...
  'grid', py.tuple({nx, ny}) ...
));
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
stepper = Stepper(options=options, grid=(nx, ny))
```
</details>

In the latter case, noticeably 
faster execution can be expected, however the resultant
stepper is less versatile as bound to the given grid size.
If number of dimensions is supplied only, the integration
might take longer, yet same instance of the
stepper can be used for different grids.  

Since creating an instance of the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) class
involves time-consuming compilation of the algorithm code,
the class is equipped with a cache logic - subsequent
calls with same arguments return references to previously
instantiated objects. Instances of [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) contain no
mutable data and are (thread-)safe to be reused.

The init method of [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) has an optional
``non_unit_g_factor`` argument which is a Boolean flag 
enabling handling of the G factor term which can be used to 
represent coordinate transformations and/or variable fluid density. 

Optionally, the number of threads to use for domain decomposition
in the first (non-contiguous) dimension during 2D and 3D calculations
may be specified using the optional ``n_threads`` argument with a
default value of ``numba.get_num_threads()``. The multi-threaded
logic of PyMPDATA depends thus on settings of numba, namely on the
selected threading layer (either via ``NUMBA_THREADING_LAYER`` env 
var or via ``numba.config.THREADING_LAYER``) and the selected size of the 
thread pool (``NUMBA_NUM_THREADS`` env var or ``numba.config.NUMBA_NUM_THREADS``).


#### Solver

Instances of the [``Solver``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html) class are used to control
the integration and access solution data. During instantiation, 
additional memory required by the solver is 
allocated according to the options provided. 

The only method of the [``Solver``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html) class besides the
init is [``advance(n_steps, mu_coeff, ...)``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html#Solver.advance) 
which advances the solution by ``n_steps`` timesteps, optionally
taking into account a given diffusion coefficient ``mu_coeff``.

Solution state is accessible through the ``Solver.advectee`` property.
Multiple solver[s] can share a single stepper, e.g., as exemplified in the shallow-water system solution in the examples package.

Continuing with the above code snippets, instantiating
a solver and making 75 integration steps looks as follows:
<details>
<summary>Julia code (click to expand)</summary>

```Julia
Solver = pyimport("PyMPDATA").Solver
solver = Solver(stepper=stepper, advectee=advectee, advector=advector)
solver.advance(n_steps=75)
state = solver.advectee.get()
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
Solver = py.importlib.import_module('PyMPDATA').Solver;
solver = Solver(pyargs('stepper', stepper, 'advectee', advectee, 'advector', advector));
solver.advance(pyargs('n_steps', 75));
state = solver.advectee.get();
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
from PyMPDATA import Solver

solver = Solver(stepper=stepper, advectee=advectee, advector=advector)
state_0 = solver.advectee.get().copy()
solver.advance(n_steps=75)
state = solver.advectee.get()
```
</details>

Now let's plot the results using `matplotlib` roughly as in Fig.&nbsp;5 in [Arabas et al. 2014](https://doi.org/10.3233/SPR-140379):

<details>
<summary>Python code (click to expand)</summary>

```Python
def plot(psi, zlim, norm=None):
    xi, yi = np.indices(psi.shape)
    fig, ax = pyplot.subplots(subplot_kw={"projection": "3d"})
    pyplot.gca().plot_wireframe(
        xi+.5, yi+.5, 
        psi, color='red', linewidth=.5
    )
    ax.set_zlim(zlim)
    for axis in (ax.xaxis, ax.yaxis, ax.zaxis):
        axis.pane.fill = False
        axis.pane.set_edgecolor('black')
        axis.pane.set_alpha(1)
    ax.grid(False)
    ax.set_zticks([])
    ax.set_xlabel('x/dx')
    ax.set_ylabel('y/dy')
    ax.set_proj_type('ortho') 
    cnt = ax.contourf(xi+.5, yi+.5, psi, zdir='z', offset=-1, norm=norm)
    cbar = pyplot.colorbar(cnt, pad=.1, aspect=10, fraction=.04)
    return cbar.norm

zlim = (-1, 1)
norm = plot(state_0, zlim)
pyplot.savefig('readme_gauss_0.png')
plot(state, zlim, norm)
pyplot.savefig('readme_gauss.png')
```
</details>

![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_gauss_0.png)    
![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_gauss.png)

#### Debugging

PyMPDATA relies heavily on Numba to provide high-performance 
number crunching operations. Arguably, one of the key advantage 
of embracing Numba is that it can be easily switched off. This
brings multiple-order-of-magnitude drop in performance, yet 
it also make the entire code of the library susceptible to
interactive debugging, one way of enabling it is by setting the 
following environment variable before importing PyMPDATA:
<details>
<summary>Julia code (click to expand)</summary>

```Julia
ENV["NUMBA_DISABLE_JIT"] = "1"
```
</details>
<details>
<summary>Matlab code (click to expand)</summary>

```Matlab
setenv('NUMBA_DISABLE_JIT', '1');
```
</details>
<details open>
<summary>Python code (click to expand)</summary>

```Python
import os
os.environ["NUMBA_DISABLE_JIT"] = "1"
```
</details>

## Contributing, reporting issues, seeking support 

Submitting new code to the project, please preferably use [GitHub pull requests](https://github.com/open-atmos/PyMPDATA/pulls) 
(or the [PyMPDATA-examples PR site](https://github.com/open-atmos/PyMPDATA-examples/pulls) if working on examples) - it helps to keep record of code authorship, 
track and archive the code review workflow and allows to benefit
from the continuous integration setup which automates execution of tests 
with the newly added code. 

As of now, the copyright to the entire PyMPDATA codebase is with the Jagiellonian
University, and code contributions are assumed to imply transfer of copyright.
Should there be a need to make an exception, please indicate it when creating
a pull request or contributing code in any other way. In any case, 
the license of the contributed code must be compatible with GPL v3.

Developing the code, we follow [The Way of Python](https://www.python.org/dev/peps/pep-0020/) and 
the [KISS principle](https://en.wikipedia.org/wiki/KISS_principle).
The codebase has greatly benefited from [PyCharm code inspections](https://www.jetbrains.com/help/pycharm/code-inspection.html)
and [Pylint](https://pylint.org) code analysis (Pylint checks are part of the
CI workflows).

Issues regarding any incorrect, unintuitive or undocumented bahaviour of
PyMPDATA are best to be reported on the [GitHub issue tracker](https://github.com/open-atmos/PyMPDATA/issues/new).
Feature requests are recorded in the "Ideas..." [PyMPDATA wiki page](https://github.com/open-atmos/PyMPDATA/wiki/Ideas-for-new-features-and-examples).

We encourage to use the [GitHub Discussions](https://github.com/open-atmos/PyMPDATA/discussions) feature
(rather than the issue tracker) for seeking support in understanding, using and extending PyMPDATA code.

Please use the PyMPDATA issue-tracking and dicsussion infrastructure for `PyMPDATA-examples` as well.
We look forward to your contributions and feedback.

## Credits:
Development of PyMPDATA was supported by the EU through a grant of the [Foundation for Polish Science](http://fnp.org.pl) (POIR.04.04.00-00-5E1C/18).

copyright: Jagiellonian University   
licence: GPL v3   

## Other open-source MPDATA implementations:
- mpdat_2d in babyEULAG (FORTRAN)
  https://github.com/igfuw/bE_SDs/blob/master/babyEULAG.SDs.for#L741
- mpdata-oop (C++, Fortran, Python)
  https://github.com/igfuw/mpdata-oop
- apc-llc/mpdata (C++)
  https://github.com/apc-llc/mpdata
- libmpdata++ (C++):
  https://github.com/igfuw/libmpdataxx
- AtmosFOAM:
  https://github.com/AtmosFOAM/AtmosFOAM/tree/947b192f69d973ea4a7cfab077eb5c6c6fa8b0cf/applications/solvers/advection/MPDATAadvectionFoam

## Other Python packages for solving hyperbolic transport equations

- PyPDE: https://pypi.org/project/PyPDE/
- FiPy: https://pypi.org/project/FiPy/
- ader: https://pypi.org/project/ader/
- centpy: https://pypi.org/project/centpy/
- mattflow: https://pypi.org/project/mattflow/
- FastFD: https://pypi.org/project/FastFD/
- Pyclaw: https://www.clawpack.org/pyclaw/

            

Raw data

            {
    "_id": null,
    "home_page": null,
    "name": "PyMPDATA",
    "maintainer": null,
    "docs_url": null,
    "requires_python": null,
    "maintainer_email": null,
    "keywords": "atmospheric-modelling, numba, numerical-integration, advection, pde-solver, advection-diffusion",
    "author": "https://github.com/open-atmos/PyMPDATA/graphs/contributors",
    "author_email": "sylwester.arabas@agh.edu.pl",
    "download_url": "https://files.pythonhosted.org/packages/b1/93/6a989073744297766a6c32ac3f8751450378b39abcb235429af2f4521f98/pympdata-1.1.6.tar.gz",
    "platform": null,
    "description": "# <img src=\"https://raw.githubusercontent.com/open-atmos/PyMPDATA/main/.github/pympdata_logo.png\" width=100 height=113 alt=\"pympdata logo\">\n\n# PyMPDATA\n\n[![Python 3](https://img.shields.io/static/v1?label=Python&logo=Python&color=3776AB&message=3)](https://www.python.org/)\n[![LLVM](https://img.shields.io/static/v1?label=LLVM&logo=LLVM&color=gold&message=Numba)](https://www.numba.org)\n[![Linux OK](https://img.shields.io/static/v1?label=Linux&logo=Linux&color=yellow&message=%E2%9C%93)](https://en.wikipedia.org/wiki/Linux)\n[![macOS OK](https://img.shields.io/static/v1?label=macOS&logo=Apple&color=silver&message=%E2%9C%93)](https://en.wikipedia.org/wiki/macOS)\n[![Windows OK](https://img.shields.io/static/v1?label=Windows&logo=Windows&color=white&message=%E2%9C%93)](https://en.wikipedia.org/wiki/Windows)\n[![Jupyter](https://img.shields.io/static/v1?label=Jupyter&logo=Jupyter&color=f37626&message=%E2%9C%93)](https://jupyter.org/)\n[![Maintenance](https://img.shields.io/badge/Maintained%3F-yes-green.svg)](https://github.com/open-atmos/PyMPDATA/graphs/commit-activity)\n[![OpenHub](https://www.openhub.net/p/atmos-cloud-sim-uj-PyMPDATA/widgets/project_thin_badge?format=gif)](https://www.openhub.net/p/atmos-cloud-sim-uj-PyMPDATA)   \n[![status](https://joss.theoj.org/papers/10e7361e43785dbb1b3d659c5b01757a/status.svg)](https://joss.theoj.org/papers/10e7361e43785dbb1b3d659c5b01757a)\n[![DOI](https://zenodo.org/badge/225610671.svg)](https://zenodo.org/badge/latestdoi/225610671)     \n[![EU Funding](https://img.shields.io/static/v1?label=EU%20Funding%20by&color=103069&message=FNP&logoWidth=25&logo=image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAC4AAAAeCAYAAABTwyyaAAAEzklEQVRYw9WYS2yUVRiGn3P5ZzozpZ3aUsrNgoKlKBINmkhpCCwwxIAhsDCpBBIWhmCMMYTEhSJ4i9EgBnSBEm81MRrFBhNXEuUSMCopiRWLQqEGLNgr085M5//POS46NNYFzHQ6qGc1i5nzP/P973m/9ztCrf7A8T9csiibCocUbvTzfxLcAcaM3cY3imXz25lT3Y34G7gQYAKV3+bFAHcATlBTPogJNADG92iY28FHW97kyPbnuW/W7xgzAhukQ9xe04PJeOT0HkQRwK0TlEeGWb/kOO9v3kdD3a8YK9GhDMfa6mg9fxunOm/lWPtcpDI4K7n/jnN8+uQbrFrUSiwU/DtSEUB/MsKKBT+zYslJqiYNgVE4JwhHkzy86wlWvrKVWDSZ/YFjZlU39yw4y/rGoyQGowWB67zl4QQue+jssMdXrQvZ/00jyeHwqCgDKwnsiJjSvkYAxsG5K9WsenYbJdqAtAjhCIxCSZt/4fK1w5A2WCvxrUAKCHwNVoA2aGmvq11jJQQapEXrgMBKqmJJugejKGWLIxXrBPFoigfv/omd675gRkU/xgqUDlAhH3UDaAAlLSqUQekAYyVTyhLs3tDMsvntlIYzOFcEcOcEGd9jx9oDbGs6QO0t/Tijxi9S4bhzxiWaVh5m94Zm0n7oui4ybo0raUlcncQnxx+g+WgDF/vLoYDmoqSl/dJUnt7XRCoTZjij0Z6Pc2LiNS4EBBkNvoeOJXN+yPWWSZeANOhwJq/98nKVwNdoL8B5AROxBKBL0gjh8DMhdCh3eJnrA0yqhLpplwmyup6IajvAOIGfKGVx3VmCRGnOMpe5QAdG0bT8CAeeep0d6z6nqjSJnQiZWEllLMWrmz6k+fE9rGk8MVqYgsGv5ZH2i1Opr+9kajzB5d74hKQ+KS3d/WVMLhtgdu1lriRiOR/4nDVunaR24x7qp3UV5Cb/fJvC83nv26W81LIK58SYNFmwq4hsGx/5BwKlzYRma2NUthgOJSew4i7ru9nJYCQF5tApb2yvjiDQKJV/IfJKh0o6qssSLKv/jcAoRKHQQzE2Lj2OMV5OkWFc4MZIpsev8uXWXRx6ZicbGk8QZLxxgwe+x/rlR3h3816+f2E7lbEU+ZDn3vKVpePCdFovzCISHqbl5EIoQOteKMPB1rto65zNyfOz+KOrGl06lHPQyi/WOohH0/T0l1MZH6A3GUEKl7Pmr2la6wBrBWWRDP2DUcqjKVKBGom9RZmABAykwnglafpSJSPQvsfiOR0EQ7ExVmazA8cY6N4K1iw6RdAXRwi4mgrheT5Dvs4LeuS81a15Ll/3dQisFVSVpnj7sf1sX/sZvhAc+6UOrQyBVUQ8gx/orFmDsZqtaw/y1qZ9zKjp5vDpenyjcNe+cLNmTiUdf/bEOddVQ0VpgsOn54ET+EYxvWKALSu+5tGG76it7MNaiZKGQ23zCIcMfUMxBnrjN3fmHHvCAlp+vJcXWx6itqoXpAEnUNLx8iMfo5Xh1i17R3PJYCpC2cZ3qK3sQ8WGEDDuXlAQuFKGHzpmopXhTNfk0bmxs7uC1w6uJul79AxFkMIiBJy5UoUWjrZLU5DCFdTARDHuDqVw+OkSwI0MCEW4gtNF2BPrBCo8fKNbtILWX9aUDqFqHnn7AAAAAElFTkSuQmCC)](https://www.fnp.org.pl/en/)\n[![PL Funding](https://img.shields.io/static/v1?label=PL%20Funding%20by&color=d21132&message=NCN&logoWidth=25&logo=image/png;base64,iVBORw0KGgoAAAANSUhEUgAAABQAAAANCAYAAACpUE5eAAAABmJLR0QA/wD/AP+gvaeTAAAAKUlEQVQ4jWP8////fwYqAiZqGjZqIHUAy4dJS6lqIOMdEZvRZDPcDQQAb3cIaY1Sbi4AAAAASUVORK5CYII=)](https://www.ncn.gov.pl/?language=en)\n[![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)\n\n[![Github Actions Build Status](https://github.com/open-atmos/PyMPDATA/actions/workflows/tests+pypi.yml/badge.svg?branch=main)](https://github.com/open-atmos/PyMPDATA/actions)\n[![Appveyor Build status](http://ci.appveyor.com/api/projects/status/github/open-atmos/PyMPDATA?branch=main&svg=true)](https://ci.appveyor.com/project/slayoo/pympdata/branch/main)\n[![Coverage Status](https://codecov.io/gh/open-atmos/PyMPDATA/branch/main/graph/badge.svg)](https://app.codecov.io/gh/open-atmos/PyMPDATA)\n\n[![PyPI version](https://badge.fury.io/py/PyMPDATA.svg)](https://pypi.org/project/PyMPDATA)\n[![API docs](https://shields.mitmproxy.org/badge/docs-pdoc.dev-brightgreen.svg)](https://open-atmos.github.io/PyMPDATA/index.html)\n\nPyMPDATA is a high-performance Numba-accelerated Pythonic implementation of the MPDATA \n  algorithm of Smolarkiewicz et al. used in geophysical fluid dynamics and beyond.\nMPDATA numerically solves generalised transport equations -\n  partial differential equations used to model conservation/balance laws, scalar-transport problems,\n  convection-diffusion phenomena.\nAs of the current version, PyMPDATA supports homogeneous transport\n  in 1D, 2D and 3D using structured meshes, optionally\n  generalised by employment of a Jacobian of coordinate transformation. \nPyMPDATA includes implementation of a set of MPDATA variants including\n  the non-oscillatory option, infinite-gauge, divergent-flow, double-pass donor cell (DPDC) and \n  third-order-terms options.\nIt also features support for integration of Fickian-terms in advection-diffusion\n  problems using the pseudo-transport velocity approach.\nIn 2D and 3D simulations, domain-decomposition is used for multi-threaded parallelism.\n\nPyMPDATA is engineered purely in Python targeting both performance and usability,\n    the latter encompassing research users', developers' and maintainers' perspectives.\nFrom researcher's perspective, PyMPDATA offers hassle-free installation on multitude\n  of platforms including Linux, OSX and Windows, and eliminates compilation stage\n  from the perspective of the user.\nFrom developers' and maintainers' perspective, PyMPDATA offers a suite of unit tests, \n  multi-platform continuous integration setup,\n  seamless integration with Python development aids including debuggers and profilers.\n\nPyMPDATA design features\n  a custom-built multi-dimensional Arakawa-C grid layer allowing\n  to concisely represent multi-dimensional stencil operations on both\n  scalar and vector fields.\nThe grid layer is built on top of NumPy's ndarrays (using \"C\" ordering)\n  using the Numba's ``@njit`` functionality for high-performance array traversals.\nIt enables one to code once for multiple dimensions, and automatically\n  handles (and hides from the user) any halo-filling logic related with boundary conditions.\nNumba ``prange()`` functionality is used for implementing multi-threading \n  (it offers analogous functionality to OpenMP parallel loop execution directives).\nThe Numba's deviation from Python semantics rendering [closure variables\n  as compile-time constants](https://numba.pydata.org/numba-doc/dev/reference/pysemantics.html#global-and-closure-variables)\n  is extensively exploited within ``PyMPDATA``\n  code base enabling the just-in-time compilation to benefit from \n  information on domain extents, algorithm variant used and problem\n  characteristics (e.g., coordinate transformation used, or lack thereof).\n\nA separate project called [``PyMPDATA-MPI``](https://github.com/open-atmos/PyMPDATA-MPI) \n  depicts how [``numba-mpi``](https://pypi.org/project/numba-mpi) can be used\n  to enable distributed memory parallelism in PyMPDATA.\n\nThe [``PyMPDATA-examples``](https://pypi.org/project/PyMPDATA-examples/) \n  package covers a set of examples presented in the form of Jupyer notebooks\n  offering single-click deployment in the cloud using [mybinder.org](https://mybinder.org)\n  or using [colab.research.google.com](https://colab.research.google.com/).\nThe examples reproduce results from several published\n  works on MPDATA and its applications, and provide a validation of the implementation\n  and its performance.\n \n## Dependencies and installation\n\nTo install PyMPDATA, one may use: ``pip install PyMPDATA`` (or \n``pip install git+https://github.com/open-atmos/PyMPDATA.git`` to get updates beyond the latest release).\nPyMPDATA depends on ``NumPy`` and ``Numba``.\n\nRunning the tests shipped with the package requires additional packages listed in the \n[test-time-requirements.txt](https://github.com/open-atmos/PyMPDATA/blob/main/test-time-requirements.txt) file\n(which include ``PyMPDATA-examples``, see below).\n\n## Examples (Jupyter notebooks reproducing results from literature):\n\nPyMPDATA examples are bundled with PyMPDATA and located in the `examples` subfolder.\nThey constitute a separate [``PyMPDATA_examples``](https://pypi.org/p/PyMPDATA-examples) Python package which is also available at PyPI.\nThe examples have additional dependencies listed in [``PyMPDATA_examples`` package ``setup.py``](https://github.com/open-atmos/PyMPDATA/blob/main/examples/setup.py) file.\nRunning the examples requires the ``PyMPDATA_examples`` package to be installed.\nSince the examples package includes Jupyter notebooks (and their execution requires write access), the suggested install and launch steps are:\n```\ngit clone https://github.com/open-atmos/PyMPDATA-examples.git\ncd PyMPDATA-examples\npip install -e .\njupyter-notebook\n```\nAlternatively, one can also install the examples package from pypi.org by using ``pip install PyMPDATA-examples``.\n  \n## Package structure and API:\n\nIn short, PyMPDATA numerically solves the following equation:\n\n$$ \\partial_t (G \\psi) + \\nabla \\cdot (Gu \\psi) + \\mu \\Delta (G \\psi) = 0 $$\n\nwhere scalar field $\\psi$ is referred to as the advectee,\nvector field u is referred to as advector, and the G factor corresponds to optional coordinate transformation.\nThe inclusion of the Fickian diffusion term is optional and is realised through modification of the\nadvective velocity field with MPDATA handling both the advection and diffusion (for discussion\nsee, e.g. [Smolarkiewicz and Margolin 1998](https://doi.org/10.1006/jcph.1998.5901), sec. 3.5, par. 4).\n\nThe key classes constituting the PyMPDATA interface are summarised below with code\nsnippets exemplifying usage of PyMPDATA from Python, Julia and Matlab.\n\nA [pdoc-generated](https://pdoc.dev/) documentation of PyMPDATA public API is maintained at: [https://open-atmos.github.io/PyMPDATA](  https://open-atmos.github.io/PyMPDATA/index.html) \n\n#### Options class\n\nThe [``Options``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/options.html) class\ngroups both algorithm variant options as well as some implementation-related\nflags that need to be set at the first place. All are set at the time\nof instantiation using the following keyword arguments of the constructor \n(all having default values indicated below):\n- ``n_iters: int = 2``: number of iterations (2 means upwind + one corrective iteration)\n- ``infinite_gauge: bool = False``: flag enabling the infinite-gauge option (does not maintain sign of the advected field, thus in practice implies switching flux corrected transport on)\n- ``divergent_flow: bool = False``: flag enabling divergent-flow terms when calculating antidiffusive velocity\n- ``nonoscillatory: bool = False``: flag enabling the non-oscillatory or monotone variant (a.k.a flux-corrected transport option, FCT)\n- ``third_order_terms: bool = False``: flag enabling third-order terms\n- ``epsilon: float = 1e-15``: value added to potentially zero-valued denominators \n- ``non_zero_mu_coeff: bool = False``: flag indicating if code for handling the Fickian term is to be optimised out\n- ``DPDC: bool = False``: flag enabling double-pass donor cell option (recursive pseudovelocities)\n- ``dimensionally_split: bool = False``: flag disabling cross-dimensional terms in antidiffusive velocity\n- ``dtype: np.floating = np.float64``: floating point precision\n\nFor a discussion of the above options, see e.g., [Smolarkiewicz & Margolin 1998](https://doi.org/10.1006/jcph.1998.5901),\n[Jaruga, Arabas et al. 2015](https://doi.org/10.5194/gmd-8-1005-2015) and [Olesik, Arabas et al. 2020](https://arxiv.org/abs/2011.14726)\n(the last with examples using PyMPDATA).\n\nIn most use cases of PyMPDATA, the first thing to do is to instantiate the \n[``Options``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/options.html) class \nwith arguments suiting the problem at hand, e.g.:\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nusing Pkg\nPkg.add(\"PyCall\")\nusing PyCall\nOptions = pyimport(\"PyMPDATA\").Options\noptions = Options(n_iters=2)\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nOptions = py.importlib.import_module('PyMPDATA').Options;\noptions = Options(pyargs('n_iters', 2));\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nfrom PyMPDATA import Options\noptions = Options(n_iters=2)\n```\n</details>\n\n#### Arakawa-C grid layer and boundary conditions\n\nIn PyMPDATA, the solution domain is assumed to extend from the\nfirst cell's boundary to the last cell's boundary (thus the\nfirst scalar field value is at $\\[\\Delta x/2, \\Delta y/2\\]$.\nThe [``ScalarField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/scalar_field.html)\nand [``VectorField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/vector_field.html) classes implement the\n[Arakawa-C staggered grid](https://en.wikipedia.org/wiki/Arakawa_grids#Arakawa_C-grid) logic\nin which:\n- scalar fields are discretised onto cell centres (one value per cell),\n- vector field components are discretised onto cell walls.\n\nThe schematic of the employed grid/domain layout in two dimensions is given below\n(with the Python code snippet generating the figure as a part of CI workflow):\n<details>\n<summary>Python code (click to expand)</summary>\n\n```Python\nimport numpy as np\nfrom matplotlib import pyplot\n\ndx, dy = .2, .3\ngrid = (10, 5)\n\npyplot.scatter(*np.mgrid[\n        dx / 2 : grid[0] * dx : dx, \n        dy / 2 : grid[1] * dy : dy\n    ], color='red', \n    label='scalar-field values at cell centres'\n)\npyplot.quiver(*np.mgrid[\n        0 : (grid[0]+1) * dx : dx, \n        dy / 2 : grid[1] * dy : dy\n    ], 1, 0, pivot='mid', color='green', width=.005,\n    label='vector-field x-component values at cell walls'\n)\npyplot.quiver(*np.mgrid[\n        dx / 2 : grid[0] * dx : dx,\n        0: (grid[1] + 1) * dy : dy\n    ], 0, 1, pivot='mid', color='blue', width=.005,\n    label='vector-field y-component values at cell walls'\n)\npyplot.xticks(np.linspace(0, grid[0]*dx, grid[0]+1))\npyplot.yticks(np.linspace(0, grid[1]*dy, grid[1]+1))\npyplot.title(f'staggered grid layout (grid={grid}, dx={dx}, dy={dy})')\npyplot.xlabel('x')\npyplot.ylabel('y')\npyplot.legend(bbox_to_anchor=(.1, -.1), loc='upper left', ncol=1)\npyplot.grid()\npyplot.savefig('readme_grid.png')\n```\n</details>\n\n![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_grid.png)\n\nThe ``__init__`` methods of\n[``ScalarField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/scalar_field.html)\nand [``VectorField``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/vector_field.html)\nhave the following signatures:\n- [``ScalarField(data: np.ndarray, halo: int, boundary_conditions)``](https://github.com/open-atmos/PyMPDATA/blob/main/PyMPDATA/scalar_field.py)\n- [``VectorField(data: Tuple[np.ndarray, ...], halo: int, boundary_conditions)``](https://github.com/open-atmos/PyMPDATA/blob/main/PyMPDATA/vector_field.py)\nThe ``data`` parameters are expected to be Numpy arrays or tuples of Numpy arrays, respectively.\nThe ``halo`` parameter is the extent of ghost-cell region that will surround the\ndata and will be used to implement boundary conditions. \nIts value (in practice 1 or 2) is\ndependent on maximal stencil extent for the MPDATA variant used and\ncan be easily obtained using the ``Options.n_halo`` property.\n\nAs an example, the code below shows how to instantiate a scalar\nand a vector field given a 2D constant-velocity problem,\nusing a grid of 24x24 points, Courant numbers of -0.5 and -0.25\nin \"x\" and \"y\" directions, respectively, with periodic boundary \nconditions and with an initial Gaussian signal in the scalar field\n(settings as in Fig.&nbsp;5 in [Arabas et al. 2014](https://doi.org/10.3233/SPR-140379)):\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nScalarField = pyimport(\"PyMPDATA\").ScalarField\nVectorField = pyimport(\"PyMPDATA\").VectorField\nPeriodic = pyimport(\"PyMPDATA.boundary_conditions\").Periodic\n\nnx, ny = 24, 24\nCx, Cy = -.5, -.25\nidx = CartesianIndices((nx, ny))\nhalo = options.n_halo\nadvectee = ScalarField(\n    data=exp.(\n        -(getindex.(idx, 1) .- .5 .- nx/2).^2 / (2*(nx/10)^2) \n        -(getindex.(idx, 2) .- .5 .- ny/2).^2 / (2*(ny/10)^2)\n    ),  \n    halo=halo, \n    boundary_conditions=(Periodic(), Periodic())\n)\nadvector = VectorField(\n    data=(fill(Cx, (nx+1, ny)), fill(Cy, (nx, ny+1))),\n    halo=halo,\n    boundary_conditions=(Periodic(), Periodic())    \n)\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nScalarField = py.importlib.import_module('PyMPDATA').ScalarField;\nVectorField = py.importlib.import_module('PyMPDATA').VectorField;\nPeriodic = py.importlib.import_module('PyMPDATA.boundary_conditions').Periodic;\n\nnx = int32(24);\nny = int32(24);\n  \nCx = -.5;\nCy = -.25;\n\n[xi, yi] = meshgrid(double(0:1:nx-1), double(0:1:ny-1));\n\nhalo = options.n_halo;\nadvectee = ScalarField(pyargs(...\n    'data', py.numpy.array(exp( ...\n        -(xi+.5-double(nx)/2).^2 / (2*(double(nx)/10)^2) ...\n        -(yi+.5-double(ny)/2).^2 / (2*(double(ny)/10)^2) ...\n    )), ... \n    'halo', halo, ...\n    'boundary_conditions', py.tuple({Periodic(), Periodic()}) ...\n));\nadvector = VectorField(pyargs(...\n    'data', py.tuple({ ...\n        Cx * py.numpy.ones(int32([nx+1 ny])), ... \n        Cy * py.numpy.ones(int32([nx ny+1])) ...\n     }), ...\n    'halo', halo, ...\n    'boundary_conditions', py.tuple({Periodic(), Periodic()}) ...\n));\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nfrom PyMPDATA import ScalarField\nfrom PyMPDATA import VectorField\nfrom PyMPDATA.boundary_conditions import Periodic\nimport numpy as np\n\nnx, ny = 24, 24\nCx, Cy = -.5, -.25\nhalo = options.n_halo\n\nxi, yi = np.indices((nx, ny), dtype=float)\nadvectee = ScalarField(\n  data=np.exp(\n    -(xi+.5-nx/2)**2 / (2*(nx/10)**2)\n    -(yi+.5-ny/2)**2 / (2*(ny/10)**2)\n  ),\n  halo=halo,\n  boundary_conditions=(Periodic(), Periodic())\n)\nadvector = VectorField(\n  data=(np.full((nx + 1, ny), Cx), np.full((nx, ny + 1), Cy)),\n  halo=halo,\n  boundary_conditions=(Periodic(), Periodic())\n)\n```\n</details>\n\nNote that the shapes of arrays representing components \nof the velocity field are different than the shape of\nthe scalar field array due to employment of the staggered grid.\n\nBesides the exemplified [``Periodic``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/periodic.html) class representing \nperiodic boundary conditions, PyMPDATA supports \n[``Extrapolated``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/extrapolated.html), \n[``Constant``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/constant.html) and\n[``Polar``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/boundary_conditions/polar.html) \nboundary conditions.\n\n#### Stepper\n\nThe logic of the MPDATA iterative solver is represented\nin PyMPDATA by the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) class.\n\nWhen instantiating the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html), the user has a choice \nof either supplying just the  number of dimensions or specialising the stepper for a given grid:\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nStepper = pyimport(\"PyMPDATA\").Stepper\n\nstepper = Stepper(options=options, n_dims=2)\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nStepper = py.importlib.import_module('PyMPDATA').Stepper;\n\nstepper = Stepper(pyargs(...\n  'options', options, ...\n  'n_dims', int32(2) ...\n));\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nfrom PyMPDATA import Stepper\n\nstepper = Stepper(options=options, n_dims=2)\n```\n</details>\nor\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nstepper = Stepper(options=options, grid=(nx, ny))\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nstepper = Stepper(pyargs(...\n  'options', options, ...\n  'grid', py.tuple({nx, ny}) ...\n));\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nstepper = Stepper(options=options, grid=(nx, ny))\n```\n</details>\n\nIn the latter case, noticeably \nfaster execution can be expected, however the resultant\nstepper is less versatile as bound to the given grid size.\nIf number of dimensions is supplied only, the integration\nmight take longer, yet same instance of the\nstepper can be used for different grids.  \n\nSince creating an instance of the [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) class\ninvolves time-consuming compilation of the algorithm code,\nthe class is equipped with a cache logic - subsequent\ncalls with same arguments return references to previously\ninstantiated objects. Instances of [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) contain no\nmutable data and are (thread-)safe to be reused.\n\nThe init method of [``Stepper``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/stepper.html) has an optional\n``non_unit_g_factor`` argument which is a Boolean flag \nenabling handling of the G factor term which can be used to \nrepresent coordinate transformations and/or variable fluid density. \n\nOptionally, the number of threads to use for domain decomposition\nin the first (non-contiguous) dimension during 2D and 3D calculations\nmay be specified using the optional ``n_threads`` argument with a\ndefault value of ``numba.get_num_threads()``. The multi-threaded\nlogic of PyMPDATA depends thus on settings of numba, namely on the\nselected threading layer (either via ``NUMBA_THREADING_LAYER`` env \nvar or via ``numba.config.THREADING_LAYER``) and the selected size of the \nthread pool (``NUMBA_NUM_THREADS`` env var or ``numba.config.NUMBA_NUM_THREADS``).\n\n\n#### Solver\n\nInstances of the [``Solver``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html) class are used to control\nthe integration and access solution data. During instantiation, \nadditional memory required by the solver is \nallocated according to the options provided. \n\nThe only method of the [``Solver``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html) class besides the\ninit is [``advance(n_steps, mu_coeff, ...)``](https://open-atmos.github.io/PyMPDATA/PyMPDATA/solver.html#Solver.advance) \nwhich advances the solution by ``n_steps`` timesteps, optionally\ntaking into account a given diffusion coefficient ``mu_coeff``.\n\nSolution state is accessible through the ``Solver.advectee`` property.\nMultiple solver[s] can share a single stepper, e.g., as exemplified in the shallow-water system solution in the examples package.\n\nContinuing with the above code snippets, instantiating\na solver and making 75 integration steps looks as follows:\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nSolver = pyimport(\"PyMPDATA\").Solver\nsolver = Solver(stepper=stepper, advectee=advectee, advector=advector)\nsolver.advance(n_steps=75)\nstate = solver.advectee.get()\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nSolver = py.importlib.import_module('PyMPDATA').Solver;\nsolver = Solver(pyargs('stepper', stepper, 'advectee', advectee, 'advector', advector));\nsolver.advance(pyargs('n_steps', 75));\nstate = solver.advectee.get();\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nfrom PyMPDATA import Solver\n\nsolver = Solver(stepper=stepper, advectee=advectee, advector=advector)\nstate_0 = solver.advectee.get().copy()\nsolver.advance(n_steps=75)\nstate = solver.advectee.get()\n```\n</details>\n\nNow let's plot the results using `matplotlib` roughly as in Fig.&nbsp;5 in [Arabas et al. 2014](https://doi.org/10.3233/SPR-140379):\n\n<details>\n<summary>Python code (click to expand)</summary>\n\n```Python\ndef plot(psi, zlim, norm=None):\n    xi, yi = np.indices(psi.shape)\n    fig, ax = pyplot.subplots(subplot_kw={\"projection\": \"3d\"})\n    pyplot.gca().plot_wireframe(\n        xi+.5, yi+.5, \n        psi, color='red', linewidth=.5\n    )\n    ax.set_zlim(zlim)\n    for axis in (ax.xaxis, ax.yaxis, ax.zaxis):\n        axis.pane.fill = False\n        axis.pane.set_edgecolor('black')\n        axis.pane.set_alpha(1)\n    ax.grid(False)\n    ax.set_zticks([])\n    ax.set_xlabel('x/dx')\n    ax.set_ylabel('y/dy')\n    ax.set_proj_type('ortho') \n    cnt = ax.contourf(xi+.5, yi+.5, psi, zdir='z', offset=-1, norm=norm)\n    cbar = pyplot.colorbar(cnt, pad=.1, aspect=10, fraction=.04)\n    return cbar.norm\n\nzlim = (-1, 1)\nnorm = plot(state_0, zlim)\npyplot.savefig('readme_gauss_0.png')\nplot(state, zlim, norm)\npyplot.savefig('readme_gauss.png')\n```\n</details>\n\n![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_gauss_0.png)    \n![plot](https://github.com/open-atmos/PyMPDATA/releases/download/tip/readme_gauss.png)\n\n#### Debugging\n\nPyMPDATA relies heavily on Numba to provide high-performance \nnumber crunching operations. Arguably, one of the key advantage \nof embracing Numba is that it can be easily switched off. This\nbrings multiple-order-of-magnitude drop in performance, yet \nit also make the entire code of the library susceptible to\ninteractive debugging, one way of enabling it is by setting the \nfollowing environment variable before importing PyMPDATA:\n<details>\n<summary>Julia code (click to expand)</summary>\n\n```Julia\nENV[\"NUMBA_DISABLE_JIT\"] = \"1\"\n```\n</details>\n<details>\n<summary>Matlab code (click to expand)</summary>\n\n```Matlab\nsetenv('NUMBA_DISABLE_JIT', '1');\n```\n</details>\n<details open>\n<summary>Python code (click to expand)</summary>\n\n```Python\nimport os\nos.environ[\"NUMBA_DISABLE_JIT\"] = \"1\"\n```\n</details>\n\n## Contributing, reporting issues, seeking support \n\nSubmitting new code to the project, please preferably use [GitHub pull requests](https://github.com/open-atmos/PyMPDATA/pulls) \n(or the [PyMPDATA-examples PR site](https://github.com/open-atmos/PyMPDATA-examples/pulls) if working on examples) - it helps to keep record of code authorship, \ntrack and archive the code review workflow and allows to benefit\nfrom the continuous integration setup which automates execution of tests \nwith the newly added code. \n\nAs of now, the copyright to the entire PyMPDATA codebase is with the Jagiellonian\nUniversity, and code contributions are assumed to imply transfer of copyright.\nShould there be a need to make an exception, please indicate it when creating\na pull request or contributing code in any other way. In any case, \nthe license of the contributed code must be compatible with GPL v3.\n\nDeveloping the code, we follow [The Way of Python](https://www.python.org/dev/peps/pep-0020/) and \nthe [KISS principle](https://en.wikipedia.org/wiki/KISS_principle).\nThe codebase has greatly benefited from [PyCharm code inspections](https://www.jetbrains.com/help/pycharm/code-inspection.html)\nand [Pylint](https://pylint.org) code analysis (Pylint checks are part of the\nCI workflows).\n\nIssues regarding any incorrect, unintuitive or undocumented bahaviour of\nPyMPDATA are best to be reported on the [GitHub issue tracker](https://github.com/open-atmos/PyMPDATA/issues/new).\nFeature requests are recorded in the \"Ideas...\" [PyMPDATA wiki page](https://github.com/open-atmos/PyMPDATA/wiki/Ideas-for-new-features-and-examples).\n\nWe encourage to use the [GitHub Discussions](https://github.com/open-atmos/PyMPDATA/discussions) feature\n(rather than the issue tracker) for seeking support in understanding, using and extending PyMPDATA code.\n\nPlease use the PyMPDATA issue-tracking and dicsussion infrastructure for `PyMPDATA-examples` as well.\nWe look forward to your contributions and feedback.\n\n## Credits:\nDevelopment of PyMPDATA was supported by the EU through a grant of the [Foundation for Polish Science](http://fnp.org.pl) (POIR.04.04.00-00-5E1C/18).\n\ncopyright: Jagiellonian University   \nlicence: GPL v3   \n\n## Other open-source MPDATA implementations:\n- mpdat_2d in babyEULAG (FORTRAN)\n  https://github.com/igfuw/bE_SDs/blob/master/babyEULAG.SDs.for#L741\n- mpdata-oop (C++, Fortran, Python)\n  https://github.com/igfuw/mpdata-oop\n- apc-llc/mpdata (C++)\n  https://github.com/apc-llc/mpdata\n- libmpdata++ (C++):\n  https://github.com/igfuw/libmpdataxx\n- AtmosFOAM:\n  https://github.com/AtmosFOAM/AtmosFOAM/tree/947b192f69d973ea4a7cfab077eb5c6c6fa8b0cf/applications/solvers/advection/MPDATAadvectionFoam\n\n## Other Python packages for solving hyperbolic transport equations\n\n- PyPDE: https://pypi.org/project/PyPDE/\n- FiPy: https://pypi.org/project/FiPy/\n- ader: https://pypi.org/project/ader/\n- centpy: https://pypi.org/project/centpy/\n- mattflow: https://pypi.org/project/mattflow/\n- FastFD: https://pypi.org/project/FastFD/\n- Pyclaw: https://www.clawpack.org/pyclaw/\n",
    "bugtrack_url": null,
    "license": "GPL-3.0",
    "summary": "Numba-accelerated Pythonic implementation of MPDATA with examples in Python, Julia and Matlab",
    "version": "1.1.6",
    "project_urls": {
        "Documentation": "https://open-atmos.github.io/PyMPDATA",
        "Source": "https://github.com/open-atmos/PyMPDATA",
        "Tracker": "https://github.com/open-atmos/PyMPDATA/issues"
    },
    "split_keywords": [
        "atmospheric-modelling",
        " numba",
        " numerical-integration",
        " advection",
        " pde-solver",
        " advection-diffusion"
    ],
    "urls": [
        {
            "comment_text": "",
            "digests": {
                "blake2b_256": "59cefa0b6ad634acd3c3d76d8f61a6ff07d6d78881c4dd9e382ee2c4dda5dece",
                "md5": "3c012a8531e5d7d5e660e5ad4710c339",
                "sha256": "0b31af52ffed8547a072f4ca41474fac41d79e87c9e4902efe7293f22910adab"
            },
            "downloads": -1,
            "filename": "PyMPDATA-1.1.6-py3-none-any.whl",
            "has_sig": false,
            "md5_digest": "3c012a8531e5d7d5e660e5ad4710c339",
            "packagetype": "bdist_wheel",
            "python_version": "py3",
            "requires_python": null,
            "size": 58375,
            "upload_time": "2024-09-21T15:03:35",
            "upload_time_iso_8601": "2024-09-21T15:03:35.480484Z",
            "url": "https://files.pythonhosted.org/packages/59/ce/fa0b6ad634acd3c3d76d8f61a6ff07d6d78881c4dd9e382ee2c4dda5dece/PyMPDATA-1.1.6-py3-none-any.whl",
            "yanked": false,
            "yanked_reason": null
        },
        {
            "comment_text": "",
            "digests": {
                "blake2b_256": "b1936a989073744297766a6c32ac3f8751450378b39abcb235429af2f4521f98",
                "md5": "4febfa33188b282437f68bccd3a5839d",
                "sha256": "a75cf38f7616b65d119bbba9c45b9817d4096c6da9f28eae00ab4136e4e84709"
            },
            "downloads": -1,
            "filename": "pympdata-1.1.6.tar.gz",
            "has_sig": false,
            "md5_digest": "4febfa33188b282437f68bccd3a5839d",
            "packagetype": "sdist",
            "python_version": "source",
            "requires_python": null,
            "size": 220308,
            "upload_time": "2024-09-21T15:03:38",
            "upload_time_iso_8601": "2024-09-21T15:03:38.436500Z",
            "url": "https://files.pythonhosted.org/packages/b1/93/6a989073744297766a6c32ac3f8751450378b39abcb235429af2f4521f98/pympdata-1.1.6.tar.gz",
            "yanked": false,
            "yanked_reason": null
        }
    ],
    "upload_time": "2024-09-21 15:03:38",
    "github": true,
    "gitlab": false,
    "bitbucket": false,
    "codeberg": false,
    "github_user": "open-atmos",
    "github_project": "PyMPDATA",
    "travis_ci": false,
    "coveralls": false,
    "github_actions": true,
    "appveyor": true,
    "lcname": "pympdata"
}
        
Elapsed time: 0.35063s