Run this tutorial

Click here to run this tutorial on mybinder.org:# Tutorial: Building a Reduced BasisÂ¶

In this tutorial we will learn more about `VectorArrays`

and how to
construct a reduced basis using pyMOR.

A reduced basis spans a low dimensional subspace of a `Model`

â€™s
`solution_space`

, in which the
`solutions`

of the `Model`

can be well approximated for all `parameter values`

. In this context,
time is treated as an additional parameter. So for time-dependent problems,
the reduced space (the span of the reduced basis) should approximate the
solution for all `parameter values`

and time instances.

An upper bound for the possible quality of a reduced space is given by the so-called Kolmogorov \(N\)-width \(d_N\) given as

In this formula \(V\) denotes the
`solution_space`

, \(\mathcal{P} \subset \mathbb{R}^p\)
denotes the set of all `parameter values`

we are interested in, and
\(u(\mu)\) denotes the `solution`

of the `Model`

for the given `parameter values`

\(\mu\).
In pyMOR the set \(\mathcal{P}\) is called the `ParameterSpace`

.

How to read this formula? For each candidate reduced space \(V_N\) we
look at all possible `parameter values`

\(\mu\) and compute the best-approximation
error in \(V_N\) (the second infimum). The supremum over the infimum
is thus the worst-case best-approximation error over all `parameter values`

of
interest. Now we take the infimum of the worst-case best-approximation errors
over all possible reduced spaces of dimension at most \(N\), and this is
\(d_N\).

So whatever reduced space of dimension \(N\) we pick, we will always find a \(\mu\) for which the best-approximation error in our space is at least \(d_N\). Reduced basis methods aim at constructing spaces \(V_N\) for which the worst-case best-approximation error is as close to \(d_N\) as possible.

However, we will only find a good \(V_N\) of small dimension \(N\) if the values \(d_N\) decrease quickly for growing \(N\). It can be shown that this is the case as soon as \(u(\mu)\) analytically depends on \(\mu\), which is true for many problems of interest. More precisely, it can be shown [BCD+11], [RD13] that there are constants \(C, c > 0\) such that

In this tutorial we will construct reduced spaces \(V_N\) for a concrete problem with pyMOR and study their error decay.

## Model setupÂ¶

First we need to define a `Model`

and a `ParameterSpace`

for which we want
to build a reduced basis. We choose here the standard
`thermal block`

benchmark
problem shipped with pyMOR (see Getting started). However, any pyMOR
`Model`

can be used, except for Section Weak Greedy Algorithm where
some more assumptions have to be made on the `Model`

.

First we import everything we need:

```
import numpy as np
from pymor.basic import *
```

```
/builds/pymor/pymor/src/pymor/core/config.py:49: UserWarning: dune-gdt bindings have been tested for version 2021.1.x (x >= 2) (installed: 2022.1.1.603).
warnings.warn('dune-gdt bindings have been tested for version 2021.1.x (x >= 2) '
/builds/pymor/pymor/src/pymor/core/config.py:57: UserWarning: dune-gdt bindings have been tested for dune-xt 2021.1.x (x >= 2) (installed: 2022.1.4.601).
warnings.warn('dune-gdt bindings have been tested for dune-xt 2021.1.x (x >= 2) '
```

Then we build a 3-by-3 thermalblock problem that we discretize using pyMORâ€™s
`builtin discretizers`

(see
Tutorial: Using pyMORâ€™s discretization toolkit for an introduction to pyMORâ€™s discretization toolkit).

```
problem = thermal_block_problem((3,3))
fom, _ = discretize_stationary_cg(problem, diameter=1/100)
```

Next, we need to define a `ParameterSpace`

of `parameter values`

for which
the solutions of the full-order model `fom`

should be approximated by the reduced basis.
We do this by calling the `space`

method
of the |parameters| of `fom`

:

```
parameter_space = fom.parameters.space(0.0001, 1.)
```

Here, `0.0001`

