Run this tutorial

Click here to run this tutorial on mybinder.org:
Please note that starting the notebook server may take a couple of minutes.

Tutorial: Binding an external PDE solver to pyMOR

One of pyMOR’s main features is easy integration of external solvers that implement the full-order model. In this tutorial we will do this step-by-step for a custom toy solver written in C++. If you use the FEniCS or NGSovle PDE solver libraries, you can find ready-to-use pyMOR bindings in the bindings package. pyMOR support for deal.II can be found in a separate repository.

Defining the PDE solver

Our solver discretizes the one-dimensional Laplace equation \(u''(x)=0\) on the interval \([\text{left},\text{right}]\) using a central differences scheme with \(h=\frac{|\text{right}-\text{left}|}{n}\). First, we need a class to store our data in and with some basic linear algebra operations declared on it.

#include <vector>


class Vector {
  friend class DiffusionOperator;
public:
  Vector(int dim, double value);
  Vector(const Vector& other);
  const int dim;
  void scal(double val);
  void axpy(double a, const Vector& x);
  double inner(const Vector& other) const;
  double* data();
private:
  std::vector<double> _data;
};

Next, we need the operator that discretizes our PDE.

class DiffusionOperator {
public:
  DiffusionOperator(int n, double left, double right);
  const int dim_source;
  const int dim_range;
  void apply(const Vector& U, Vector& R) const;
private:
  const double h;
  const double left;
  const double right;
};

Together with some header guards, these two snippets make up our model.hh.

The definitions for the Vector class are pretty straightforward:

Vector::Vector(int dim, double value) : _data(dim, value), dim(dim) {}

Vector::Vector(const Vector& other) : _data(other._data), dim(other.dim) {}

void Vector::scal(double val) {
  for (int i = 0; i < dim; i++) {
    _data[i] *= val;
  }
}

void Vector::axpy(double a, const Vector& x) {
  assert(x.dim == dim);
  for (int i = 0; i < dim; i++) {
    _data[i] += a * x._data[i];
  }
}

double Vector::inner(const Vector& other) const {
  assert(other.dim == dim);
  double result = 0;
  for (int i = 0; i < dim; i++) {
    result += _data[i] * other._data[i];
  }
  return result;
}

double* Vector::data() {
  return _data.data();
}

Just like the diffusion operator that computes a central differences stencil:

DiffusionOperator::DiffusionOperator(int n, double left, double right)
  : dim_source(n+1), dim_range(n+1), h(1./n), left(left), right(right) {}

void DiffusionOperator::apply(const Vector& U, Vector& R) const {
  const double one_over_h2 = 1. / (h * h);
  for (int i = 1; i < dim_range - 1; i++) {
    if ((i * h > left) && (i * h <= right)) {
      R._data[i] = -(U._data[i-1] - 2*U._data[i] + U._data[i+1]) * one_over_h2;
    }
  }
}

This completes all the C++ code needed for the toy solver itself. Next, we will make this code usable from Python. We utilize the pybind11 library to create a Python extension module named model, that allows us to manipulate instances of the C++ Vector and DiffusionOperator classes.

Compiling the PDE solver as a shared library and creating Python bindings for it using pybind11, Cython or ctypes is the preferred way of integrating external solvers, as it offers maximal flexibility and performance. For instance, in this example we will actually completely implement the Model in Python using a time stepper from pyMOR to solve the Model.

When this is not an option, RPC-based approaches are possible as well. For small to medium-sized linear problems, another option it to import system matrices and snapshot data into pyMOR via file exchange and to use NumPy-based Operators and VectorArrays to represent the full-order model.

Binding the solver to Python

All of the C++ code related to the extension module is defined inside a scope started with

