Tutorial: Using pyMOR’s discretization toolkit
==============================================

.. include:: jupyter_init.txt

pyMOR’s discretization toolkit allows to quickly build parameterized
full-order models based on the NumPy/SciPy software stack. Currently
supported are stationary or instationary scalar PDEs of up to second
order with possibly nonlinear advection or reaction terms in one or two
spatial dimensions. Unstructured meshes can be imported in the Gmsh file
format.

In this tutorial we will focus on elliptic equations of the form

.. math::

   -\nabla \cdot \big(\sigma(x, \mu) \nabla u(x, \mu) \big) = f(x, \mu),\quad x \in \Omega,

on the domain :math:`\Omega:= (0, 1)^2 \subset \mathbb{R}^2` with data
functions :math:`f(\cdot, \mu) \in L^2(\Omega)`,
:math:`\sigma(\cdot, \mu) \in L^\infty(\Omega)`.


A first equation without parameters
-----------------------------------

First, let us assume that the source :math:`f(x, \mu)` is an indicator
function of a circular disk with radius :math:`0.3` and that
:math:`\sigma(x, \mu)` is constant:

.. math::

   f(x, \mu) :=
   \begin{cases}
      1, & |x - (0.5, 0.5)| < 0.3, \\
      0, & \text{otherwise},
   \end{cases} \quad\text{and}\quad
   \sigma(x, \mu) :\equiv 1.

We start by importing commonly used pyMOR classes and methods from the
:mod:`~pymor.basic` module:

.. jupyter-execute::

   from pymor.basic import *

To specify the problem at hand using pyMOR’s discretization toolkit, we
first need to specify the computational domain :math:`\Omega`. Multiple
classes are available to define such domains in the
:mod:`~pymor.analyticalproblems.domaindescriptions` module,
which all derive from the |DomainDescription| interface class.

In our case, we can use a |RectDomain|:

.. jupyter-execute::

    domain = RectDomain([[0.,0.], [1.,1.]])

Data functions are defined using classes which derive from
the |Function| interface. We specify the constant diffusivity :math:`\sigma`
using a |ConstantFunction|:

.. jupyter-execute::

    diffusion = ConstantFunction(1, 2)

Here, the first argument is the function’s constant value. The second
argument is the spatial dimension of the domain the problem is defined
on.

For the definition of the source term :math:`f` we use an
|ExpressionFunction| which is given an arbitrary Python expression
used to evaluate the function. In this expression, the coordinates at
which the function shall be evaluated are given as the variable `x`.
Many NumPy functions can be used directly. The entire NumPy module is
available under the name `np`.

Thus, to define :math:`f` we could write ::

   (sqrt((x[0]-0.5)**2 + (x[1]-0.5)**2) <= 0.3) * 1.

However, pyMOR |Functions| are required to be vectorized with respect
to the coordinate `x`. In the case of |ExpressionFunction| this
means that `x` can be an arbitrary dimensional NumPy array of
coordinates where the last array index specifies the spacial dimension.
Therefore, the correct definition of :math:`f` is:

.. jupyter-execute::

   rhs = ExpressionFunction('(sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 1.', 2, ())

Similarly to |ConstantFunction|, the second argument is the dimension
of the computational domain. As the shape of the return value cannot be
easily inferred from the given string expression, it has to be provided
as a third argument to |ConstantFunction|. For scalar functions we
provide the empty tuple `()`, for functions returning
three-dimensional vectors we would specify `(3,)`, and for functions
returning :math:`2\times 2` matrices we would specify `(2,2)`.

Finally, the computational domain and all data functions are collected
in a |StationaryProblem|:

.. jupyter-execute::

   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       rhs=rhs,
   )

This problem definition can now be handed over to discretization methods
(see :mod:`pymor.discretizers.builtin`) which assemble corresponding
discrete full-order models. For finite elements, we use
:func:`~pymor.discretizers.builtin.cg.discretize_stationary_cg`,
which receives the maximum mesh element diameter via the `diameter`
argument:

.. jupyter-execute::

   m, data = discretize_stationary_cg(problem, diameter=1/4)

The resulting |Model| can be :meth:`solved <pymor.models.interface.Model.solve>`,
returning a |VectorArray| with the solution data:

.. jupyter-execute::

   U = m.solve()

Finally, we visualize the solution:

.. jupyter-execute::

   m.visualize(U)

In case a specific grid type shall be used (|RectGrid| or
|TriaGrid|), the corresponding class has to be passed to the
discretizer as the `grid_type` argument. By using |RectGrid| we get
bilinear finite elements:

.. jupyter-execute::

   m, data = discretize_stationary_cg(problem, diameter=1/4, grid_type=RectGrid)
   m.visualize(m.solve())

We get a finite volume model using
:func:`~pymor.discretizers.builtin.fv.discretize_stationary_fv`:

.. jupyter-execute::

   m, data = discretize_stationary_fv(problem, diameter=1/4, grid_type=TriaGrid)
   m.visualize(m.solve())


