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.

Example: custom Python functions as C++ (grid) functions¶

So far, every coefficient function in these examples has been either a constant (ConstantFunction) or a parsed expression string (ExpressionFunction, e.g. 'exp(x[0]*x[1])'). The expression parser is convenient, but it is also a fragile input channel: it understands only what its grammar implements, and it silently mishandles some inputs (e.g. it rejects subnormal float literals).

dune.xt.functions.GenericFunction and GenericGridFunction close that gap: instead of a string to be parsed, they take a plain Python callable and invoke it directly from the C++ grid walk. This means any Python code – piecewise definitions, random coefficients, calls into numpy/scipy – can serve as a dune-gdt coefficient function.

GenericFunction is evaluated in global (physical) coordinates, like ExpressionFunction.
GenericGridFunction is evaluated locally, in reference-element coordinates, and additionally supports a post_bind callback invoked once per element during the grid walk – so the function can depend on the element itself (its index, its center, …), not just on the point being evaluated.

Both take an order (the polynomial degree used to pick a quadrature rule) exactly like ExpressionFunction, since the callback is opaque to the C++ side; there is no way to infer the polynomial degree of a Python function automatically.

# wurlitzer: display dune's output in the notebook
%load_ext wurlitzer

import time

import numpy as np

from dune.xt.grid import Dim

try:
    from dune.xt.functions import GenericFunction, GenericGridFunction

    HAVE_GENERIC_FUNCTION = True
except ImportError:
    HAVE_GENERIC_FUNCTION = False
    print(
        "This wheel was not built with the GenericFunction/GenericGridFunction bindings, so the\n"
        "cells below are skipped instead of failing."
    )

1: `GenericFunction` vs. `ExpressionFunction`¶

We solve the same stationary heat equation as in example__ipdg_stationary_heat_equation.md, once with the diffusion and source given as ExpressionFunction (a parsed string), and once with the same math given as a GenericFunction (a Python lambda). The two must produce the same discrete solution, since they describe the same function through two independent code paths.

d = 2
omega = ([0, 0], [1, 1])

if HAVE_GENERIC_FUNCTION:
    from dune.xt.functions import ConstantFunction, ExpressionFunction

    kappa_expression = ConstantFunction(
        dim_domain=Dim(d), dim_range=Dim(1), value=[1.0], name="kappa"
    )
    f_expression = ExpressionFunction(
        dim_domain=Dim(d),
        variable="x",
        expression="exp(x[0]*x[1])",
        order=3,
        name="f",
    )

    # the same source term, but evaluated by a Python callable instead of the expression parser
    def f_evaluate(x, mu):
        return [float(np.exp(x[0] * x[1]))]

    f_generic = GenericFunction(
        dim_domain=Dim(d), dim_range=Dim(1), order=3, evaluate=f_evaluate, name="f_generic"
    )

from dune.xt.grid import Simplex, make_cube_grid

if HAVE_GENERIC_FUNCTION:
    grid = make_cube_grid(Dim(d), Simplex(), lower_left=omega[0], upper_right=omega[1], num_elements=[2, 2])
    grid.global_refine(1)  # we need to refine once to obtain a symmetric grid
    print(f"grid has {grid.size(0)} elements, {grid.size(d - 1)} edges and {grid.size(d)} vertices")

grid has 16 elements, 28 edges and 13 vertices

if HAVE_GENERIC_FUNCTION:
    from discretize_elliptic_ipdg import discretize_elliptic_ipdg_dirichlet_zero

    # penalty_parameter is given explicitly (as in example__ipdg_stationary_heat_equation.md) to
    # avoid the LAPACKE-dependent auto-estimation path (estimate_combined_inverse_trace_inequality_
    # constant), which is not available in every build
    u_h_expression = discretize_elliptic_ipdg_dirichlet_zero(
        grid, kappa_expression, f_expression, penalty_parameter=16
    )
    u_h_generic = discretize_elliptic_ipdg_dirichlet_zero(
        grid, kappa_expression, f_generic, penalty_parameter=16
    )

    dofs_expression = np.array(u_h_expression.dofs.vector, copy=False)
    dofs_generic = np.array(u_h_generic.dofs.vector, copy=False)
    max_difference = np.abs(dofs_expression - dofs_generic).max()
    print(f"max |dofs_expression - dofs_generic| = {max_difference:.3e}")
    assert max_difference < 1e-13

max |dofs_expression - dofs_generic| = 0.000e+00

if HAVE_GENERIC_FUNCTION:
    from dune.gdt import visualize_function

    _ = visualize_function(u_h_generic)

2: `GenericGridFunction` with `post_bind` – an element-dependent coefficient¶

A GenericFunction can only see the (global) evaluation point, so it cannot depend on which element it is being evaluated on – e.g. to implement a piecewise-constant “checkerboard” coefficient without going through dune.xt.functions.CheckerboardFunction. GenericGridFunction’s post_bind callback is invoked once whenever the grid walk moves to a new element, before any evaluate calls on that element, so we can look up a per-element value there and stash it in a Python closure for evaluate to read:

