Scikit-fdiff in short

Scikit-fdiff is a python library that aim to solve partial derivative equations without pain. As its name says, it uses finite difference method to discretize the spatial derivative. As its name does not say, it is based on *method of lines* where all the dimension of the PDE but the last (the time) is discretized. That turns the PDE in a high-dimension ODE that can be solved with standard numerical integration of ODEs.

The discretization is done in a symbolic way using sympy, and the exact Jacobian matrix associated with the resulting ODE is also obtained via symbolic derivation.

As final step, the symbolic description of the ODE is transformed in numerical routines thanks to a backend (for now, numpy as well as numba are available), and these routine can feed a ODE solver. The standard are available via scipy, and more specialized ones, that can take advantage of the fast computation of the Jacobian matrix have been written.

A simulation handler is also available, with some useful goodies as real-time display and on-disk persistence of the simulation results, handy interface for post-processing and more.

Should I use scikit-fdiff?

Scikit-fdiff is a robust way to use finite-difference in order to solve system of PDE. It should be able to deal with almost every model you will provide, but never in an optimal way. That make this tool perfect for model prototyping, where you want quick answer to your question. But, for a well established model, you will certainly find a better way to solve it using a specialized scheme.

Be aware that, using finite difference method, you will have all its drawback. Especially, you will be limited to simple geometry (rectangular geometry). The only way to solve problem in a more complex domain is to use variable change and penalisation method, and neither of them are trivial.

On the bright side, scikit-fdiff gives you access to high-order implicit solvers with adaptive time-stepping. It also allow to control the spatial scheme derivative per derivative, letting you play with the numerical solving of your problem, all of that without hard dependencies to install.

In short, scikit-fdiff is a great tool for yet unknown family of problem, when the solution are stiff, when you want to quickly access to a (non-optimal but accurate) solution as first approach, as long as your domain is simple enough.

Rational and software history

That tool was first written after we notice that no robust tool existed to solve a problem via finite-difference through the full modelling and solve process.

A lot of them exists for finite volume (see Clawpack or FiPy) and finite element method (see SfePy or the FEniCS Project) thanks to the increasing interest for them in the engineering field. But these method do not suit all the problem, or not without cumbersome algebraic work on the model.

As alternative, some interesting tools appear using pseudo spectral methods, as the Dedalus project. But no serious challenger using finite difference method.

A lot of different people write finite difference routine in their hand-made code, writing the same schemes again and again. This take a lot of time and can lead to error that are hard to catch.

For our internal use first, then as full library, we began to write a tool named Triflow. As this tool evolved and loose its first goal (solve TRansient Instable FLOW) and become a general PDE solver, we decided to change its name and join the scikit family.

Installation

The software is available on PyPI and conda on the scikit-fdiff channel. It should migrate soon in the conda forge.

To install scikit-fdiff, use

pip install scikit-fdiff

You can add some of the extra dependencies with

pip install scikit-fdiff[interactive,numba]

for jupyter interactive tools as well as high-performance numba back-end. For that later (and it’s true for most of the future backends), you will need to install “hard” dependencies as gcc backend or other: see the page dedicated to backends Backends.

Quickstart

As a complete example, we will show the usage of scikit-fdiff to solve the 2D shallow water equations.

We first import the library we need

>>> from skfdiff import Model, Simulation
>>> import pylab as pl
>>> import numpy as np
>>> from scipy.signal.windows import gaussian

We then write the model. We will use a simple periodic boundary condition

>>> model = Model(["-(dx((H + h) * u) + dy((H + h) * v))",
...                "-(u * dxu + v * dyu) - g * dxh + nu * (dxxu + dyyu)",
...                "-(u * dxv + v * dyv) - g * dyh + nu * (dxxv + dyyv)"],
...                ["h(x, y)", "u(x, y)", "v(x, y)"],
...                parameters=["H(x, y)", "nu", "g"],
...                boundary_conditions="periodic")

We set the initial condition

>>> L = 10
>>> x = y = np.linspace(-L / 2, L / 2, 56)
>>> xx, yy = np.meshgrid(x, y, indexing="ij")
>>> h = (gaussian(x.size, x.size // 20)[:, None] *
...      gaussian(y.size, y.size // 20)[None, :]) + 1
>>> h = np.roll(h, 12, axis=0)
>>> h = np.roll(h, 12, axis=1)
>>> H = np.zeros_like(h)
>>> u = np.zeros_like(h)
>>> v = np.zeros_like(h)
>>> initial_fields = model.Fields(x=x, y=y, h=h, u=u, v=v,
...                               H=H, g=9.81, nu=0)

We can run the simulation itself

>>> simulation = Simulation(model, initial_fields, dt=.1, tmax=1)
>>> container = simulation.attach_container()
>>> tmax, final_fields = simulation.run()

If we are in a notebook, we can use holoviews to plot the solution for each steps.

import holoviews as hv
hv.notebook_extension("bokeh")
hv.Dataset(container.data.h).to(hv.Image, ["x", "y"])