PYBIND11_MODULE(model, m)
{

This tells pybind11 to make the contained symbols accessible in the module instance m that will be importable by the name model. Now we create a new pybind11 class\_ object that wraps the DiffusionOperator. Note that the module instance is passed to the constructor alongside a name for the Python class and a docstring. The second line shows how to define an init function for the Python object by using the special py:init object to forward arguments to the C++ constructor.

    py::class_<DiffusionOperator> op(m, "DiffusionOperator", "DiffusionOperator Docstring");
    op.def(py::init<int, double, double>());

Next, we define read-only properties on the Python side named after and delegated to the members of the C++ class.

    op.def_readonly("dim_source", &DiffusionOperator::dim_source);
    op.def_readonly("dim_range", &DiffusionOperator::dim_range);

The last DiffusionOperator-related line exposes the function call to apply in the same way:

    op.def("apply", &DiffusionOperator::apply);

This is everything that is needed to expose the operator to Python. We will now do the same for the Vector, with a few more advanced techniques added.

    py::class_<Vector> vec(m, "Vector", "Vector Docstring", py::buffer_protocol());
    vec.def(py::init([](int size, double value) { return std::make_unique<Vector>(size, value); }));
    vec.def(py::init([](const Vector& other) { return std::make_unique<Vector>(other); }));

Again we define a py:class\_ with appropiate name and docstring, but now we also indicate to pybind11 that this class will implement the buffer protocol, which basically exposes direct access to the chunk of memory associated with a Vector instance to Python. We also see how we can dispatch multiple init functions by using py:init objects with C++ lambda functions. Note that direct memory access to the vector data from Python is not required to integrate a solver with pyMOR. It is, however, useful for debugging and quickly modifying or extending the solver from within Python. For instance, in our toy example we will use the direct memory access to quickly define a visualization of the solutions and to construct the right-hand side vector for our problem.

    vec.def_readonly("dim", &Vector::dim);
    vec.def("scal", &Vector::scal);
    vec.def("axpy", &Vector::axpy);
    vec.def("inner", &Vector::inner);
    vec.def("data", &Vector::data);

    vec.def_buffer([](Vector& vec) -> py::buffer_info {
        return py::buffer_info(
            vec.data(), sizeof(double), py::format_descriptor<double>::format(), 1, {vec.dim}, {sizeof(double)});
    });
} // end PYBIND11_MODULE(model, m)

This completes the model.cc.

This extension module needs to be compiled to a shared object that the Python interpreter can import. We use a minimal CMake project that generates makefiles for us to achieve this.

First we make sure pybind11 can be used:

cmake_minimum_required(VERSION 3.1)

set(Python3_EXECUTABLE /usr/local/bin/python3 CACHE PATH "")
project(minimal-cpp-demo)

