Conversion of complex ufl.algebra to fem.function

variational formulation
1

Dear all,
I am working on dolfinx in complex. The program finds two complex-valued functions (psi1 and psi2), then it computes the following quantity

vel
vel1450×200 11.9 KB

which is a real function. The command to produce u_k is

u_k = hbar * ufl.real(- imaginaryUnit * ( ufl.grad(psi1)[i] * ufl.conj(u1) + ufl.grad(psi2)[i] * ufl.conj(u2) ) )

The variable type for uk is ufl.algebra.Product. I need to convert this quantity into a fem.Function type, in order to visualize it or to put into the RHS of a variational formulation.
Do you have any advice on how to achieve this?

Ps: I try to follow the tutorial on Linear Elasticity, in which a ufl.algebra.sqrt is converted into a fem.Function, however with same procedure fem.Expression produces an error saying

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-93-cee90e875b8a> in <module>()
      1 ConversionSpace = FunctionSpace(domain, ("CG", deg))
----> 2 tmpExpression = fem.Expression(uTildey, ConversionSpace.element.interpolation_points)
      3 

21 frames
/usr/local/lib/python3.7/dist-packages/ffcx/codegeneration/C/format_value.py in format_float(x, precision)
     25             s = "{:.{prec}}".format(float(x), prec=precision)
     26     else:
---> 27         s = repr(float(x))
     28     for r, v in _subs:
     29         s = r.sub(v, s)

TypeError: can't convert complex to float

Thanks in advance

2

We would need a code that reproduces the error, as the error message is quite general and hard to debug without a minimal example.
I would suggest start having a look at: Interpolate expression with complex litteral - #4 by dokken

3

Here’s the code

import numpy as np
from time import process_time

import ufl
from dolfinx import plot, fem
from dolfinx.fem import Function, FunctionSpace, assemble_scalar, form, locate_dofs_geometrical, dirichletbc, form
from dolfinx.fem.petsc import LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, CellType
from ufl import dx, grad, inner, div, conj, real, imag, Dx

from mpi4py import MPI
from petsc4py import PETSc


if np.issubdtype(PETSc.ScalarType, np.complexfloating):
    imaginaryUnit = PETSc.ScalarType(0 + 1j)
    
# Parameters
dt = 1/50;
hbar = 0.03;
inlet_velocity = np.array([1.0, 0.0]);
T = 0.5;
numsteps = T/dt

# Mesh Discretization
dimensions = np.array([10,4]);
nx = 250;
ny = 100;

domain = create_rectangle(MPI.COMM_WORLD, [np.array([0, 0]), dimensions],
                                       [nx, ny], CellType.triangle)
tdim = domain.topology.dim;

# Schrodinger: definition of the functional space 

deg = 2; # linear = 1, parabolic = 2 Lagrangian FE are used
V = FunctionSpace(domain, ("CG", deg))

Psi1 = Function(V)
Psi2 = Function(V)

Psi1.name = "WaveFunction - 1"
Psi2.name = "WaveFunction - 2"

PsiTilde1 = Function(V)
PsiTilde2 = Function(V)

psi1old = Function(V)
psi2old = Function(V)

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)

# Schrodinger: defining the inlet boundary condition

dofs_inlet = locate_dofs_geometrical(V, lambda x: np.isclose(x[0],0))
inletPlaneWave = Function(V)
kvec = inlet_velocity/hbar
inletPlaneWave.interpolate(lambda x: np.exp(imaginaryUnit * (kvec[0]*x[0] + kvec[1]*x[1])))
bc_Inlet = dirichletbc(inletPlaneWave, dofs_inlet)

# Schrodinger: initialization

def initial_condition1(x):
    return 1 + 0j + 0 * x[0] + 1 * x[1]
def initial_condition2(x, c = 0.01 ):
    return c + 0j + 0 * x[0] + 1 * x[1]

psi1old.interpolate(initial_condition1)
psi2old.interpolate(initial_condition2)

theta = 0.5 # Crank-Nicolson is used - theta method is implemented

a = inner(u, v) * dx - imaginaryUnit * dt * hbar / 2 * (1-theta) * inner(grad(u), grad(v)) * dx
L1 = inner(psi1old, v) * dx + imaginaryUnit * dt * hbar / 2 * theta * inner(grad(psi1old), grad(v)) * dx
L2 = inner(psi2old, v) * dx + imaginaryUnit * dt * hbar / 2 * theta * inner(grad(psi2old), grad(v)) * dx

schrodinger_bilinear = form(a)
schrodinger_RHS1 = form(L1)
schrodinger_RHS2 = form(L2)

schrodinger_A = fem.petsc.assemble_matrix(schrodinger_bilinear, bcs = [bc_Inlet])
schrodinger_A.assemble()
schrodinger_L1 = fem.petsc.create_vector(schrodinger_RHS1)
schrodinger_L2 = fem.petsc.create_vector(schrodinger_RHS2)

schrodinger_solve = PETSc.KSP().create(domain.comm)
schrodinger_solve.setOperators(schrodinger_A)
schrodinger_solve.setType(PETSc.KSP.Type.CG)


# Defining a function to compute velocity field from wave function
def wave_derivative_i(u1, u2, i):
    return hbar * (- imaginaryUnit * ( grad(u1)[i] * conj(u1) + grad(u2)[i] * conj(u2) ) )

uTildex = real(wave_derivative_i(PsiTilde1, PsiTilde2, 0))
uTildey = real(wave_derivative_i(PsiTilde1, PsiTilde2, 1))

sourceTerm = grad(uTildex)[0] + grad(uTildey)[1]

ConversionSpace = FunctionSpace(domain, ("DG", 0))  

ConversionSpace = FunctionSpace(domain, ("DG", 0))  
tmpExpression = fem.Expression(sourceTerm, ConversionSpace.element.interpolation_points)
# other option
tmpExpression = fem.Expression(uTildex, ConversionSpace.element.interpolation_points)
tmpFunction = Function(ConversionSpace)
tmpFunction.interpolate(tmpExpression)
4

Thanks Mr. Dokken, after several trials I managed to tackle this problem. Your advice has been very useful.
Thanks a lot