DOLFIN

Complex-valued formulation returns NaNs

Asked by Nico Schlömer

I would like to convert the curl-curl demo into something that can deal with complex values, so as advised in several answers around here, I split the problem into real an imaginary part and pose the problem in the product space.
As a first attempt, I tried converting the exact demo problem into product space:

======================= *snip* =======================
from dolfin import *

parameters["form_compiler"]["representation"] = "quadrature"

mesh = UnitSphere(8)

PN = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
W = PN * PN

bc = DirichletBC(W, (0, 0, 0, 0, 0, 0), "on_boundary")

rhs = project(Constant([0.0, 0.0, 1.0, 0.0, 0.0, 1.0]), W)

A = Function(W)
v0 = TestFunction(W)
u0 = TrialFunction(W)
ur = u0[0]
ui = u0[1]
vr = v0[0]
vi = v0[1]
a = inner(curl(ur), curl(vr))*dx - inner(curl(ui), curl(vi))*dx \
  + inner(curl(ur), curl(vi))*dx + inner(curl(ui), curl(vr))*dx
L = -inner(rhs[0],vr)*dx + inner(rhs[1], vi)*dx \
  - inner(rhs[0],vi)*dx - inner(rhs[1], vr)*dx

solve(a == L, A, bc)

Ar, Ai = A.split()

P1 = VectorFunctionSpace(mesh, "Lagrange", 1)
v1 = TestFunction(P1)
u1 = TrialFunction(P1)
J = Function(P1)
solve(inner(v1, u1)*dx == dot(v1, curl(Ar))*dx, J)

plot(J, vmin=-1.0, vmax=1.0)

file=File("current_density.pvd")
file << J

interactive()
======================= *snap* =======================

The density equation for J is solved with real(A) as a right hand side, so I expected to get the same result as in the vanilla real-valued demo. Hoever, no data is plotted and the color bar indicates a number of NaNs.

What could be the problem?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Nico Schlömer
Solved:
Last query:
Last reply:
Revision history for this message
Nico Schlömer (nschloe) said : 		#1

I just repeated the same exercise with the Poisson problem and everything seems fine: real and imaginary part are solved orthogonally and spit out the same solution.
Something has to be wrong with this curl-curl code, so I'd highly appreciate comments here.

Revision history for this message
Nico Schlömer (nschloe) said : 		#2

I now broke this down to a small piece of code. Consider

====================== *snip* ======================
from dolfin import *

# Create mesh
mesh = UnitCube(1, 1, 1)

# real-valued rhs
V = FunctionSpace(mesh, "Nedelec 1st kind H(curl)", 1)
f = Constant([0.0, 0.0, 1.0])
v = TestFunction(V)
L = -inner(f, v)*dx
b = assemble(L)
info(b, True)

# complex-valued rhs
W = V * V
f = Constant([0.0, 0.0, 1.0, 0.0, 0.0, 1.0])
fr = f[0]
fi = f[1]
v = TestFunction(W)
vr = v[0]
vi = v[1]
L = -inner(fr, vr)*dx + inner(fi, vi)*dx \
    -inner(fr, vi)*dx - inner(fi, vr)*dx
b = assemble(L)
info(b, True)
====================== *snap* ======================

This constructs the right hand side of a typical curl-curl problem, once for real values and once for complex. The output shows that the complex formulations contains only zeros!
May that have to do with how f is defined? I'm not sure at all how to do this correctly.

Revision history for this message
Nico Schlömer (nschloe) said : 		#3

I think I have it.

What needs to be done is to replace

v = TestFunction(W)
vr = v[0]
vi = v[1]

by

(vr, vi) = TestFunctions(W)

and likewise for the trial function u. Also, the right hand side f can be defined via

fr = Constant([0,0,1])
fi = Constant([0,0,1])

this is different from what's been advised before, but actually works for vector elements as well.