Solving linear systems

Consider a simple operator equation of the form

\[A\phi = f\]

Here, \(f\) and \(\phi\) are grid functions and \(A\) is a boundary integral operator.

The weak form of the above equation is given by

\[\langle Au, v\rangle = <f, v>\]

for all \(v\) in the dual space to the range space.

The standard way to solve this system in BEM++ consists of the following steps.

  1. Compute the weak form of \(A\) using discrete_op = A.weak_form().

  2. Compute the projection of the right-hand side \(f\) onto the dual space by p = f.projections(dual_to_range), where dual_to_range is the name of the dual space.

  3. Solve using an iterative solver from Scipy (e.g.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 BEM++ implements within the bempp.api.linalg module 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. Hence, we can simply call:

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

The functions in the 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. BEM++ has its own discretization mode that automatically applies the inverse of the mass matrix by computing its LU factorization. In this case we would use:

from scipy.sparse.linalg import gmres
discrete_op = A.strong_form()
coeffs = f.coefficients
x, info = gmres(A, coeffs)

Please note that instead of the projections on the right-hand side we now need to use the basis coefficients themselves. This typically converges faster than the standard weak form discretization.

If the linear system should be solved by a dense LU decomposition we need a representation of the Galerkin discretization as Numpy array. This is achieved via the command numpy_mat = bempp.api.as_matrix(A). If H-Matrix discretization was enabled this operation can be very slow. It should only be used for very small problems.

Function and class reference

bempp.api.linalg.iterative_solvers.cg(A, b, tol=1e-05, maxiter=None, use_strong_form=False, return_residuals=False, return_iteration_count=False)

Interface to the scipy.sparse.linalg.cg function.

This function behaves like the scipy.sparse.linalg.cg function. But instead of a linear operator and a vector b it takes a boundary operator and a grid function. The result is returned as a grid function in the correct space.

bempp.api.linalg.iterative_solvers.gmres(A, b, tol=1e-05, restart=None, maxiter=None, use_strong_form=False, return_residuals=False, return_iteration_count=False)

Interface to the scipy.sparse.linalg.gmres function.

This function behaves like the scipy.sparse.linalg.gmres function. But instead of a linear operator and a vector b it takes a boundary operator and a grid function or a blocked operator and a list of grid functions. The result is returned as a grid function or as a list of grid functions in the correct spaces.

bempp.api.linalg.direct_solvers.lu(A, b, lu_factor=None)

Simple direct solver interface.

This function takes an operator and a grid function, converts the operator into a dense matrix and solves the system via LU decomposition. The result is again returned as a grid function.

Parameters:
  • A (bempp.api.BoundaryOperator) – The left-hand side boundary operator
  • b (bempp.api.GridFunction) – The right-hand side grid function
  • lu_decomp (tuple) – Optionally pass the tuple (lu, piv) obtained by the scipy method scipy.linalg.lu_factor
bempp.api.as_matrix(operator)

Convert a discrete operator into a dense matrix.

Parameters:operator (scipy.sparse.linalg.interface.LinearOperator) – The linear operator to be converted into a dense matrix.

Notes

Note that this function may be slow depending on how the original discrete operator was stored. In the case of a dense assembly simple the underlying NumPy matrix is returned. Otherwise, the operator needs to be converted to an array, which can take a long time.