GSOC_review.md

Complex Number support in FEniCS

This blog post contains the overview of the work done for the Google Summer of Code 2018.

The goal of this project was to open up the possibility of the solution of large-scale complex-valued PDEs using FEniCS. To achieve this goal, contributions of varying degree have been made on some FEniCS components, namely DOLFIN-X, FFC-X and UFL.

To better understand these contributions, a chronological list (oldest first) of the Pull Requests (PR) created during the project will be provided below. A short summary will accompany each PR, but more information can be found on their respective links.

Pull Requests List

DOLFIN-X

• PR #72 - Make la::PETScMatrix and la::PETScVector compatible with PETSc-complex - merged on May 10
  • This was my first Pull Request during the GSoC. The purpose of this PR was to make the classes PETScMatrix and PETScVector from the linear algebra (la) module compatible with PETSc compiled with complex number support.
• PR #81 - Add some complex compatibility changes - merged on May 16
  • This PR can be considered an extension of the previous one, with modifications on the class VectorSpaceBasis and addition of the MPI definition for complex double (MPI_DOUBLE_COMPLEX).
• PR #87 - Add complex suppport to la::SLEPcEigenSolver - merged on May 21
  • Although SLEPc separates the real and imaginary part of eigenvalues and eigenvectors, when PETSc is compiled with complex support this behavior changes. In this PR, changes were made to create a common interface between real and complex modes.
• PR #93 - Update python la to handle complex numbers - merged on May 29
  • With the previous three PRs, the linear algebra module, within the C ++ layer, was already compatible with complex PETSc. In this PR, the la python interface was also updated to handle complex numbers.
• PR #96 - Define ufc_scalar_t for a complex UFC interface - merged on Jun 6
  • This PR should be considered together with FFC-x PR #31 (more information on the FFC-X section).
• PR #104 - Update the function module to handle complex numbers - merged on Jun 26
  • The changes that have been made in this PR allow the creation of complex functions (functions, expressions and constants). In addition, the literal j, representing the imaginary unit, was defined to facilitate the creation of complex expressions using C ++ syntax.
• PR #112 - Reorder complex values - merged on Jul 10
  • The purpose of this PR is to allow the reordering of complex values according to explicit global indices. Later, these changes were used to represent complex values in XDMF.
• PR #120 - Fix some tests and demos to work in both real and complex modes - merged on Jul 10
  • Basically, the changes of this PR included the exchange of the operator * with inner in the forms expressed in the tests, to reflect the changes in UFL, since the emergence of complex support. Since now the operator inner takes the conjugate of the second term and and the former does not.
• PR #127 - More changes in demos and tests to work in complex mode - merged on Jul 11
  • This PR can be considered a sequence of the previous one, with the addition of some demos.
• PR #128 - Update assembler to handle complex-valued equations - merged on Jul 13
  • In this PR, the goal was to allow assembling of complex-valued forms and boundary conditions. A set of unit tests was also added.
• PR #134 - Add Complex support to XDMFFile - merged on Jul 27
  • This PR added support for writing and reading complex functions using the class XDMFFile. The real and imaginary components are split into two different datasets to allow visualization using some visualization application, such as paraview.
• PR #141 - Uptade XDMF tests to cover complex mode - merged on Aug 4
  • The aim, in this PR, was to increase the coverage of the current xdmf tests to treat the complex case.
• PR #144 - Get dolfinx to compile with petsc-complex (remaining changes) - merged on Aug 4
  • This PR involved the remaining modifications to get dolfinx to compile with PETSc-complex, in preparation to run nightly tests for the complex mode on master branch.
• PR #146 - Set nightly builds for complex mode - merged on Aug 4
  • To test all changes to the code base considering the complex mode, CircleCI nightly builds were set. The complex build, was configured to run every day at 05:00 am UTC.
• PR #149 - Update the scalar type according to the mode (real or complex) - merged on Aug 4
  • In this PR, the scalar type parameter, to be passed to the form compiler, was update according to the mode. For example, if PETSc is compiled with complex support we pass a complex scalar to the form compiler.
• PR #153 - Set nightly build branch to master - merged on Aug 13
  • This is a one-line PR with small fix to set the nightly build branch to run only on the master branch.

FFC-X

• PR #26 - Add #include<complex.h> to the include list - closed on May 14
  • This PR was closed. And its idea of including suitable headers for the complex mode in the generated code was incorporated on PR #31.
• PR #31 - Create interface for complex scalars - merged on Jun 21
  • The aim of this PR was to create a more flexible interface for the generated code, by using the ufc_scalar_t datatype that could be real or complex depending on the mode. This should apply to tabulate_tensor and transform_values functions. Also, the interface with the TSFC representation was updated, allowing the generation of code from complex-valued forms.
