### PDE complex numbers

631

views

1

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

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

### 3 Answers

4

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)
a = - inner(grad(ur), grad(vr))*dx + (k**2)*inner(ur, vr)*dx\
- inner(grad(ui), grad(vi))*dx + (k**2)*inner(ui, vi)*dx \
+ 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]
```

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
7 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!

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
7 months ago by
caterinabig

0

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

I've tried

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 ???

```
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)
```

See: https://www.allanswered.com/post/kzqk/auto-adaptive-poisson-equation/

written
7 months ago by
Adam Janecka

This is exactly what I posted above, except for the

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

`'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-pAnyway, 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
7 months ago by
caterinabig

1

Well, it is not exactly the same:

To recover the real and imaginary parts,

`TrialFunction`

vs. `TrialFunctions`

& `TestFunction`

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
7 months ago by
Adam Janecka

Please login to add an answer/comment or follow this question.

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]