Technical Overview¶
Three Central Classes¶
From a bird’s eye perspective, pyMOR is a collection of generic algorithms operating on objects of the following types:
VectorArrays
Vector arrays are ordered collections of vectors. Each vector of the array must be of the same
dimension
. Vectors can becopied
to a new array,appended
to an existing array ordeleted
from the array. Basic linear algebra operations can be performed on the vectors of the array: vectors can bescaled
in-place, the BLASaxpy
operation is supported andinner products
between vectors can be formed. Linear combinations of vectors can be formed using thelincomb
method. Moreover, various norms can be computed and selecteddofs
of the vectors can be extracted forempirical interpolation
. To act on subsets of vectors of an array, arrays can beindexed
with an integer, a list of integers or a slice, in each case returning a newVectorArray
which acts as a modifiable view onto the respective vectors in the original array. As a convenience, many of Python’s math operators are implemented in terms of the interface methods.Note that there is not the notion of a single vector in pyMOR. The main reason for this design choice is to take advantage of vectorized implementations like
NumpyVectorArray
which internally store the vectors as two-dimensionalNumPy
arrays. As an example, the application of a linear matrix based operator to an array via theapply
method boils down to a call toNumPy
’s optimizeddot
method. If there were only lists of vectors in pyMOR, the above matrix-matrix multiplication would have to be expressed by a loop of matrix-vector multiplications. However, when working with external solvers, vector arrays will often be given as lists of individual vector objects. For this use-case we provideListVectorArray
, aVectorArray
based on a Python list of vectors.Associated to each vector array is a
VectorSpace
which acts as a factory for new arrays of a given type. New vector arrays can be created using thezeros
andempty
methods. To wrap the raw objects of the underlying linear algebra backend into a newVectorArray
,make_array
is used.The data needed to define a new
VectorSpace
largely depends on the implementation of the underlying backend. ForNumpyVectorSpace
, the only required datum is the dimension of the contained vectors.VectorSpaces
for other backends could, e.g., hold a socket for communication with a specific PDE solver instance. Additionally, eachVectorSpace
has a stringid
, defaulting toNone
, which is used to signify the mathematical identity of the given space.Two arrays in pyMOR are compatible (e.g. can be added) if they are from the same
VectorSpace
. If aVectorArray
is contained in a givenVectorSpace
can be tested with thein
operator.Operators
The main property of operators in pyMOR is that they can be
applied
toVectorArrays
resulting in a newVectorArray
. For this operation to be allowed, the operator’ssource
VectorSpace
must be identical with theVectorSpace
of the given array. The result will be a vector array from therange
space. An operator can belinear
or not. Theapply_inverse
method provides an interface for (linear) solvers.Operators in pyMOR are also used to represent bilinear forms via the
apply2
method. A functional in pyMOR is simply an operator withNumpyVectorSpace(1)
asrange
. Dually, a vector-like operator is an operator withNumpyVectorSpace(1)
assource
. Such vector-like operators are used in pyMOR to representParameter
-dependent vectors such as the initial data of anInstationaryModel
. For linear functionals and vector-like operators, theas_vector
method can be called to obtain a vector representation of the operator as aVectorArray
of length 1.Linear combinations of operators can be formed using a
LincombOperator
. When such a linear combination isassembled
,_assemble_lincomb
is called to ensure that, for instance, linear combinations of operators represented by a matrix lead to a new operator holding the linear combination of the matrices.For many interface methods default implementations are provided which may be overridden with operator-specific code. Base classes for
NumPy
-based operators can be found inpymor.operators.numpy
. Several methods for constructing new operators from existing ones are contained inpymor.operators.constructions
.Models
Models in pyMOR encode the mathematical structure of a given discrete problem by acting as container classes for operators. Each model object has
operators
,products
dictionaries holding theOperators
which appear in the formulation of the discrete problem. The keys in these dictionaries describe the role of the respective operator in the discrete problem.Apart from describing the discrete problem, models also implement algorithms for
solving
the given problem, returningVectorArrays
from thesolution_space
. The solution can becached
, s.t. subsequent solving of the problem for the sameparameter values
reduces to looking up the solution in pyMOR’s cache.While special model classes may be implemented which make use of the specific types of operators they contain (e.g. using some external high-dimensional solver for the problem), it is generally favourable to implement the solution algorithms only through the interfaces provided by the operators contained in the model, as this allows to use the same model class to solve high-dimensional and reduced problems. This has been done for the simple stationary and instationary models found in
pymor.models.basic
.Models can also implement
estimate
andvisualize
methods to estimate the discretization or model reduction error of a computed solution and create graphic representations ofVectorArrays
from thesolution_space
.
Base Classes¶
While VectorArrays
are mutable objects, both Operators
and Models
are immutable in pyMOR: the application of an Operator
to the same
VectorArray
will always lead to the same result, solving a Model
for the same parameter will always produce the same solution array. This has two
main benefits:
If multiple objects/algorithms hold references to the same
Operator
orModel
, none of the objects has to worry that the referenced object changes without their knowledge.The return value of a method of an immutable object only depends on its arguments, allowing reliable
caching
of these return values.
A class can be made immutable in pyMOR by deriving from ImmutableObject
,
which ensures that write access to the object’s attributes is prohibited after
__init__
has been executed. However, note that changes to private attributes
(attributes whose name starts with _
) are still allowed. It lies in the
implementors responsibility to ensure that changes to these attributes do not
affect the outcome of calls to relevant interface methods. As an example, a call
to enable_caching
will set the
objects private __cache_region
attribute, which might affect the speed of a
subsequent solve
call, but not its result.
Of course, in many situations one may wish to change properties of an immutable
object, e.g. the number of timesteps for a given model. This can be
easily achieved using the
with_
method every immutable
object has: a call of the form o.with_(a=x, b=y)
will return a copy of o
in which the attribute a
now has the value x
and the attribute b
the
value y
. It can be generally assumed that calls to
with_
are inexpensive. The
set of allowed arguments can be found in the
with_arguments
attribute.
All immutable classes in pyMOR and most other classes derive from
BasicObject
which, through its meta class, provides several convenience
features for pyMOR. Most notably, every subclass of BasicObject
obtains its
own logger
instance with a class
specific prefix.
Creating Models¶
pyMOR ships a small (and still quite incomplete) framework for creating finite
element or finite volume discretizations based on the NumPy/Scipy software stack. To end up with an appropriate
Model
, one starts by instantiating an analytical problem
which
describes the problem we want to discretize. analytical problems
contain
Functions
which define the analytical data functions associated with the
problem and a DomainDescription
that provides a geometrical definition of the
domain the problem is posed on and associates a boundary type to each part of
its boundary.
To obtain a Model
from an analytical problem
we use a
discretizer
. A discretizer will first mesh the
computational domain by feeding the DomainDescription
into a
domaindiscretizer
which will return the Grid
along with a BoundaryInfo
associating boundary
entities with boundary types. Next, the Grid
, BoundaryInfo
and the various
data functions of the analytical problem
are used to instatiate
finite element
or
finite volume
operators.
Finally these operators are used to instatiate one of the provided
Model
classes.
In pyMOR, analytical problems
, Functions
, DomainDescriptions
,
BoundaryInfos
and Grids
are all immutable, enabling efficient
disk caching
for the resulting Models
, persistent over various
runs of the applications written with pyMOR.
While pyMOR’s internal discretizations are useful for getting started quickly
with model reduction experiments, pyMOR’s main goal is to allow the reduction of
models provided by external solvers. In order to do so, all that needs
to be done is to provide VectorArrays
, Operators
and Models
which
interact appropriately with the solver. pyMOR makes no assumption on how the
communication with the solver is managed. For instance, communication could take
place via a network protocol or job files. In particular it should be stressed
that in general no communication of high-dimensional data between the solver
and pyMOR is necessary: VectorArrays
can merely hold handles to data in the
solver’s memory or some on-disk database. Where possible, we favor, however, a
deep integration of the solver with pyMOR by linking the solver code as a Python
extension module. This allows Python to directly access the solver’s data
structures which can be used to quickly add features to the high-dimensional
code without any recompilation. A minimal example for such an integration using
pybind11 can be found in the
src/pymordemos/minimal_cpp_demo
directory of the pyMOR repository.
Bindings for FEnicS and
NGSolve packages are available in the
bindings.fenics
and
bindings.ngsolve
modules.
The pymor-deal.II repository contains
bindings for deal.II.
Parameters¶
pyMOR classes implement dependence on a parameter by deriving from the
ParametricObject
base class. This class gives each instance a
parameters
attribute describing the
Parameters
the object and its relevant methods (apply
, solve
, evaluate
, etc.)
depend on. Each Parameter
in pyMOR has a name and a fixed dimension, i.e. the
number of scalar components of the Parameter
. Scalar parameters are simply
represented by one-dimensional Parameters
. To assign concrete values to Parameters
the specialized dict-like class Mu
is used.
In particular, it ensures, that all of its values are one-dimensional NumPy arrays
.
The Parameters
of a ParametricObject
are usually automatically derived
as the union of all Parameters
of the objects that are passed to it’s __init__
method.
For instance, an Operator
that implements the L2-product with some user-provided
Function
will automatically inherit all Parameters
of that Function
.
Additional Parameters
can be easily added by setting the
parameters_own
attribute.
Defaults¶
pyMOR offers a convenient mechanism for handling default values such as solver
tolerances, cache sizes, log levels, etc. Each default in pyMOR is the default
value of an optional argument of some function. Such an argument is made a
default by decorating the function with the defaults
decorator:
@defaults('tolerance')
def some_algorithm(x, y, tolerance=1e-5)
...
Default values can be changed by calling set_defaults
.
By calling print_defaults
a summary of all defaults
in pyMOR and their values can be printed. A configuration file with all defaults
can be obtained with write_defaults_to_file
. This file can
then be loaded, either programmatically or automatically by setting the
PYMOR_DEFAULTS
environment variable.
As an additional feature, if None
is passed as value for a function argument
which is a default, its default value is used instead of None
. This allows
writing code of the following form:
def method_called_by_user(U, V, tolerance_for_algorithm=None):
...
algorithm(U, V, tolerance=tolerance_for_algorithm)
...
See the defaults
module for more information.
RuleTables¶
Many algorithms in pyMOR can be seen as transformations acting on trees of
Operators
. One example is the structure-preserving (Petrov-)Galerkin
projection of Operators
performed by the project
method. For instance, a
LincombOperator
is projected by replacing all its children (the Operators
forming the affine decomposition) with projected Operators
.
During development of pyMOR, it turned out that using inheritance for selecting
the action to be taken to project a specific operator (i.e. single dispatch
based on the class of the to-be-projected Operator
) is not sufficiently
flexible. With pyMOR 0.5 we have introduced algorithms which are based on
RuleTables
instead of inheritance. A RuleTable
is simply an ordered list of
rules
, i.e. pairs of conditions to match
with corresponding actions. When a RuleTable
is applied
to an object (e.g. an Operator
),
the action associated with the first matching rule in the table is executed. As
part of the action, the RuleTable
can be easily applied recursively
to the children of the given
object.
This approach has several advantages over an inheritance-based model:
Rules can match based on the class of the object, but also on more general conditions, i.e. the name of the
Operator
or being linear and non-parametric
.The entire mathematical algorithm can be specified in a single file even when the definition of the possible classes the algorithm can be applied to is scattered over various files.
The precedence of rules is directly apparent from the definition of the
RuleTable
.Generic rules (e.g. the projection of a linear non-
parametric
Operator
by simply applying the basis) can be easily scheduled to take precedence over more specific rules.Users can implement or modify
RuleTables
without modification of the classes shipped with pyMOR.
The Reduction Process¶
The reduction process in pyMOR is handled by so called reductors
which take arbitrary Models
and additional data (e.g. the reduced
basis) to create reduced Models
. If proper offline/online
decomposition is achieved by the reductor, the reduced Model
will
not store any high-dimensional data. Note that there is no inherent distinction
between low- and high-dimensional Models
in pyMOR. The only
difference lies in the different types of operators, the Model
contains.
This observation is particularly apparent in the case of the classical
reduced basis method: the operators and functionals of a given discrete problem
are projected onto the reduced basis space whereas the structure of the problem
(i.e. the type of Model
containing the operators) stays the same.
pyMOR reflects this fact by offering with GenericRBReductor
a generic algorithm which can be used to RB-project any model available to pyMOR.
It should be noted however that this reductor is only able to efficiently
offline/online-decompose affinely Parameter
-dependent linear problems.
Non-linear problems or such with no affine Parameter
dependence require
additional techniques such as empirical interpolation
.
If you want to further dive into the inner workings of pyMOR, we
recommend to study the source code of GenericRBReductor
and to step through calls of it’s reduce
method with a Python debugger, such as
ipdb.