• PR #51 - Add scalar_type to ffc default parameters - merged on Aug 04
  • This PR adds the scalar_type parameter to the list of ffc default parameters. By default scalar_type is set to double.
• PR #47 - Add scalar_type to command line options - opened on Jul 13
  • Due to a disagreement on how the command line flags for the complex mode should be specified, this PR remains open.

UFL

• PR #95 - Extend UFL to the complex domain - merged on Jul 23
  • Fundamentally, the changes that allow the extension to complex values came from the Firedrake/UFL branch, with minor fixes to make it compatible with dolfin, and some improvements on the documentation.

Demonstration

For demonstration purposes we will consider the problem of an acoustic plane-wave scattering by an infinitely long cylinder. A detailed description of the continuous variational problem we are interested in can be found here.

Compile a complex problem described by a UFL file with FFC-X

The first step is to clone and install FFC-X from source (please refer to https://github.com/FEniCS/ffcx). To use the complex mode, some additional dependencies from the TSFC representation backend are needed:

pip3 install six singledispatch pulp networkx
pip3 install git+https://github.com/blechta/COFFEE
pip3 install git+https://github.com/blechta/FInAT
pip3 install git+https://github.com/blechta/tsfc

The UFL specification serves as the input of a form compiler. The UFL file (Helmholtz.ufl) for this demonstrantion can be found here or below:

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)

f = Coefficient(element)
g = Coefficient(element)
k = Coefficient(element)

a = -inner(grad(u), grad(v))*dx + k**2*inner(u, v)*dx - 1j*k*inner(u, v)*ds
L = inner(f, v)*dx + inner(g, v)*ds

Now to compile this form, the following command line can be used:

python3 -m ffc -l dolfin -f 'scalar_type' 'double complex' -r 'tsfc' Helmholtz.ufl

Solve a complex-valued PDE with DOLFIN-X

Building dolfin-x from source is somewhat more involving than FFC-X. It is therefore recommended to use an image that already contains petsc installed with support for complex numbers.

A (non-official) docker image containing all dependencies, can be pulled using: docker pull quay.io/igorbaratta/dolfinx-complex

And then just donwload the python code describing our model problem. The complete code can be found here and below:

from dolfin import (Constant, DirichletBC, Expression, FacetNormal, Function,
                    FiniteElement, FunctionSpace, TrialFunction, TestFunction,
                    dx, ds, dot, grad, inner)
from dolfin.io import XDMFFile
from dolfin.fem.assembling import assemble_system
from dolfin.la import PETScKrylovSolver, PETScOptions, PETScVector
import geometry
import numpy as np

# Generate mesh and subdomains
mesh, boundaries = geometry.circleHole(h=1e-2)

# Problem Parameters
k = 5*np.pi         # Wavenumber
theta = np.pi       # Wave direction
degree = 6          # Polynomial degree

# Define finite element and function space
P = FiniteElement("Lagrange", mesh.ufl_cell(), degree)
V = FunctionSpace(mesh, P)

# Compute incoming wave and boundary source term
n = FacetNormal(mesh)
ex = Expression("exp(-j*k*(cos(theta)*x[0]+sin(theta)*x[1]))",
                k=k, theta=theta, degree=10, domain=mesh.ufl_domain())
g = dot(grad(ex), n) + 1j*k*ex


# Apply Dirichlet Boundary condition (on cylinder surface)
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundaries, 1)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
a = -inner(grad(u), grad(v))*dx + k**2*inner(u, v)*dx - 1j*k*inner(u, v)*ds
L = inner(g, v)*ds

# Assemble system
A, b = assemble_system(a, L, [bc])

# Solve linear system
solver = PETScKrylovSolver(mesh.mpi_comm())
PETScOptions.set("ksp_type", "preonly")
PETScOptions.set("pc_type", "lu")
solver.set_from_options()
solver.set_operator(A)
x = PETScVector()
n_iter = solver.solve(x, b)

# Write computed solution
u = Function(V)
u.vector().vec().array = x.vec().array
with XDMFFile(mesh.mpi_comm(), "scattering.xdmf") as xdmf:
    xdmf.write(u)

The Scattering.py example can be run using:

cd path/to/file/
python3 Scattering.py

Then, simply open the scattering.xdmf on paraview to check the solution.

Real part of the total field:

Imaginary part of the total field:

Future Work

• Improve the support for complex-valued equations in FFC-X - Currently the compilation of complex problems is done exclusively with a representation (TSFC) that requires several external dependencies. One possible way forward to achieve this improvement could be to add support for complex numbers in the UFLACS representation.

• Increase the number and scope of complex tests, since they now run separately from real-mode tests on CircleCI.