Solving linear systems

Consider a simple operator equation of the form

\mathsf{A}\phi = f

where \displaystyle f and \displaystyle \phi are grid functions and \displaystyle \mathsf{A} is a boundary integral operator.

The weak form of the above equation is given by

\langle \mathsf{A}u, v\rangle = \langle f, v\rangle

for all \displaystyle v in the dual_to_range space.

The standard way to solve this system in Bempp consists of the following steps:

  1. Compute the weak form of \displaystyle \mathsf{A} using discrete_op = A.weak_form().
  2. Compute the projection of the right-hand side \displaystyle f onto the dual space by = f.projections(dual_to_range), where dual_to_range is the name of the dual space.
  3. Solve using an iterative solver from SciPy (eg. CG or GMRES), potentially together with a suitable preconditioner, that is x, info = gmres(discrete_op, p).
  4. Turn the result x back into a grid function using the command sol_fun = bempp.api.GridFunction(A.domain, coefficients=x).

In order to make these steps shorter, Bempp implements wrappers that look like the corresponding functions in Scipy but automatically do the unpacking of the operators and wrapping of the solution into a grid function. These can be found in the module bempp.api.linalg module. Hence, we can simply call:

from bempp.api.linalg.iterative_solvers import gmres
phi = gmres(A, f)

The functions in the bempp.api.linalg module accept all the usual parameters that are also accepted by the corresponding Scipy routines.

A frequently used form of preconditioning is via the mass matrix. Bempp has its own discretisation mode that automatically applies the inverse of the mass matrix by computing its LU factorisation. In this case we would use:

from bempp.api.linalg.iterative_solvers import gmres
phi = gmres(A, f, use_strong_from=True)

This typically converges faster than the standard weak form discretisation.

If the linear system should be solved by a dense LU decomposition (or another direct method), we need a representation of the Galerkin discretisation as a NumPy array. This is achieved via the command numpy_mat = bempp.api.as_matrix(A). If H-matrix discretisation was enabled this operation can be very slow and should only be used for very small problems.