API for schemes and operators

Abstract concepts

In the following we give a short mathematical description of the operator concepts used in this module. This can be skipped and the description of the methods further down should still be understandable.

The operator classes represent operators $L\colon V_h\to W_h^*$ where the domain and range space $V_h,W_h$ are discrete function spaces and $W_h^*$ is the dual of $W_h$. We will in general not differentiate between duals and the actual function because the range space is finite dimensional. In the most common example these operators are given by forms $a(v_h,w_h)$ which are linear in the second argument (test function) but can be nonlinear in the first argument (trial function). In this case $L[v_h] := a(v_h,\cdot)\in W_h^*$ or for any $w_h\in W_h\colon <L[v_h],w_h> := a(v_h,w_h)$. This is provided by dune.fem.operator.galerkin for forms defined in ufl.

Mainly operators can be evaluated, i.e., we can compute $w_h:=L[v_h]$ using the __call__ method: operator(v_h,w_h). In addition it is possible to compute the linearization $DL[\bar{v}_h]$ around some function $\bar{v}_h\in V_h$. This is done using the jacobian method on the operator. The linearization is now a bilinear form $\bar{a}(v_h,w_h;\bar{v_h})$.

Note

both the __call__ and the jacobian method can be called using grid functions instead of discrete functions $v_h$ and $\bar{v}_h$, respectively.

Operator without DBC:

In the following we focus on the operators returned by dune.fem.operator.galerkin. The minimal requirements for constructing such an operator is to provide a ufl form where the domain and range space is extracted from the trial and test function. In this case these have to be constructed based on discrete function spaces $V_h,W_h$, respectively, so one can not use the abstract dune.ufl.Space. It is instead possible to provide discrete spaces as extra arguments domainSpace,rangeSpace to override the spaces uses to construct the trial and test functions. Many examples like the one below can be found throughout the tutorial:

u = ufl.TrialFunction(domainSpace)
v = ufl.TestFunction(rangeSpace)
op = dune.fem.operator.galerkin(u**2/2*v*dx)
v_h = rangeSpace.interpolate(...)

# In the following gf is some grid function,
# e.g., a function over the domain space:
op(gf,v_h)
# To linearize firs construct a linear operator
linOp = op.linear()
op.jacobian(gf,linOp) # gf is the function around which to linearize

The call to dune.fem.operator.galerkin could as mentioned above include the domain and ranges spaces to use.

domainSpace: the discrete domain space $V_h$ for the operator.
rangeSpace: the discrete range space $W_h$ for the operator.
linear: return a operator class that can store the linearization $DL$. This contains an underlying (sparse) matrix structure that will depend on the storage argument provided to the spaces.
__call__: evaluate the operator.
jacobian: perform the linearization constructing the matrix inside the provided linear operator. There are two versions for this call: op.jacobian(gf,linOp) and op.jacobian(gf,linOp,b_h). The first version was described above, the second version fills in an addition discrete function b_h. This provides the negative of the constant term in the Taylor expansion $<L[v_h],w_h> = -<b_h,w_h> + <DL[\bar{v_h}]v_h,w_h>$. We provide the negative here but that term is often the right hand side needed to solve the problem. For example if $L[v_h]$ is linear then $DL[0]v_h = b_h$ is the linear system that defines the solution to $L[v_h]=0$.
model:
setCommunicate:
setQuadratureOrders:
gridSizeInterior:

Operator with DBC:

If the operator was constructed not only from a form but including DirichletBC classes which are provided as a tuple/list in the first argument during construction of the operator:

u = ufl.TrialFunction(domainSpace)
v = ufl.TestFunction(rangeSpace)
op = dune.fem.operator.galerkin( [u**2/2*v*dx, dune.ufl.DirichletBC(rangeSpace,g)] )

The resulting operator will have the same methods as before but both the __call__ and the jacobian method will include handing of the Dirichlet boundary conditions:

__call__: this method will as before compute $w_h=L[v_h]$ (where $v_h$ can be a general grid function as pointed out above). But in addition all components in $w_h$ associated with the Dirichlet boundary will be of the form $w_i = v_i - g_i$ where $g_i$ is the given boundary data evaluate at the corresponding degree of freedom.

Note

a more mathematical explanation is to consider all functionals $\lambda$ which are associated with the boundaries, i.e., depend on their argument restricted to the boundary: $\lambda(g_1)=\lambda(g_2)$ if $g_1=g_2$ on the dirichlet boundary. After $w_h=L[v_h]$ the following holds $\lambda(w_h) = \lambda(v_h-g) = \lambda(v_h)-\lambda(g)$ for each of these functionals.
jacobain: the matrix assembled will contain a unit row for each row associated with the Dirichlet constraints. If the version with the right hand side argument is used, the components for the Dirichlet degrees of freedom will contain $g_i$.

Note

this corresponds to the correct linearization of $w_i=u_i-g_i$.

In addition to the methods from the standard operator these operators have additional methods:

setConstraints: in it’s simplest version this method takes a single discrete function $w_h$. This changes the components of $w_h$ at the Dirichlet boundary, i.e., $w_i=g_i$ where the boundary data is taken from the DirichletBC instances. Another version takes a general grid function as a first argument op.setConstrains(v,w_h) which leads to $w_i=v_i$ (i.e. $\lambda(w)=\lambda(v)$ with the boundary degree functionals). Finally one can provide a single number as first argument, e.g., op.setConstraint(0,w_h) which leads to $w_i=0$.
subConstraints: there is only one version of this method op.subConstrains(v,w_h) leading to $w_i = w_i+v_i$.
dirichletIndices: This method returns the indices of dofs on the Dirichlet boundary. By default all indices are returned but a boundary id can be provided to only return a subset of indices.

dirichletBlocks: this method provides indices of the degrees of freedom on the Dirichlet boundary defined by the DirichletBC instances. dblocks = op.dirichletBlocks in blocked form. The indices as returned by the dirichletIndices can be obtained using:

for i,block in enumerate(dblock):         # a block is of size `dimRange` (for Lagrange)
  for j,b in enumerate(block):            # iterate over the block
      # 'b' is the `id` of the boundary as described below
      if b > 0:                           # if b>0 it's a Dirichlet boundary
          # do something

Scheme without DBC:

In this project we use the term scheme to describe a special type of operators with the main additional functionality of solving the system $L[v_h]=0$ (or $L[v_h]=b_h$). A consequence of this is that the domain and the range space need to be identical - Petrov Galerkin method which only require that the dimensions match is not currently available.

All methods from operators are available on schemes. In addition we have

space: this is a convenience method made available since the domain and the range space are identical.
solve: this is the central new method available in schemes. A call scheme.solve(target=v_h) returns a solution to the problem $L[v_h]=0$. This is by default based on a Newton-Krylov solver.
parameters:
parameterHelp:
preconditioning:
setErrorMeasure:
inverseLinearOperator: not sure if this is used

Scheme with DBC:

If a scheme is constructed including boundary conditions then the additional method available on the operator are available on the scheme as well.

Additional

virtualization:

Unclear

dimRange: is available on schemes. The following should be removed or on operator (same as scheme.rangeSpace.dimRange or even scheme.space.dimRange=scheme.domainSpace.dimRange)
gridSizeInterior: is missing on scheme