Getting started with josiepy

A Python framework to simulate Partial Differential Equations focused on usability

(without sacrifying performances)

Ruben Di Battista


M. Massot Research Team

Mathematical modeling

  • Mathematical modeling of two-phase injection in rocket engines

  • Mathematical modeling of plasma flows for heliophysics and hall thrusters

Numerical Calculus for High Performance Computing (HPC)

  • Development of numerical schemes and algorithms adapted to HPC


A C++ wrapper on p4est, a quad/octo-trees for parallel AMR computation, to simulate PDEs w/ AMR

  • C++ is notoriously hard to master

  • p4est is an advanced library w/ features not immediately needed when exploring numerical schemes


An easy to use playground to easily implement numerical schemes on 1D, 2D (and maybe 3D in the future) structured (and maybe unstructured and non-conformal meshes… dreaming is still free 😇) meshes, in Python

Design Choices

  • Use Python instead of a Domain Specific Language (DSL), like the one OpenFOAM uses, or another scripting language, like lua in canoP to describe your case configuration. “You code your case configuration”

  • Modern Python (Python >= 3.7) features allows static checking (through the usage of type hinting) of the code that predates some advantages of the compiled languages

  • The numpy library and its API allows to write algorithms once while being able to run them on different architectures and different programming paradigms

    • CPU Shared memory (numpy + OpenMP-enabled BLAS library)

    • CPU distributed memory (Dask)

    • Nvidia GPUs with CUDA (cupy)

    • (Soon) AMD GPUs with HIP (cupy)

    • … FPGA? Anyone? It’s a cool project too…

    Everything mostly at the same computational speed as C and C++, or faster (GPU for example)

Numpy for josiepy 101

Vectorized operations

What you would generally do in a C++-like compiled code is:

for( cell : mesh.cells() ) {

What we generally will do in numpy is operating on an array containing all the cells values at the same time.

Let’s imagine we have the velocity vector stored in an array of size N_\text{cells} \times \text{dim}$, where $\text{dim} is the dimensionality of the problem (e.g. 2D)

cell_values = np.random.random((30, 2))
U = cell_values[:, 0]
V = cell_values[:, 1]

Let’s imagine now to have to compute a normal product between all the velocity values in each cell, and a singular normal vector

normal = np.array([1, 0])

What we could do is repeating the normal array N_\text{cells} times and multipliying it by cell_values

… but there’s a better way! The operation we want to do is

\left\{\mathbf{U}\right\} \cdot \mathbf{n} = U_{ij}n_j

normal_velocities = np.einsum("ij,j->i", cell_values, normal)

np.array_equal(U, normal_velocities)


numpy.einsum is an implementation of the Einstein’s tensor summation on steroids 💪. You specify the input vectors dimensions and the output dimension you want to achieve, and it figures it out automatically.


U_{ij} n_j \rightarrow u_i

np.einsum("ij,j->i", cell_values, normal)
  • Repeated indices on the inputs are summed as a result if they’re not present in the result (after the ->): cell_values[i][0] * normal[0] + cell_values[i][1] * normal[1]

np.einsum avoids memory copies and if we provide the optimize=True keyword argument

A = np.random.random((1000, 1000))
B = np.random.random((1000, 1000))

%timeit np.einsum("ij,ij->", A, B)
799 µs ± 4.62 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
A = np.random.random((1000, 1000))
B = np.random.random((1000, 1000))

%timeit np.einsum("ij,ij->", A, B, optimize=True)
531 µs ± 43 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

it is gonna exploit parallelization of operations like and np.tensordot that are parallelized natively with OpenMP if the BLAS library which numpy is compiled on allows it

(hint: the classical pip install numpy pre-compiled wheel, does not)

Vectorized conditionals

Sometimes, for example to initialize the domain at the beginning of a simulation, or for an upwind scheme, you need to apply some conditionals of the type:

for (cell : mesh.cells()) {
    x, y = cell.centroid()
    if(x > 0.5) 

What you would do in vectorized form, is:

values = np.zeros(30)
x = np.linspace(0, 1, np.size(values))
[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0.]
values[np.where(x > 0.5)] = 1

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1.
 1. 1. 1. 1. 1. 1.]

Basics of josiepy

Problem statement

\partial_t \mathbf{q} + 
    \nabla \cdot \left(\underline{\underline{\mathbf{F}}}(\mathbf{q}) 
        + \underline{\underline{\underline{\underline{D(\mathbf{q})}}}} \cdot \nabla \mathbf{q} \right) +
    \underline{\underline{\underline{\mathbf{B}(\mathbf{q})}}} \cdot \nabla\mathbf{q} = 

\partial_t q_p + 
    \partial_{x_r} \left(F_{pr}(q_p) + D_{pqrs}(q_p)\partial_{x_s}q_q\right) + 
    B_{pqr}\partial_{x_r}q_q = 

The State object

\partial_t \mathbf{q} + 
    \nabla \cdot \left(\underline{\underline{\mathbf{F}}}(\mathbf{q}) 
        + \underline{\underline{\underline{\underline{D(\mathbf{q})}}}} \cdot \nabla \mathbf{q} \right) +
    \underline{\underline{\underline{\mathbf{B}(\mathbf{q})}}} \cdot \nabla\mathbf{q} = 

Represents the \mathbf{q}

from josie.state import State
from josie.fluid.fields import FluidFields
class EulerFields(FluidFields):
    rho = 0
    rhoU = 1
    rhoV = 2
    rhoE = 3
    U = 4
    V = 5
    rhoe = 6
    p = 7
    c = 8
class EulerState(State):
    fields = EulerFields
rnd_state = np.random.random(len(EulerFields)).view(EulerState)
[0.78036773 0.8141633  0.48401207 0.4379539  0.55518386 0.39926222
 0.59759767 0.93960397 0.64511782]
U1 = rnd_state[EulerFields.U]
U2 = rnd_state[4]
U3 = rnd_state[-5]
U4 = rnd_state[rnd_state.fields.U]

assert np.array_equal(U1, U2) and np.array_equal(U2, U3) and np.array_equal(U3, U4)

State can be multidimensional! (e.g. a 2D 100x100 mesh, each cell containing an Euler state)

rnd_state = np.random.random((100, 100, len(EulerFields)))
U = rnd_state[..., EulerFields.U]
(100, 100)

The Ellipsis object (...) allows to index on the last axis, whatever dimension the array is

rnd_state = np.random.random((100, 100, 100, len(EulerFields)))
U = rnd_state[..., EulerFields.U]
(100, 100, 100)

Let’s compute the normal velocity again…

UV = rnd_state[..., EulerFields.U : EulerFields.V + 1]

(100, 100, 100, 2)
U_norm = np.einsum("...ij,j->...i", UV, normal)

assert np.array_equal(U_norm, U)
(100, 100, 100)

The Mesh object

A mesh object is used to create an manage a mesh. A mesh is built from BoundaryCurves:

from josie.boundary import Line, CircleArc

left = Line([0, 0], [0, 1])
bottom = CircleArc([0, 0], [1, 0], [0.5, 0.5])
right = Line([1, 0], [1, 1])
top = Line([0, 1], [1, 1])

for curve in [left, bottom, right, top]:


You can (and need to) attach BoundaryConditions to each BoundaryCurve of a domain

from josie.bc import Dirichlet, Neumann

Q_zero = np.zeros(len(EulerState.fields)).view(EulerState)
dQ_zero = Q_zero

left.bc = Dirichlet(Q_zero)
top.bc = Neumann(dQ_zero)
bottom.bc = Dirichlet(Q_zero)
right.bc = Dirichlet(Q_zero)

You can also provide mixed Boundary conditions per each field

from josie.bc import BoundaryCondition

Q_bc = EulerState(rho=Dirichlet(1), rhoU=Neumann(0), rhoV=Dirichlet(1), 
                  rhoE=Dirichlet(1), rhoe=Dirichlet(3), U=Neumann(0), 
                  V=Neumann(0), p=Dirichlet(1), c=Dirichlet(10), e=Dirichlet(1))

left.bc = BoundaryCondition(Q_bc)

Then you can effectively generate the mesh object

from josie.mesh import Mesh
from josie.mesh.cell import SimpleCell

mesh = Mesh(left, bottom, right, top, SimpleCell)
points = mesh.interpolate(20, 20)  # Gen. internal points
mesh.generate() # Generate connectivity


PS: You can create 1D meshes, creating a rectangular mesh with just 1 cell in the vertical direction

The Problem object.

\partial_t \mathbf{q} + 
    \nabla \cdot \left(\underline{\underline{\mathbf{F}}}(\mathbf{q}) 
        + \underline{\underline{\underline{\underline{D(\mathbf{q})}}}} \cdot \nabla \mathbf{q} \right) +
    \underline{\underline{\underline{\mathbf{B}(\mathbf{q})}}} \cdot \nabla\mathbf{q} = 

Problem is the “continuous” representation of our problem. It implements the \mathbf{F}, \mathbf{D}, \mathbf{B}, \mathbf{s} operators

from josie.problem import Problem

from josie.euler.problem import EulerProblem


The Scheme object

It represents the discrete representation of our problem, after being discretized with Finite-Volume-Method

V_i\partial_t \mathbf{q}_i &+
    \sum_{f \in \text{faces}} \left| \mathbf{F} \cdot \hat{\mathbf{n}} \right|_f S_f + 
    \sum_{f \in \text{faces}} \left| \left(\mathbf{D}(\mathbf{q})\cdot\nabla\mathbf{q} \right) \cdot \hat{\mathbf{n}} \right|_f S_f + \\
    &+ \left<\mathbf{B}(\mathbf{q})\right>_i \cdot \sum_{f \in \text{faces}} \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f
    - \left< \mathbf{s}(\mathbf{q})\right>_i V_i = \mathbf{0}

Schemes are organized in spatial and time schemes.

Spatial Schemes

Spatial schemes implement the different discretizations for each operators w.r.t. a single cell and its neighbour, whatever that neighbour is (no knowledge of left/right, top/bottom, etc…)


and are organized in specific classes:

from josie.scheme.nonconservative import NonConservativeScheme
from josie.scheme.convective import ConvectiveScheme
from josie.scheme.source import SourceScheme
# from josie.scheme.diffusion import DiffusionScheme <to be implemented>



It implements, in its ConvectiveScheme.F method, the discretization:

\mathbf{F}_\frac{1}{2}(\mathbf{q})= \left| \mathbf{F} \cdot \hat{\mathbf{n}} \right|_f S_f

on one face


Example Rusanov for Euler system


It implements, in its Rusanov.F method, the discretization:

            \left(\mathbf{F}(\mathbf{q}_i) + \mathbf{F}(\mathbf{q}_{i+1})\right) \cdot \hat{\mathbf{n}} - 
            \left|\text{max}(\lambda_p)\right|\left( \mathbf{q}_{i+1} - \mathbf{q}_i\right)

on one face

from josie.euler.schemes import Rusanov




A discretization of the type:

\left<\mathbf{B}(\mathbf{q})\right>_i \cdot \sum_{f \in \text{faces}} \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f

is assumed, and the NonConservativeScheme.G implements the term:

\mathbf{G}_{\frac{1}{2}}(\mathbf{q}) = \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f

Example Upwind for Baer-Nunziato non-conservative term


It implements, in its Upwind.G method, the discretization:

\mathbf{\alpha}_\frac{1}{2}(\mathbf{u}_i) \hat{\mathbf{n}} = 
        \left(\mathbf{u}_{i_f} \cdot \hat{\mathbf{n}} > 0\right) \alpha_i + 
        \left(\mathbf{u}_{i_f} \cdot \hat{\mathbf{n}} < 0\right) \alpha_{i+1}
    \right] \hat{\mathbf{n}}

