Help understanding correct equation form in complex mode

Errors
,
1

I am a bit confused on which functions are not compatible (and how to modify them) for dolfinx complex mode. In below example, all 3 versions work for real mode, but only version 3 works for complex mode. I would like to know how to modify version 1 or version 2 to make them work in complex mode. Thank you for any help or suggestions! (This example runs with dolfinx v0.8.0)

import dolfinx as dfx
from dolfinx import fem
import dolfinx.fem.petsc
import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
import ufl
import basix

if MPI.COMM_WORLD.rank == 0:
    print("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")
    if np.issubdtype(dfx.default_scalar_type, np.complexfloating):
        print(f"........................... DOLFINx version {dfx.__version__} is running in COMPLEX mode ...........................")
    else:
        print(f"........................... DOLFINx version {dfx.__version__} is running in REAL mode ...........................")
    print("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++")


Lx, Ly, Lz = 1.0,1.0,1.0
nx, ny, nz = 10, 10, 1

domain = dfx.mesh.create_rectangle(MPI.COMM_WORLD, [[0.0,0.0], [Lx,Ly]],[nx,ny], dfx.mesh.CellType.quadrilateral)
metadata_={'quadrature_degree':4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata_)
dx = ufl.Measure('dx', domain=domain, metadata={'quadrature_degree': metadata_['quadrature_degree'], 'quadrature_rule':'default'})
Vu_el = basix.ufl.element("CG", domain.basix_cell(), 1, shape=(domain.topology.dim,))
Vu = fem.functionspace(domain, Vu_el)
u = fem.Function(Vu)
u_test = ufl.TestFunction(Vu)
du = ufl.TrialFunction(Vu)

E = PETSc.ScalarType(7.1e10)
nu = PETSc.ScalarType(0.35)
mu = fem.Constant(domain, E / (2.0*(1.0 + nu)))
lmbda = fem.Constant(domain, E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))) 

eval_version = 1 # --> {1 , 2, or 3}
#...................................................................................................

if eval_version == 1:
    # >>>>>>>>>>>> Version 1 does not work ---
    epsilon_u = ufl.sym(ufl.grad(u))
    Psi = mu * ufl.inner(epsilon_u, epsilon_u) + lmbda/2.0 * ufl.tr(epsilon_u) * ufl.tr(epsilon_u)
    Res = ufl.derivative(Psi, u, u_test)*dx

if eval_version == 2:
    # >>>>>>>>>>>> Version 2 does not work ---
    I = ufl.variable(ufl.Identity(len(u)))
    F  = ufl.variable(I + ufl.grad(u))
    epsilon_u = ufl.variable(ufl.sym(F - I))
    Psi = mu * ufl.inner(epsilon_u, epsilon_u) + lmbda/2.0 * ufl.tr(epsilon_u) * ufl.tr(epsilon_u)
    PK1 = ufl.diff(Psi, F) 
    Res = ufl.inner(PK1, ufl.grad(u_test))*dx

if eval_version == 3:
    # >>>>>>>>>>>>> Version 3 works
    I = ufl.variable(ufl.Identity(len(u)))
    PK1 = 2.0 * mu * ufl.sym(ufl.grad(u)) + lmbda * ufl.tr(ufl.sym(ufl.grad(u))) * I
    Res = ufl.inner(PK1, ufl.grad(u_test))*dx

#...................................................................................................

M0 = ufl.inner(du, u_test)
M0_form = fem.form(M0*dx)
M_mat = fem.petsc.create_matrix(M0_form)
fem.petsc.assemble_matrix(M_mat, M0_form)
M_mat.assemble()

K0 = ufl.derivative(Res, u, du)
K0_form = fem.form(K0)
K_mat = fem.petsc.create_matrix(K0_form)
fem.petsc.assemble_matrix(K_mat, K0_form)
K_mat.assemble()

if MPI.COMM_WORLD.rank == 0:
    print("M & K mat assembly complete")
2

In complex-valued mode, you have to be incredibly careful when it comes to conjugates.
See for instance: The Poisson problem with complex numbers — FEniCSx tutorial for a brief intro
Version 1 I cannot get to work (would need to dig deep into ufl for that one).
Version 2 can be fixed by using:

    I = ufl.variable(ufl.Identity(len(u)))
    F = ufl.variable(I + ufl.grad(u))
    epsilon_u = ufl.variable(ufl.sym(F - I))
    i, j = ufl.indices(2)
    Psi = mu * epsilon_u[i, j] * epsilon_u[i, j] + lmbda / 2.0 * ufl.tr(
        epsilon_u
    ) * ufl.tr(epsilon_u)
    PK1 = ufl.diff(Psi, F)
    Res = ufl.inner(PK1, ufl.grad(u_test)) * dx
1 Like