if(${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.18.0")

Next, we define a new library with our model.cc as the single source file and let pybind11 set the proper compile flags.

  find_package(Python3 COMPONENTS Interpreter Development REQUIRED)
endif()

find_package(pybind11 CONFIG REQUIRED HINTS ${PYBIND11_DIR} ${PYBIND11_ROOT}

That is all that is needed for CMakeLists.txt. In the next step, we will switch to a bash terminal and actually compile this module.

After creating a build directory for the module, we let cmake initialize the build and call make to execute the compilation.

%%bash
mkdir -p minimal_cpp_demo/build
cmake -B minimal_cpp_demo/build -S minimal_cpp_demo
make -C minimal_cpp_demo/build

-- The C compiler identification is GNU 10.2.1

-- The CXX compiler identification is GNU 10.2.1

-- Detecting C compiler ABI info

-- Detecting C compiler ABI info - done

-- Check for working C compiler: /usr/lib/ccache/cc - skipped

-- Detecting C compile features

-- Detecting C compile features - done

-- Detecting CXX compiler ABI info

-- Detecting CXX compiler ABI info - done

-- Check for working CXX compiler: /usr/lib/ccache/c++ - skipped

-- Detecting CXX compile features

-- Detecting CXX compile features - done

-- Found Python3: /usr/local/bin/python3 (found version "3.9.16") found components: Interpreter Development.Module 

-- Performing Test HAS_FLTO

-- Performing Test HAS_FLTO - Success

-- Found pybind11: /usr/local/lib/python3.9/site-packages/pybind11/include (found version "2.9.2")

-- Found pybind11 v2.9.2: /usr/local/lib/python3.9/site-packages/pybind11/include;/usr/local/include/python3.9

-- Configuring done

-- Generating done

-- Build files have been written to: /builds/pymor/pymor/docs/source/minimal_cpp_demo/build

make[2]: Entering directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

make[3]: Entering directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

make[4]: Entering directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

make[4]: Leaving directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

make[4]: Entering directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

[ 50%] Building CXX object CMakeFiles/model.dir/model.cc.o

[100%] Linking CXX shared module model.cpython-39-x86_64-linux-gnu.so

make[4]: Leaving directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

[100%] Built target model

make[3]: Leaving directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

make[2]: Leaving directory '/builds/pymor/pymor/docs/source/minimal_cpp_demo/build'

To be able to use this extension module we need to insert the build directory into the path where the Python interpreter looks for things to import. Afterwards we can import the module and create and use the exported classes.

import sys
sys.path.insert(0, 'minimal_cpp_demo/build')

import model
mymodel = model.DiffusionOperator(10, 0, 1)
myvector = model.Vector(10, 0)
mymodel.apply(myvector, myvector)
dir(model)

['DiffusionOperator',
 'Vector',
 '__doc__',
 '__file__',
 '__loader__',
 '__name__',
 '__package__',
 '__spec__']

Using the exported Python classes with pyMOR

All of pyMOR’s algorithms operate on VectorArray and Operator objects that all share the same programming interface. To be able to use our Python model.Vector and model.DiffusionOperator in pyMOR, we have to provide implementations of VectorArray, VectorSpace and Operator that wrap the classes defined in the extension module and translate calls to the interface methods into operations on model.Vector and model.DiffusionOperator.

Instead of writing a full implementaion of a VectorArray that manages multiple model.Vector instances, we can instead implement a wrapper WrappedVector for a single model.Vector instance based on CopyOnWriteVector which will be used to create ListVectorArrays via a ListVectorSpace-based WrappedVectorSpace.

The CopyOnWriteVector base class manages a reference count for us and automatically copies data when necessary in methods scal and axpy. To use this, we need to implement _scal and _axpy in addition to all the abstract methods from CopyOnWriteVector. We can get away with using just a stub that raises an NotImplementedError in some methods that are not actually called in our example.

from pymor.operators.interface import Operator
from pymor.vectorarrays.list import CopyOnWriteVector, ListVectorSpace

import numpy as np
import math
from model import Vector, DiffusionOperator


class WrappedVector(CopyOnWriteVector):

    def __init__(self, vector):
        assert isinstance(vector, Vector)
        self._impl = vector

    @classmethod
    def from_instance(cls, instance):
        return cls(instance._impl)

    def to_numpy(self, ensure_copy=False):
        # Note how this uses the buffer protocol setup to allow efficient
        # data access as a Numpy Vector
        result = np.frombuffer(self._impl, dtype=np.float64)
        if ensure_copy:
            result = result.copy()
        return result

    def _copy_data(self):
        # This uses the second exposed 'init' signature to delegate to the C++ copy constructor
        self._impl = Vector(self._impl)

    def _scal(self, alpha):
        self._impl.scal(alpha)

    def _axpy(self, alpha, x):
        self._impl.axpy(alpha, x._impl)

    def inner(self, other):
        return self._impl.inner(other._impl)

    def norm(self):
        return math.sqrt(self.inner(self))

    def norm2(self):
        return self.inner(self)

    def sup_norm(self):
        raise NotImplementedError

    def dofs(self, dof_indices):
        raise NotImplementedError

    def amax(self):
        raise NotImplementedError

The implementation of the WrappedVectorSpace is very short as most of the necessary methods of VectorSpace are implemented in ListVectorSpace.

class WrappedVectorSpace(ListVectorSpace):

    def __init__(self, dim):
        self.dim = dim

    def zero_vector(self):
        return WrappedVector(Vector(self.dim, 0))

    def make_vector(self, obj):
        # obj is a `model.Vector` instance
        return WrappedVector(obj)

    def __eq__(self, other):
        return type(other) is WrappedVectorSpace and self.dim == other.dim

Wrapping the model.DiffusionOperator is straightforward as well. We just need to attach suitable VectorSpaces to the class and implement the application of the operator on a VectorArray as a sequence of applications on single vectors.

class WrappedDiffusionOperator(Operator):
    def __init__(self, op):
        assert isinstance(op, DiffusionOperator)
        self.op = op
        self.source = WrappedVectorSpace(op.dim_source)
        self.range = WrappedVectorSpace(op.dim_range)
        self.linear = True

    @classmethod
    def create(cls, n, left, right):
        return cls(DiffusionOperator(n, left, right))

    def apply(self, U, mu=None):
        assert U in self.source

        def apply_one_vector(u):
            v = Vector(self.range.dim, 0)
            self.op.apply(u._impl, v)
            return v

        return self.range.make_array([apply_one_vector(u) for u in U.vectors])

Putting it all together

As a demonstration, we will use our toy Laplace solver to compute an approximation for the transient diffusion equation

\[\frac{\partial u}{\partial t} = {\alpha_\mu} \frac{\partial^2 u}{\partial x^2},\]

with explicit timestepping provided by pyMOR, with a parameterized, block-wise defined, diffusion coefficient \(\alpha_\mu\).

First up, we implement a discretize function that uses the WrappedDiffusionOperator and WrappedVectorSpace to assemble an InstationaryModel.

from pymor.algorithms.pod import pod
from pymor.algorithms.timestepping import ExplicitEulerTimeStepper
from pymor.discretizers.builtin.gui.visualizers import OnedVisualizer
from pymor.models.basic import InstationaryModel
from pymor.discretizers.builtin import OnedGrid
from pymor.operators.constructions import VectorOperator, LincombOperator
from pymor.parameters.functionals import ProjectionParameterFunctional
from pymor.reductors.basic import InstationaryRBReductor


def discretize(n, nt, blocks):
    h = 1. / blocks
    ops = [WrappedDiffusionOperator.create(n, h * i, h * (i + 1)) for i in range(blocks)]
    pfs = [ProjectionParameterFunctional('diffusion_coefficients', blocks, i) for i in range(blocks)]
    operator = LincombOperator(ops, pfs)

    initial_data = operator.source.zeros()

    rhs_vec = operator.range.zeros()
    rhs_data = rhs_vec.vectors[0].to_numpy()
    rhs_data[:] = np.ones(len(rhs_data))
    rhs_data[0] = 0
    rhs_data[len(rhs_data) - 1] = 0
    rhs = VectorOperator(rhs_vec)

    # we can re-use pyMOR's builtin grid and visualizer for our demonstration
    grid = OnedGrid(domain=(0, 1), num_intervals=n)
    visualizer = OnedVisualizer(grid)

    time_stepper = ExplicitEulerTimeStepper(nt)

    fom = InstationaryModel(T=1, operator=operator, rhs=rhs, initial_data=initial_data,
                            time_stepper=time_stepper, num_values=20,
                            visualizer=visualizer, name='C++-Model')
    return fom

Now we can build a reduced basis for our model. Note that this code is not specific to our wrapped classes. Those wrapped classes are only directly used in the discretize call.

%matplotlib inline
# discretize
fom = discretize(50, 10000, 4)
parameter_space = fom.parameters.space(0.1, 1)

# generate solution snapshots
snapshots = fom.solution_space.empty()
for mu in parameter_space.sample_uniformly(2):
    snapshots.append(fom.solve(mu))

# apply POD
reduced_basis = pod(snapshots, modes=4)[0]

# reduce the model
reductor = InstationaryRBReductor(fom, reduced_basis, check_orthonormality=True)
rom = reductor.reduce()

# stochastic error estimation
mu_max = None
err_max = -1.
for mu in parameter_space.sample_randomly(10):
    U_RB = reductor.reconstruct(rom.solve(mu))
    U = fom.solve(mu)
    err = np.max((U_RB-U).norm())
    if err > err_max:
        err_max = err
        mu_max = mu

# visualize maximum error solution
U_RB = (reductor.reconstruct(rom.solve(mu_max)))
U = fom.solve(mu_max)
fom.visualize((U_RB, U), title=f'mu = {mu}', legend=('reduced', 'detailed'))

Once Loop Reflect

Once Loop Reflect

As you can see in this comparison, we get a good approximation of the full-order model here and the error plot confirms it:

fom.visualize((U-U_RB), title=f'mu = {mu}', legend=('error'))

Once Loop Reflect

You can download this demonstration plus the wrapper definitions as a notebook tutorial_external_solver.ipynb or as Markdown text tutorial_external_solver.md.