if HAVE_GENERIC_FUNCTION:
    def make_checkerboard_diffusion(grid, blocks, values):
        """A piecewise-constant, element-wise isotropic diffusion tensor kappa(x) = values[cell(x)].

        `values` is indexed by the coarse "checkerboard" cell (row-major over `blocks`) the
        element's center falls into.
        """
        d = grid.dimension
        blocks = np.asarray(blocks[:d])
        state = {"value": values[0]}

        def post_bind(element):
            cell = np.minimum((np.asarray(element.center) * blocks).astype(int), blocks - 1)
            state["value"] = values[int(np.ravel_multi_index(tuple(cell), tuple(blocks)))]

        def evaluate(x, mu):
            value = float(state["value"])
            return [[value if i == j else 0.0 for j in range(d)] for i in range(d)]

        # dim_range disambiguates the scalar/vector/matrix GenericGridFunction overloads sharing
        # this factory name -- it is not otherwise inferable from the post_bind/evaluate callables
        return GenericGridFunction(
            grid,
            order=0,
            post_bind=post_bind,
            evaluate=evaluate,
            dim_range=(Dim(d), Dim(d)),
            name="checkerboard_kappa",
        )

    rng = np.random.default_rng(0)
    checkerboard_values = rng.uniform(0.1, 10.0, size=4)
    kappa_checkerboard = make_checkerboard_diffusion(grid, blocks=(2, 2), values=checkerboard_values)

    u_h_checkerboard = discretize_elliptic_ipdg_dirichlet_zero(
        grid, kappa_checkerboard, f_expression, penalty_parameter=16
    )
    _ = visualize_function(u_h_checkerboard)

3: overhead of a Python callback vs. a parsed expression¶

Every GenericFunction/GenericGridFunction evaluation calls back into Python, which means pybind11 has to acquire the GIL for each call (this happens automatically as part of the std::function <-> Python-callable conversion pybind11 generates – no code in the bindings manages the GIL explicitly). ExpressionFunction never leaves C++. We measure the per-call overhead by assembling the same Laplace matrix with each and dividing the wall-clock difference by the number of evaluate calls (elements times quadrature points times repetitions):

if HAVE_GENERIC_FUNCTION:
    from dune.gdt import (
        BilinearForm,
        ContinuousLagrangeSpace,
        LocalElementIntegralBilinearForm,
        LocalLaplaceIntegrand,
        MatrixOperator,
        make_element_sparsity_pattern,
    )
    from dune.xt.functions import GridFunction as GF

    def assemble_laplace(grid, diffusion, repetitions=5):
        space = ContinuousLagrangeSpace(grid, order=1)
        start = time.perf_counter()
        for _ in range(repetitions):
            form = BilinearForm(grid)
            form += LocalElementIntegralBilinearForm(
                LocalLaplaceIntegrand(GF(grid, diffusion, dim_range=(Dim(d), Dim(d))))
            )
            op = MatrixOperator(
                grid,
                source_space=space,
                range_space=space,
                sparsity_pattern=make_element_sparsity_pattern(space),
            )
            op.append(form)
            op.assemble()
        return time.perf_counter() - start

    fine_grid = make_cube_grid(Dim(d), Simplex(), lower_left=omega[0], upper_right=omega[1], num_elements=[8, 8])
    fine_grid.global_refine(1)
    num_elements = fine_grid.size(0)
    # order-1 CG on a 2d simplex: a 2nd order quadrature rule has 3 points per element
    quadrature_points_per_element = 3
    repetitions = 5

    time_expression = assemble_laplace(fine_grid, kappa_expression, repetitions)

    # scalar, matching kappa_expression exactly (GF's dim_range=(Dim(d), Dim(d)) tag in
    # assemble_laplace promotes both identically) -- so the timing difference below isolates the
    # cost of the Python callback itself, not an unrelated scalar-vs-matrix evaluation difference
    def kappa_evaluate(x, mu):
        return [1.0]

    kappa_generic = GenericFunction(
        dim_domain=Dim(d), dim_range=Dim(1), order=0, evaluate=kappa_evaluate, name="kappa_generic"
    )
    time_generic = assemble_laplace(fine_grid, kappa_generic, repetitions)

    total_calls = num_elements * quadrature_points_per_element * repetitions
    overhead_per_call = (time_generic - time_expression) / total_calls
    print(f"ExpressionFunction assembly: {time_expression:.4f} s")
    print(f"GenericFunction assembly:    {time_generic:.4f} s")
    print(f"~{total_calls} evaluate() calls total")
    print(f"estimated overhead per Python callback: {overhead_per_call * 1e6:.2f} us")

ExpressionFunction assembly: 0.0082 s
GenericFunction assembly:    0.0102 s
~3840 evaluate() calls total
estimated overhead per Python callback: 0.53 us

This overhead (typically low single-digit microseconds per call on CPython, dominated by the GIL acquisition and the argument/return-value conversion) is the price of the flexibility gained: any Python code, however irregular, becomes usable as a coefficient function.

Download the code: example__custom_python_functions.md example__custom_python_functions.ipynb

Example: custom Python functions as C++ (grid) functions¶

1: GenericFunction vs. ExpressionFunction¶

2: GenericGridFunction with post_bind – an element-dependent coefficient¶

3: overhead of a Python callback vs. a parsed expression¶

1: `GenericFunction` vs. `ExpressionFunction`¶

2: `GenericGridFunction` with `post_bind` – an element-dependent coefficient¶