### Vector-valued PDE with complex numbers

305
views
1
4 months ago by
Hello,

I've a question about solving a Vector PDE in the Complex plane. The question is related to my previous question https://www.allanswered.com/post/vnebp/pde-complex-numbers/. In that case I was solving a scalar PDE though.

For a scalar problem I'm doing the following:
Vr = FiniteElement('P', mesh.ufl_cell(), pdim)
Vi = FiniteElement('P', mesh.ufl_cell(), pdim)
Vcomplex = Vr * Vi
V = FunctionSpace(mesh, Vcomplex)
u_r, u_i = TrialFunctions(V)
phi_r, phi_i = TestFunctions(V)

What should I do if instead I'd like to have a space that has 2 dimensions?, i.e. I'd like to have both u_r and u_i as vectors of dimensions 2. I've tried replacing this line
V = VectorFunctionSpace(mesh, Vcomplex, dim=2)​
but I get the following error " TypeError: VectorFunctionSpace() takes at least 3 arguments (3 given)"

Any idea why?

Thank you,
Caterina
Community: FEniCS Project

3
4 months ago by
You should create a MixedElement

from dolfin import *

mesh = UnitSquareMesh(1, 1)
Vr = VectorElement("CG", mesh.ufl_cell(), 1)
Vi = VectorElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, MixedElement([Vr, Vi]))
ur, ui = TrialFunctions(V)
phi_r, phi_i = TestFunctions(V)
​
Thank you!

I still have one question.

I'm trying to solve the same problem in time and in space. In time I'm doing
u1,u1_proj = Function(V),Function(V)
while t+dt < T:
{....
solve(A, u1.vector(), b)
u1_proj.vector()[:] = project(u1, V).vector()

exact_point = Point(np.array((xs,ys)))
u1value = u1_proj(exact_point)
signal_x[tstep] = u1value[0]
signal_y[tstep] = u1value[1]​
...
}
How can I do the same kind of projection with vector valued real and imaginary part? Here is what I've tried:
u_real,u_imag,uR_proj,uI_proj = Function(V),Function(V),Function(V),Function(V)
...
for kk in range(0,N):
{....
solve(A, u.vector(), b)
(u_real, u_imag) = u.split(True)

uR_proj.vector()[:] = project(u_real, V).vector()
uI_proj.vector()[:] = project(u_imag, V).vector()

# Save the solution at the sensors location
# -----------------------------------------
exact_point = Point(np.array((xs,ys)))
uRvalue = uR_proj(exact_point)
uIvalue = uI_proj(exact_point)
F_values_x[k_count,tstep] = uRvalue[0]+1j*uIvalue[0]
F_values_y[k_count,tstep] = uRvalue[1]+1j*uIvalue[1]

...
}​
But I get an error: shape mismatch...
written 4 months ago by caterinabig
1
V is the total function space. You need to extract the real and imaginary subspaces and manipulate accordingly.
written 4 months ago by Nate
I've tried to avoid projection and use directly the splitted real and imaginary solution...

....
solve(A, u.vector(), b)
(u_real, u_imag) = u.split(True)
...
# Save the solution at the sensors location
# -----------------------------------------
exact_point = Point(np.array((xs,ys)))
uRvalue = u_real(exact_point)
uIvalue = u_imag(exact_point)
F_values_x[k_count,tstep] = uRvalue[0]+1j*uIvalue[0]
F_values_y[k_count,tstep] = uRvalue[1]+1j*uIvalue[1]
​

However, I still get different result from the solution in time.

Do you think I'm using the total space V and its solutions correctly?
written 4 months ago by caterinabig
Nate, do you know how things work with Mixed elements?

I mean when I solve and split the problem:
solve(A, u.vector(), b)
(u_real, u_imag) = u.split(True)​

can I be sure that u_real contains indeed the real part of the first component and the second component? How does "it" know that it's not the real and imaginary part of the first component instead???

Also, given V_r and V_i teo VectorElement(...), where are we actually imposing that the solution is a complex vector of type
(u_xR + i u_xI; u_yR + i u_yI)???
written 4 months ago by caterinabig
1
You have to impose that in a variational formulation of your PDE.
written 4 months ago by Jan Blechta
How should I extract the real and imaginary 2D vectors? Is the following correct? I mean, I believe that that u_real containes a vector solution of the real parts of the first and second component (displacement in the x and y directions respectively in my problem). Is this correct?

Here's the code:
....
solve(A, u.vector(), b)
(u_real, u_imag) = u.split(True)
...
# Save the solution at the sensors location
# -----------------------------------------
exact_point = Point(np.array((xs,ys)))
uRvalue = u_real(exact_point)
uIvalue = u_imag(exact_point)
F_values_x[k_count,tstep] = uRvalue[0]+1j*uIvalue[0]
F_values_y[k_count,tstep] = uRvalue[1]+1j*uIvalue[1]
​​

Also, is it enough to write Vr = VectorElement("CG", mesh.ufl_cell(), 1) ? Where do we impose that it's a 2D vector and not a 3D one for example? Is this information hided in the mesh.ufl_cell() ?

Thanks a lot :)

written 4 months ago by caterinabig
1
VectorElement has by default dimension equal to geometric dimension of underlying cell. That can be overridden by providing dim keyword argument, e.g., VectorElement(..., dim=6).

Regarding your code to debug: it is impossible for users of this board to debug incomplete code. I would suggest to opening a new question with complete code which others can run.

Also note that there is no native support for complex numbers in FEniCS yet, although it is on high-priority plan. Complex problems can be dealt with by splitting them into real and imaginary parts but that does not happen automatically. You have to reformulate your PDE from complex plane to  $R^2$R2. In order to learn how to deal with mixed/vector-space problems, complex problems are not the best place to start, I believe. I would suggest to go through documented demos or FEniCS tutorial, see https://fenicsproject.org/documentation/ for resources.
written 4 months ago by Jan Blechta
OK, thanks

I've a code working in R^2 (time domain), which I'm trying to translate to C^2 (frequency domain). I've already done the matching for the scalar case and it works just fine.
written 4 months ago by caterinabig