Defining boundary conditions
----------------------------

As the vigilant reader will already have noticed, we did not specify any
boundary conditions when defining and solving our problem. When no
boundary conditions are specified, pyMOR’s discretization toolkit will
assume that homogeneous Dirichlet conditions are implied over the entire
boundary of :math:`\Omega`.

As the next example, let us now assume that the data functions are given
by

.. math::

   f(x, \mu) :\equiv 0 \quad\text{and}\quad
   \sigma(x, \mu) :=
   \begin{cases}
      0.001, & |x - (0.5, 0.5)| < 0.3, \\
      1, & \text{otherwise},
   \end{cases}

and that we have the following mixed boundary conditions

.. math::

   \begin{align}
    -\sigma(x, \mu) \nabla u(x, \mu) \cdot n &= g_N(x), && x \in (0,1) \times \{0\} =: \Omega_N \\
    u(x, \mu) &= 0, && x \in \partial\Omega \setminus \Omega_N,
   \end{align}

with :math:`g_N(x) \equiv -1`.

Before solving this problem, let us first silence pyMOR’s verbose log
messages for the rest of this tutorial using the :func:`~pymor.core.logger.set_log_levels`
method:

.. jupyter-execute::

   set_log_levels({'pymor': 'WARN'})

To impose the correct boundary conditions we need to declare which type of
boundary condition should be active on which part of
:math:`\partial\Omega` when defining the computational domain:

.. jupyter-execute::

   domain = RectDomain(bottom='neumann')

Then all we need is to pass the Neumann data function :math:`g_N` to the
|StationaryProblem|. Here, we can use again a |ConstantFunction|.
The diffusivity can be defined similarly as above:

.. jupyter-execute::

   neumann_data = ConstantFunction(-1., 2)

   diffusion = ExpressionFunction('1. - (sqrt( (x[...,0]-0.5)**2 + (x[...,1]-0.5)**2) <= 0.3) * 0.999' , 2, ())

   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=neumann_data
   )

Finally, we discretize and solve:

.. jupyter-execute::

   m, data = discretize_stationary_cg(problem, diameter=1/32)
   m.visualize(m.solve())


Another example
---------------

Even with a single |ExpressionFunction| we can build many different examples.
For instance, to let :math:`\sigma` be given by a periodic pattern of
:math:`K\times K` circular disks of radius :math:`0.3/K` we can use the
following definition:

.. jupyter-execute::

   diffusion = ExpressionFunction(
       '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
       2, (),
       values={'K': 10}
   )


Here, we have used the `values` parameter of |ExpressionFunction| to
make `K` available as an additional constant in the defining
expression. In particular, we can easily change `K` programatically
without having to resort to string manipulations. The solution looks
like this:

.. jupyter-execute::

   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=neumann_data
   )


   m, data = discretize_stationary_cg(problem, diameter=1/100)
   m.visualize(m.solve())


Data functions defined from pixel graphics
------------------------------------------

|BitmapFunction| uses the Python Imaging Library (PIL) to read gray
scale images in various image file formats. The resulting
two-dimensional NumPy array of pixel values defines a piecewise constant
data function on a rectangular domain, where the range of the function
(from black to white) is specified via the `range` parameter. For
instance, when using a |BitmapFunction| for :math:`\sigma` with the
following graphic stored in `RB.png`:

.. image:: RB.png

and a range of `[0.001 1]` we obtain:

.. jupyter-execute::

   diffusion = BitmapFunction('RB.png', range=[0.001, 1])
   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=neumann_data
   )

   m, data = discretize_stationary_cg(problem, diameter=1/100)
   m.visualize(m.solve())

The displayed warning results from the fact that the used image file has
an additional channel of transparency values (alpha channel) and can be
ignored.


A parametric problem
--------------------

Now, let us consider the Neumann data function:

.. math::

   g_N((x_0, x_1), \mu_{neum}) := -\cos(\pi \cdot x_0)^2 \cdot\mu_{neum}

with a single |Parameter| :math:`\mu_{neum} \in \mathbb{R}`.

In pyMOR, |parameter values| are specified as a dictionary of one-dimensional
|NumPy arrays|. Each value in the dictionary must have a correct size specified by the
|Parameters| of the |ParametricObject|. In this example we have a single scalar
valued parameter which we call `'neum'`. Thus, the |Parameters| of the function
will be ::

   {'neum': 1}

We can then make the following definition of the Neumann data:

.. jupyter-execute::

   neumann_data = ExpressionFunction('-cos(pi*x[...,0])**2*neum[0]', 2, (), parameters= {'neum': 1})

Similar to the range of the function, pyMOR cannot infer from the given
string expression the parameters used in the expression, so these
|Parameters| have to be provided as the `parameters` argument.
The individual parameters are then available as variables in
the expression.

We can then proceed as usual and automatically obtain a parametric
|Model|:

.. jupyter-execute::

   diffusion = ExpressionFunction(
       '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * 0.999',
       2, (),
       values={'K': 10}
   )
   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=neumann_data
   )

   m, data = discretize_stationary_cg(problem, diameter=1/100)
   m.parameters

