 # Solving complex-valued PDE

Hi,
I’m trying to solve the helmholtz equation in its form Laplace§+k^2*p = 0, where p is a complex-valued quantity. I have a 3D Box mesh and following boundary conditions.
boundary_conditions = {
0: {“Neumann”: a_val * v_s}, # 0 left boundary
1: {“Neumann”: a_val * v_s}, # 1 lower boundary
2: {“Robin”: (a_val * Y_right, -1 * Z_right * v_s)}, # 2 right boundary
3: {“Neumann”: a_val * v_s}, # 3 upper boundary
4: {“Neumann”: a_val * v_s}, # 4 front boundary
5: {“Neumann”: a_val * v_s}, # 5 back boundary
}
where Y_right and Z_right are real valued quantities, a_val and the list v_s has a complex values. I’m then collecting my Neumann and Robin boundary conditions in lists and evalaute them.
bc = []
integrals_n = [] # list for NBC
integrals_ra = [] # list for RBC
integrals_rl = [] # list for RBC
j = complex(0, 1)
for i in boundary_conditions:
# collect Neumann integrals
if “Neumann” in boundary_conditions[i]:
if boundary_conditions[i][“Neumann”] != 0:
g = boundary_conditions[i][“Neumann”]
integrals_n.append(df.inner(g, v) * df.ds(i))
# Collect Robin integrals, can be split up in two lists
if “Robin” in boundary_conditions[i]:
if boundary_conditions[i][“Robin”]:
r, s = boundary_conditions[i][“Robin”]
integrals_ra.append(r * p * v * df.ds(i))
integrals_rl.append(r * s * v * df.ds(i))
Afterwards i’m building my left and right side and solve it.
left_side = -1 * (df.inner(df.grad§, df.grad(v)) - (k ** 2) * v * p) * df.dx - sum(integrals_ra)
right_side = df.Constant(0) * v * df.dx + sum(integrals_n) - sum(integrals_rl)

u = df.Function(function_space)

df.solve(left_side == right_side, u, bc)
Since p is a complex valued quantity, i’m expecting a complex solution, but i’m only getting the following error message:
ufl.log.UFLException: Unexpected complex value in real expression.
I’m not able to figure out why this is happening. I was trying to do it like in here:

But none of this is working and i wasn’t able to find anything on this forum, so i’m hoping anyone of you can help me.
There are also only Neumann and Robin boundary conditions, therefore i’m not sure how to handle it properly.
Any little help is highly appreciated!

Which version of fenics are you using?

1 Like

I’m using FEniCS version 2019.1.0, sorry I forgot to mention it.

The example you’re following comes from dolfin-x, the development version of dolfin. dolfin 2019.1 does not support complex numbers. You’ll need to build dolfin-x from source.

(an alternative of course is to split the real and imaginary parts and solve the simultaneous system of 2 real equations using a MixedFunctionSpace)

1 Like

Thank you very much, is there an example of how to calculate the real and imaginary parts using a MixedFunctionSpace?

You’ll need to use your own maths to split the equations. But there are examples of MixedFunctionSpaces

e.g. Stokes equations (solving for u and p),
https://fenicsproject.org/docs/dolfin/2019.1.0/python/demos/stokes-iterative/demo_stokes-iterative.py.html
/usr/share/dolfin/demo-python/documented/stokes-iterative/

also
/usr/share/dolfin/demo-python/undocumented/mixed-poisson-sphere/demo_mixed-poisson-sphere.py

1 Like

Thank you very much that really helped me out a lot, i’m currently trying to figure out how to build dolfin-x from source, but i will try it with the MixedFunctionSpace.

Hey Julian,

I work with the same PDE and I use Fenics 2019.1.0 as well.

We use the `FunctionSpace` this way

``````P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

# # Define 2D-Real function space to substitute for Complex space
V = FunctionSpace(mesh, P1 * P1)   # NOTE: V.sub(0) and V.sub(1) are the subspaces
``````

and take the Trial and test functions as tuples from V i.e

``````(u_re, u_im) = TrialFunctions(V)
(v_re, v_im) = TestFunctions(V)
``````

All the operations over ufl forms can then be simplified by defining a Complex class to handle complex objects.

Cheers, Hope this helps!

1 Like

Thank you very much, I really appreciate it.
How have you implemented your boundary conditions on this problem?

Dirichlet boundary conditions are imposed on each subspace :

``````bcs = [DirichletBC(V.sub(0), u0_re, boundary_markers, 1),
DirichletBC(V.sub(1), u0_im, boundary_markers, 1)]
``````

The Neumann conditions can be easily applied by dividing the integrals for real and imaginary parts.

``````g_re * v_re * ds(1) + g_im * v_im * ds(1)
``````

PS: This is a based on a simple (yet, neat) trick of representing the complex space of test functions by `v_re + i * 0` and `0 + i * v_im` functions.

1 Like

Wow thank you very very much, it really helped me out a lot!
Can you post your weak formulation so I can compare it?
I’m not sure if I have done it right.

Maybe one question here:
I am having a similar problem and have my FunctionSpace set up like in your examples above.
However, I am using SLEPcEigenSolver and it looks like the complex part of the eigenvalue is always zero. (I’ll have to sit for an MWE…)

Has anyone here ever successfully obtained complex eigenvalues from a real-valued problem with `dolfin.SLEPcEigenSolver`?