on one face

from import Upwind


Time schemes

Time schemes simply implement the time update for the discretized equation:

V_i\partial_t \mathbf{q}_i &+
    \sum_{f \in \text{faces}} \left| \mathbf{F} \cdot \hat{\mathbf{n}} \right|_f S_f + 
    \sum_{f \in \text{faces}} \left| \left(\mathbf{D}(\mathbf{q})\cdot\nabla\mathbf{q} \right) \cdot \hat{\mathbf{n}} \right|_f S_f + \\
    &+ \left<\mathbf{B}(\mathbf{q})\right>_i \cdot \sum_{f \in \text{faces}} \left| \mathbf{q}\hat{\mathbf{n}} \right|_f S_f
    - \left< \mathbf{s}(\mathbf{q})\right>_i V_i = \mathbf{0}

In compact form:

\partial_t \mathbf{q}_i = - \frac{1}{V_i}\mathbf{f}\left(\mathbf{q}_i, t\right)

where \mathbf{f} contains all the discretized space operators. Integrating in time:

q_i^{k+1} = q_i^k - \frac{1}{V_i}\int_t^{t + \Delta t} \mathbf{f}\left(\mathbf{q}_i, \tilde{t}\right)\; \text{d}t


It implements the term \frac{1}{V_i}\int_t^{t+\Delta t} \mathbf{f}\left(\mathbf{q}_i, \tilde{t}\right)\; \text{d}t

from josie.scheme.time import TimeScheme


Example ExplicitEuler

\frac{1}{V_i} \int_t^{t + \Delta t} \mathbf{f}\left(\mathbf{q}_i, \tilde{t}\right)\; \text{d}t \sim
    \frac{1}{V_i}\mathbf{f}\left(\mathbf{q}_i^k, t_k\right) \Delta t

from josie.general.schemes.time import ExplicitEuler