and `1.`

are the common lower and upper bounds for the
individual components of all |parameters| of `fom`

. In our case `fom`

expects a single parameter `'diffusion'`

of 9 values:

```
fom.parameters
```

```
Parameters({diffusion: 9})
```

If `fom`

were to depend on multiple |parameters|, we could also call the
`space`

method with a dictionary
of lower and upper bounds per parameter.

The main use of `ParameterSpaces`

in pyMOR is that they allow to easily sample
`parameter values`

from their domain using the methods
`sample_uniformly`

and
`sample_randomly`

.

## Computing the snapshot dataÂ¶

Reduced basis methods are snapshot-based, which means that they build
the reduced space as a linear subspace of the linear span of solutions
of the `fom`

for certain `parameter values`

. The easiest
approach is to just pick these values randomly, what we will do in the
following. First we define a training set of 25 parameters:

```
training_set = parameter_space.sample_randomly(25)
print(training_set)
```

```
[Mu({'diffusion': array([0.37460266, 0.95071923, 0.73202074, 0.59869862, 0.15610304,
0.15607892, 0.0581778 , 0.86618953, 0.6011549 ])}), Mu({'diffusion': array([0.70810177, 0.02068244, 0.96991286, 0.8324594 , 0.21241788,
0.18190678, 0.18348617, 0.30431182, 0.52480396])}), Mu({'diffusion': array([0.43200182, 0.29130002, 0.61189171, 0.13957991, 0.29221543,
0.36642521, 0.45612438, 0.78519744, 0.19975381])}), Mu({'diffusion': array([0.51428301, 0.59245533, 0.04654577, 0.6075841 , 0.17060707,
0.06514509, 0.94889065, 0.96563547, 0.80841651])}), Mu({'diffusion': array([0.30468331, 0.09776235, 0.6842646 , 0.44020848, 0.12212603,
0.49522739, 0.03448508, 0.90932947, 0.2588541 ])}), Mu({'diffusion': array([0.66255603, 0.3117799 , 0.52011601, 0.54675561, 0.18493597,
0.96958767, 0.77515531, 0.93950499, 0.89483787])}), Mu({'diffusion': array([0.59794019, 0.92188205, 0.08858365, 0.19606326, 0.04532277,
0.3253978 , 0.38873842, 0.2714219 , 0.82875464])}), Mu({'diffusion': array([0.35681765, 0.28100642, 0.54274181, 0.14101013, 0.80221676,
0.07464319, 0.98688825, 0.77226754, 0.19879581])}), Mu({'diffusion': array([0.00562156, 0.81547988, 0.70688666, 0.72903427, 0.77129322,
0.07413725, 0.35852988, 0.11595747, 0.86311712])}), Mu({'diffusion': array([0.6233358 , 0.33096494, 0.06365199, 0.31105122, 0.3252508 ,
0.72963322, 0.63759372, 0.88722402, 0.4722677 ])}), Mu({'diffusion': array([0.11968229, 0.71327346, 0.76080897, 0.56132107, 0.77099008,
0.49384622, 0.52278056, 0.42759826, 0.02551658])}), Mu({'diffusion': array([0.10798064, 0.03152604, 0.63644677, 0.31442455, 0.50861983,
0.90757572, 0.2493673 , 0.41044188, 0.75557558])}), Mu({'diffusion': array([0.22887529, 0.07707221, 0.28982248, 0.16130517, 0.92970468,
0.80813957, 0.63344042, 0.87147344, 0.80369171])}), Mu({'diffusion': array([0.1866514 , 0.89256974, 0.53938831, 0.80745941, 0.89610169,
0.31807167, 0.11014092, 0.22801237, 0.42716508])}), Mu({'diffusion': array([0.81803296, 0.86074451, 0.00705144, 0.51079623, 0.41746926,
0.2221856 , 0.11995338, 0.33768141, 0.94291541])}), Mu({'diffusion': array([0.32327061, 0.51883874, 0.70304866, 0.36369324, 0.9717849 ,
0.96245105, 0.25185712, 0.49729878, 0.30094822])}), Mu({'diffusion': array([0.28491201, 0.03698326, 0.60960338, 0.50272876, 0.0515736 ,
0.2787186 , 0.90827506, 0.23963793, 0.14498038])}), Mu({'diffusion': array([0.48950382, 0.98565189, 0.24213107, 0.67216833, 0.76164345,
0.23771378, 0.72824353, 0.36784635, 0.6323426 ])}), Mu({'diffusion': array([0.63356636, 0.53582111, 0.09038074, 0.83531897, 0.32084799,
0.18659986, 0.04087106, 0.59093385, 0.67759661])}), Mu({'diffusion': array([0.01668617, 0.51214185, 0.22657313, 0.64520827, 0.17444899,
0.69096864, 0.38679667, 0.93673632, 0.13760719])}), Mu({'diffusion': array([0.34113224, 0.11356217, 0.92470115, 0.87735162, 0.25801583,
0.66001805, 0.81724048, 0.55524529, 0.52969761])}), Mu({'diffusion': array([0.24192811, 0.09319346, 0.89722604, 0.90042802, 0.63313815,
0.33909589, 0.34927465, 0.72598308, 0.89712055])}), Mu({'diffusion': array([0.88709772, 0.77989756, 0.64206744, 0.08423155, 0.16171255,
0.89856433, 0.60646842, 0.00929613, 0.1015614 ])}), Mu({'diffusion': array([0.66353542, 0.00516108, 0.16089197, 0.54877892, 0.69192601,
0.65199606, 0.22434688, 0.712208 , 0.23732536])}), Mu({'diffusion': array([0.32546716, 0.74651676, 0.64966794, 0.84923849, 0.65764713,
0.56835177, 0.0937654 , 0.36777903, 0.26527585])})]
```

