How to use complex numbers + iterative solvers?

Question

I would like to solve large wave propagation problems involving absorbing media. For this reason i need (1) iterative solvers and (2) complex valued functions.

The use of iterative solvers is demonstrated in the official tutorial example "d2_p2D.py" and while complex numbers are not supported by FEniCS directly, the limitation can be circumvented by using a two-valued FunctionSpace to which the (manually separated) real and imaginary parts of the equation(s) are assigned.

Hence in isolation, I am able to use both iterative solvers and complex valued functions. However, it turns out that when the complex value workaround is applied, the stiffness matrix is no longer positive definite (not even positive semi definite), and the iterative solvers thus no longer converge.

Is there anyone who can enlighten me to why this happens and/or have an idea of how to solve the problem?

Below a MWE crafted from the "d2_p2D.py" tutorial example is attached. The example is changed to deal with a complex field where the real/imaginary components "bend" in opposite directions. After the modification, only the direct solvers work.

from dolfin import *
import numpy as np

# region Helper functions

def is_pd(x):
    return np.all(np.linalg.eigvals(x) > 0)


def is_psd(x, tol=1e-5):
    E, V = np.linalg.eigh(x)
    return np.all(E > -tol)

# endregion

# Create mesh and define function space
mesh = UnitSquareMesh(15, 15)
V = FunctionSpace(mesh, 'Lagrange', 1)
V2 = V * V
# Define boundary conditions
u0 = Expression(('1 + x[0]*x[0] + 2*x[1]*x[1]', '-1 - x[0]*x[0] - 2*x[1]*x[1]'))
bc = DirichletBC(V2, u0, lambda x, ob: ob)
# Create test and trail functions
u_r, u_i = TrialFunction(V2)
v_r, v_i = TestFunction(V2)
# Define variational problem
a = + (inner(nabla_grad(u_r), nabla_grad(v_r)) - inner(nabla_grad(u_i), nabla_grad(v_i))) * dx \
    + (inner(nabla_grad(u_r), nabla_grad(v_i)) + inner(nabla_grad(u_i), nabla_grad(v_r))) * dx
L = (v_r * Constant(-6.0) + v_i * Constant(6.0)) * dx

# Compute solution
u = Function(V2)
# solve(a == L, u, bc)  # Direct: Working
solve(a == L, u, bc, solver_parameters={'linear_solver': 'cg', 'preconditioner': 'ilu'})  # Iterative: NOT working
u_r, u_i = u.split(True)  # Split in real/imag parts

# Inspect stiffness matrix
A = assemble(a).array()
print('A is {}'.format('PD' if is_pd(A) else ('PSD' if is_psd(A) else 'not PSD')))

# Plot solution
plot(u_r)
plot(u_i)
interactive()

Kent-Andre Mardal · Answer 1 · Feb 4, 2016

It seems symmetric but indefinite. Hence you should use the Minimal Residual method.
It is called tfqmr in PETSc. However, you may want to avoid using the above form in
the construction of the preconditioner and rather use something like:

p = + (inner(nabla_grad(u_r), nabla_grad(v_r)) + inner(nabla_grad(u_i), nabla_grad(v_i))) * dx

Look at the stokes-iterative demo for how to use a different form in the construction
of the preconditioner. MinRes requires a positive preconditioner.

How to use complex numbers + iterative solvers?

1 Answer