PDE complex numbers

735
views
1
9 months ago by
Hello,

Can anyone confirm that complex numbers are still not supported on fenics?

Is there a good example on how to solve the wave equation in 2D in the frequency domain (or something similar) somewhere?

I've tried to follow this https://fenicsproject.org/qa/6775/gmres-convergence-stronly-depends-on-frequency, but I get problems (TypeError: unsupported operand type(s) for *: 'FunctionSpace' and 'FunctionSpace') with the following:

V = FunctionSpace(mesh, 'Lagrange', 1)
V2 = V * V
Community: FEniCS Project

4
9 months ago by

To solve the Helmholtz equation with 1st order absorbing boundary condition, for example,  you can use something like:

Vr = FiniteElement('P', mesh.ufl_cell(), degree)
Vi = FiniteElement('P', mesh.ufl_cell(), degree)
Vcomplex = Vr * Vi
V = FunctionSpace(mesh, Vcomplex)

#Real and Imaginary parts of the trial and test functions
ur, ui = TrialFunctions(V)
vr, vi = TestFunctions(V)

+ k*(ur*vi*ds) - k*(ui*vr*ds)

# Source fr - real fi - imaginary
L = (fr*vr + fi*vi)*dx

A = assemble(a)
b = assemble(L)

u =  Function(V)

solve(A,u.array(),b)

# Split u in real and imaginary parts
ur , ui = u.split()

However, If you want to get the complex matrix, it's just a matter of finding the real and imaginary dofs:

# dofs Real Part
dofs_real = V.sub(0).dofmap().dofs()
# dofs Imaginary Part
dofs_imag = V.sub(1).dofmap().dofs()

To correctly extract the real and imaginary matrices, in one of scipy sparse matrix classes:

# Slice Matrix - Real Part
A_real = A[dofs_real,:]
A_real = A_real[:,dofs_real]
b_real = b [dofs_real]

# Slice Matrix - Real Part
A_imag = A[dofs_imag,:]
A_imag = A_imag[:,dofs_real]
b_imag = b [dofs_imag]

Does k represent the velocity of the wave?

Where there are no mixed terms ur, vi and ui, vr? http://pe.org.pl/articles/2017/3/40.pdf at equation (10) they say that the system to solve for complex values is [A_r -A_i; Ai A_r][u_r, u_i] = [v_r, v_i]
written 9 months ago by caterinabig

Hello,
k is the wavenumber.

In the code above, vte mixed terms are due to the abc (absorbing boundary conditions):
a = ....
k*(ur*vi*ds) - k*(ui*vr*ds)
​

Actually, this system is correct only if the imaginary dofs are numbered after the real ones.
written 9 months ago by Igor Baratta
OK thanks!

I've tried to write down in a more extensive form what my problem is here: https://www.allanswered.com/post/wlkmr/wave-equation-in-2d-solved-in-the-frequency-domain/

I've tried to change the function space definition with what you suggested, and I get the exact same result..

Any chance you can spot what I'm doing wrong?

Thanks!
written 9 months ago by caterinabig
0
9 months ago by
You might want to use
mixed_element = MixedElement(V.ufl_element(),V.ufl_element())
V2 = FunctionSpace(mesh,mixed_element)​

I believe fenics does not support complex numbers so you will have to solve for the real and imaginary part equations separately.

0
9 months ago by
I've tried
V2 = VectorFunctionSpace(mesh, 'CG', pdim, dim=2)
(u_r, u_i) = TrialFunction(V2)
(phi_r, phi_i) = TestFunction(V2)​

...

u = Function(V2)

A_tot = A_block(u_r,phi_r,omega_) - A_block(u_i,phi_i,omega_) + A_block(u_r,phi_i,omega_) + A_block(u_i,phi_r,omega_)
b_tot = b_block(source_termR,phi_r) + b_block(source_termI,phi_i)

A = assemble(A_tot)
b = assemble(b_tot)
[bc.apply(A,b) for bc in bcs]
solve(A, u.vector(), b)

(u_real, u_imag) = u.split(True)



but I don't get the expected result.. and I can't tell if it is a bug somewhere else or if this "split" doesn't work properly.

Ideally I'd like to do something like in this paper  http://pe.org.pl/articles/2017/3/40.pdf

Does anyone has an example ???
This should work:
V = VectorFunctionSpace(mesh, 'P', pdim, dim=2)

(u_r, u_i) = TrialFunctions(V)
(phi_r, phi_i) = TestFunctions(V)​
written 9 months ago by Adam Janecka
This is exactly what I posted above, except for the 'P', which I think is equivalent even though not reported on the official manual. Someone already asked about it here: https://fenicsproject.org/qa/10141/which-element-family-is-p

Anyway, in the example https://www.allanswered.com/post/kzqk/auto-adaptive-poisson-equation/ there is no mention on how to recover the real and imaginary parts once we solve for u=Function(V). Do you know how to do it? Is split good?
written 9 months ago by caterinabig
1
Well, it is not exactly the same: TrialFunction vs. TrialFunctionsTestFunction vs. TestFunctions, there should be a plural. The corresponding help entries are:
TrialFunctions(V)
UFL value: Create a TrialFunction in a mixed space, and return a
tuple with the function components corresponding to the subelements.​
and
TestFunctions(V)
UFL value: Create a TestFunction in a mixed space, and return a
tuple with the function components corresponding to the subelements.​

To recover the real and imaginary parts, (u_r, u_i) = u.split() should work just fine.
written 9 months ago by Adam Janecka