Then we `solve`

the full-order model
for all `parameter values`

in the training set and accumulate all
solution vectors in a single `VectorArray`

using its
`append`

method. But first
we need to create an empty `VectorArray`

to which the solutions can
be appended. New `VectorArrays`

in pyMOR are always created by a
corresponding `VectorSpace`

. Empty arrays are created using the
`empty`

method. But what
is the right `VectorSpace`

? All `solutions`

of a `Model`

belong to its `solution_space`

,
so we write:

```
U = fom.solution_space.empty()
for mu in training_set:
U.append(fom.solve(mu))
```

Note that `fom.solve`

returns a `VectorArray`

containing a single vector.
This entire array (of one vector) is then appended to the `U`

array.
pyMOR has no notion of single vectors, we only speak of `VectorArrays`

.

What exactly is a `VectorSpace`

? A `VectorSpace`

in pyMOR holds all
information necessary to build `VectorArrays`

containing vector objects
of a certain type. In our case we have

```
fom.solution_space
```

```
NumpyVectorSpace(20201, id='STATE')
```

which means that the created `VectorArrays`

will internally hold
`NumPy arrays`

for data storage. The number is the dimension of the
vector. We have here a `NumpyVectorSpace`

because pyMORâ€™s builtin
discretizations are built around the NumPy/SciPy stack. If `fom`

would
represent a `Model`

living in an external PDE solver, we would have
a different type of `VectorSpace`

which, for instance, might hold a
reference to a discrete functions space object inside the PDE solver
instead of the dimension.

After appending all solutions vectors to `U`

, we can verify that `U`

now really contains 25 vectors:

```
len(U)
```

```
25
```

Note that appending one `VectorArray`

`V`

to another array `U`

will append copies of the vectors in `V`

to `U`

. So
modifying `U`

will not affect `V`

.

Letâ€™s look at the solution snapshots we have just computed using
the `visualize`

method of `fom`

.
A `VectorArray`

containing multiple vectors is visualized as a
time series:

```
fom.visualize(U)
```