When solving the model, we now need to specify appropriate
|parameter values|:

.. jupyter-execute::

   m.visualize(m.solve({'neum': [1.]}))

For the :meth:`~pymor.models.interface.Model.solve` method, the
parameter value can also be specified as a single number:

.. jupyter-execute::

   m.visualize(m.solve(-100))


Multiple parameters
-------------------

Next we also want to parameterize the diffusivity in the
:math:`K \times K` circular disks by a scalar factor
:math:`\mu_{diffu}`. To this end we define:

.. jupyter-execute::

   diffusion = ExpressionFunction(
       '1. - (sqrt( (np.mod(x[...,0],1./K)-0.5/K)**2 + (np.mod(x[...,1],1./K)-0.5/K)**2) <= 0.3/K) * (1 - diffu[0])',
       2, (),
       values={'K': 10},
       parameters= {'diffu': 1}
   )

We proceed as usual:

.. jupyter-execute::

   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=neumann_data
   )

   m, data = discretize_stationary_cg(problem, diameter=1/100)
   m.parameters

As we can see, pyMOR automatically derives that in this case the model
depends on two |Parameters|, and we have to provide two values
when solving the model:

.. jupyter-execute::

   m.visualize(m.solve({'diffu': 0.001, 'neum': 1}))

For :meth:`~pymor.models.interface.Model.solve` we can also
simply pass a list of parameter values, in which case
pyMOR assumes an alphabetical ordering of the parameters:

.. jupyter-execute::

   m.visualize(m.solve([1, -1]))


Parameter-separability
----------------------

For the generation of online-efficient reduced-order models, it is often
crucial that data functions are parameter separable. We call a
parametric function `f(x, \mu)` parameter separable if it admits a
decomposition

.. math::

   f(x, \mu) = \sum_{q=1}^Q f_q(x) \cdot \theta_q(\mu)

where the :math:`f_q` are non-parametric and the *parameter
functionals* :math:`\theta_q` map parameters to real numbers.

To model such a decomposition using pyMOR’s discretization toolkit, we
specify the :math:`f_q` as non-parametric |Functions|, the
:math:`\theta_q` are represented by |ParameterFunctionals| and the
given sum decomposition is represented by a |LincombFunction| of these
objects.

As an example let us go back to the case where the diffusivity is
represented by indicator functions of point sets forming the letters
`RB`. We now want to parameterize the diffusivity in the individual
letters. This admits a decomposition of the form

.. math::

   \sigma(x,y) = 1 + f_R \cdot (\mu_R - 1) + f_B \cdot (\mu_L - 1)

Again, we define :math:`f_R` and :math:`f_L` as |BitmapFunctions| for
the following image files:

.. image:: R.png

.. image:: B.png

.. jupyter-execute::

   f_R = BitmapFunction('R.png', range=[1, 0])
   f_B = BitmapFunction('B.png', range=[1, 0])

Next we need to define the |ParameterFunctionals|

.. math::

   \theta_R(\mu) = \mu_R - 1 \quad\text{and}\quad \theta_B(\mu) = \mu_B - 1.

Similar to an |ExpressionFunction|, we can use
|ExpressionParameterFunctionals| for that:

.. jupyter-execute::

   theta_R = ExpressionParameterFunctional('R[0] - 1', {'R': 1})
   theta_B = ExpressionParameterFunctional('B[0] - 1', {'B': 1})

Note that the second argument is again the |Parameters|
that are used in the expression. Finally, we form the linear
combination using a |LincombFunction| which is given a list of
|Functions| as the first and a corresponding list of
|ParameterFunctionals| or constants as the second argument:

.. jupyter-execute::

   diffusion = LincombFunction(
       [ConstantFunction(1., 2), f_R, f_B],
       [1., theta_R, theta_B]
   )
   diffusion.parameters

Again, pyMOR automatically derives that the evaluation of `diffusion`
depends on the two |Parameters| `'B'` and `'R'`. Now, we can
proceed as usual:

.. jupyter-execute::

   problem = StationaryProblem(
       domain=domain,
       diffusion=diffusion,
       neumann_data=ConstantFunction(-1, 2)
   )
   m, data = discretize_stationary_cg(problem, diameter=1/100)
   m.visualize((m.solve([1., 0.001]), m.solve([0.001, 1])))

Looking at the |Model| `m`, we can see that the decomposition of
:math:`\sigma` has been preserved by the discretizer:

.. jupyter-execute::

   m.operator

The |LincombFunction| has become a |LincombOperator|, with the same
linear coefficients but the |BitmapFunctions| replaced by
corresponding stiffness matrices. Note that an additional summand
appears which ensures correct enforcement of Dirichlet boundary values
for all possible parameter value combinations.

Download the code:
:jupyter-download:script:`tutorial_builtin_discretizer`
:jupyter-download:notebook:`tutorial_builtin_discretizer`