Nonlinear complex-value problem

variational formulation
1

Hello,
I am trying to solve a nonlinear complex-value problem with item of 1j*|u|^2*u. When I create the weak form in axi-symmetric for this item as:
(1j)*(ufl.inner(u,u))*u*ufl.conj(v)*xm[0]*ufl.dx
which bring error of “ArityMismatch: Adding expressions with non-matching form arguments (‘v_1’,) vs (‘conj(v_1)’,)”

Does any one know how to solve this issue.

Below is the full code:

# %%

import dolfinx

from mpi4py import MPI

import numpy as np

from dolfinx.fem import functionspace, Function,Constant

from dolfinx.mesh import create_unit_cube,create_unit_square,create_rectangle

import matplotlib.pyplot as plt

from ufl import TrialFunction, TestFunction, inner,grad

import ufl

from petsc4py import PETSc

import dolfinx.fem.petsc

from dolfinx.nls.petsc import NewtonSolver

comm = MPI.COMM_WORLD

comm_rank = comm.Get_rank()

numProc=comm.Get_size()

# %%

nEle=20

#msh=create_unit_cube(comm,nEle,nEle,nEle)

msh = create_rectangle(comm,points=((0.0, 0.0), (1.0, 2.0)), n=(16, 32))

xm = ufl.SpatialCoordinate(msh)

V = functionspace(msh, ("Lagrange", 1))

#A = Function(V, dtype=np.complex128)

u_c = Function(V, dtype=np.complex128)

u_c.interpolate(lambda x:(x[1]+(x[0]**2 + x[1])*1j))

#u = TrialFunction(V)

u=Function(V, dtype=np.complex128)

v = TestFunction(V)

f = Function(V, dtype=np.complex128)

f.interpolate(lambda x:(1j+5-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]*1j-x[0]**2-x[1])))

k=Constant(msh,PETSc.ScalarType(((1,0),(0,0))))

# %%

a =(u.dx(1)*ufl.conj(v))*xm[0]*ufl.dx+(1j)*ufl.inner(u.dx(0),v.dx(0)) *xm[0]* ufl.dx+(1j)*(ufl.inner(u,u))*u*ufl.conj(v)*xm[0]*ufl.dx

L = ufl.inner(f, v)*xm[0] * ufl.dx

F=a-L

# %%

msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)

boundary_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)

boundary_dofs = dolfinx.fem.locate_dofs_topological(V, msh.topology.dim-1, boundary_facets)

bc = dolfinx.fem.dirichletbc(u_c, boundary_dofs)

problem = dolfinx.fem.petsc.NonlinearProblem(F, u, bcs=[bc])

solver = NewtonSolver(comm, problem)

# Set Newton solver options

solver.atol = 1e-8

solver.rtol = 1e-8

solver.convergence_criterion = "incremental"

# We can customize the linear solver used inside the NewtonSolver by

# modifying the PETSc options

ksp = solver.krylov_solver

opts = PETSc.Options() # type: ignore

option_prefix = ksp.getOptionsPrefix()

opts[f"{option_prefix}ksp_type"] = "gmres"

opts[f"{option_prefix}pc_type"] = "lu"

ksp.setFromOptions()

r = solver.solve(u)

print("Number of iterations: %d, converged %d"%(r[0],r[1]))

# %%

u.name='numerical'

u_c.name='exact'

file_results= dolfinx.io.XDMFFile(MPI.COMM_WORLD, "./fenicsx_complex/example.xdmf","w" )

file_results.write_mesh(msh)

file_results.write_function(u,0)

file_results.write_function(u_c,0)

file_results.close()

# %%

max_error = (np.max(np.abs(u_c.x.array-u.x.array)))

idx= (np.argmax(np.abs(u_c.x.array-u.x.array)))

print(max_error)

print(u_c.x.array[idx],u.x.array[idx])

fig=plt.figure(1)

ax=plt.subplot(1,1,1)

ax.plot(u_c.x.array.imag[:], 'o', label='exact')

ax.plot(u.x.array.imag[:], '+', label='Fenicsx')

ax.legend()

fig=plt.figure(2)

ax=plt.subplot(1,1,1)

ax.plot(u_c.x.array.real[:], 'o', label='exact')

ax.plot(u.x.array.real[:], '+', label='Fenicsx')

ax.legend()
2

Please edit your post to format properly the code using “```”, otherwise it is difficult to read it.

Can you write your form without explicitly using ufl.conj? If you consistently use ufl.inner everwhere, it will take care on its own to put a ufl.conj where appropriate.

3

Thanks for reply @francesco-ballarin.
I edit the post with code block.

I also remove the ufl.conj and use ufl.inner as

a =ufl.inner(u.dx(1),v)*xm[0]*ufl.dx+(1j)*ufl.inner(u.dx(0),v.dx(0)) *xm[0]* ufl.dx+(1j)*(ufl.inner(u,u))*ufl.inner(u,v)*xm[0]*ufl.dx # ufl.imag(u)**2+ufl.real(u)**2
L = ufl.inner(f, v)*xm[0] * ufl.dx
F=a-L

But I still got the error of
ArityMismatch: Adding expressions with non-matching form arguments ('v_1',) vs ('conj(v_1)',).
if I remove the (1j)*(ufl.inner(u,u))*ufl.inner(u,v)*xm[0]*ufl.dx, the error gone.

4

Should this be inner or dot(u, u). The latter suggestion compiles.

This is because one has to compute the Jacobian of F to solve the non-linear problem, and using ufl.derivative(...) on the term above gives you contributions in both parts of the inner product.

5

Thanks @dokken.
yes, dot(u,u) compiles. But I would like to calculate |u|^2=u*\bar{u}. Does dot(u,u) evaluate u*u or u*\bar{u}? I guess the first one.
BTW, I also tried (1j)*(ufl.imag(u)**2+ufl.real(u)**2)*ufl.inner(u,v)*xm[0]*ufl.dx
but got error of
ArityMismatch: Applying nonlinear operator Imag to expression depending on form argument v_1.

6

You would have to carefully write out the derivative of your non-linear form by hand. Using automatic differentiation on non-linear complex valued problems is a bit problematic, as the forms are sesquilinear, and there has be a consistent usage of trial and test functions.

I think you could reproduce your error by calling ufl.derivative(F, u).

7

Thanks @dokken.

I also tried to solve it using mixed-element. I am solving the manufacturing problem of
\frac{\partial u}{\partial z}=i\left(\frac{1}{r} \frac{\partial u}{\partial r}+\frac{\partial^2 u}{\partial r^2}\right)+i|u|^2u+f
Assume the exact solution is u=z+i(r^2 + z)
so f=5+(z^2+(r^2+z)^2)(r^2+z)+i(1-z(z^2+(r^2+z)^2))

using mixed-element:

\frac{\partial u_{real}}{\partial z}=-\left(\frac{1}{r} \frac{\partial u_{imag}}{\partial r}+\frac{\partial u^2_{imag}}{\partial r^2}\right)-(u_{real}^2+u_{imag}^2)u_{imag}+ f_{real}
\frac{\partial u_{imag}}{\partial z}=\left(\frac{1}{r} \frac{\partial u_{real}}{\partial r}+\frac{\partial u^2_{real}}{\partial r^2}\right)+(u_{real}^2+u_{imag}^2)u_{real}+f_{imag}

But it can not converge. Did I miss something?
Below is the full code:

# %%
import dolfinx
from mpi4py import MPI
import numpy as np
from dolfinx.fem import functionspace, Function,Constant
from dolfinx.mesh import create_unit_cube,create_unit_square,create_rectangle
import matplotlib.pyplot as plt
from ufl import TrialFunction, TestFunction, inner,grad
import ufl
from petsc4py import PETSc
import dolfinx.fem.petsc
from basix.ufl import element, mixed_element
from dolfinx.nls.petsc import NewtonSolver
comm = MPI.COMM_WORLD
comm_rank = comm.Get_rank()
numProc=comm.Get_size()

# %%
nEle=20
#msh=create_unit_cube(comm,nEle,nEle,nEle)
msh = create_rectangle(comm,points=((0.0, 0.0), (1.0, 2.0)), n=(16, 32))
xm = ufl.SpatialCoordinate(msh)
P1 = element("Lagrange", msh.basix_cell(), 2)
ME = functionspace(msh, mixed_element([P1, P1]))

q, v = ufl.TestFunctions(ME)
u=Function(ME) # solution
uc=Function(ME) #exact solution
u_re,u_im=ufl.split(u)
uc_re,uc_im=ufl.split(uc)
uc.sub(0).interpolate(lambda x: (x[1])) # real
uc.sub(1).interpolate(lambda x: (x[0]**2 + x[1])) # imag
uc.x.scatter_forward()
V_re, map_re = ME.sub(0).collapse()
V_im, map_im = ME.sub(1).collapse()
# %%

f = Function(ME)
f_re,f_im =ufl.split(u)
f.sub(0).interpolate(lambda x: (5+(x[1]**2+(x[0]**2+x[1])**2)*(x[0]**2+x[1]))) # real
f.sub(1).interpolate(lambda x:(1-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]))) # imag
f.x.scatter_forward()
# k=Constant(msh,PETSc.ScalarType(((1,0),(0,0))))\
# %%
F_re=(u_re.dx(1)*q)*xm[0]*ufl.dx-(u_im.dx(0)*q.dx(0))*xm[0]* ufl.dx-(f_re*q)*xm[0]* ufl.dx \
    + (u_re**2+u_im**2)*u_im*q*xm[0]*ufl.dx
F_im=(u_im.dx(1)*v)*xm[0]*ufl.dx+(u_re.dx(0)*v.dx(0))*xm[0]* ufl.dx-(f_im*v)*xm[0]* ufl.dx \
    - (u_re**2+u_im**2)*u_re*v*xm[0]*ufl.dx
F=F_re+F_im
Jac = ufl.derivative(F, u) # jacobian
# %%
sdim=msh.topology.dim
msh.topology.create_connectivity(sdim-1, sdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(ME, msh.topology.dim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(uc, boundary_dofs)

problem = dolfinx.fem.petsc.NonlinearProblem(F, u,[bc])
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
solver = NewtonSolver(comm, problem)
# Set Newton solver options
#solver.convergence_criterion = "incremental"
solver.error_on_nonconvergence = False
solver.convergence_criterion = "residual" #"residual" #"incremental"
solver.atol = 1e-6
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
r = solver.solve(u)
print("Number of iterations: %d, converged %d"%(r[0],r[1]))


# %%


fig=plt.figure(1)
ax=plt.subplot(1,1,1)
ax.plot(uc.x.array[map_im], 'o', label='exact')
ax.plot(u.x.array[map_im], '+', label='Fenicsx')
ax.legend()


fig=plt.figure(2)
ax=plt.subplot(1,1,1)
ax.plot(uc.x.array[map_re], 'o', label='exact')
ax.plot(u.x.array[map_re], '+', label='Fenicsx')
ax.legend()
8

As a test, set the initial guess for the nonlinear solver to be the interpolation of the exact solution. If it converges in a handful of iterations, then it may mean that the initial guess you currently provide (zero, judging from your code) is not good enough. If it does not converge, or takes forever to do so, you probably are not solving the problem you think, i.e. that solution is not really a solution of the problem as you wrote it in the code.

9

Thanks @francesco-ballarin.
I set the u as exact solution before solving as u.x.array[:]=uc.x.array

it does not work. But after changed

f = Function(ME)
f_re,f_im =ufl.split(u)
f.sub(0).interpolate(lambda x: (5+(x[1]**2+(x[0]**2+x[1])**2)*(x[0]**2+x[1]))) 
f.sub(1).interpolate(lambda x:(1-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]))) 
f.x.scatter_forward()

to

x = ufl.SpatialCoordinate(msh)
f_re=(5+(x[1]**2+(x[0]**2+x[1])**2)*(x[0]**2+x[1]))
f_im=(1-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]))

It works. Do you know why? Here f_re and f_im are used in

F_re=(u_re.dx(1)*q)*x[0]*ufl.dx-(u_im.dx(0)*q.dx(0))*x[0]* ufl.dx-(f_re*q)*x[0]* ufl.dx \
    + (u_re**2+u_im**2)*u_im*q*x[0]*ufl.dx
F_im=(u_im.dx(1)*v)*x[0]*ufl.dx+(u_re.dx(0)*v.dx(0))*x[0]* ufl.dx-(f_im*v)*x[0]* ufl.dx \
    - (u_re**2+u_im**2)*u_re*v*x[0]*ufl.dx
F=F_re+F_im

Below is the full code

# %%
import dolfinx
from mpi4py import MPI
import numpy as np
from dolfinx.fem import functionspace, Function,Constant
from dolfinx.mesh import create_unit_cube,create_unit_square,create_rectangle
import matplotlib.pyplot as plt
from ufl import TrialFunction, TestFunction, inner,grad
import ufl
from petsc4py import PETSc
import dolfinx.fem.petsc
from basix.ufl import element, mixed_element
from dolfinx.nls.petsc import NewtonSolver
comm = MPI.COMM_WORLD
comm_rank = comm.Get_rank()
numProc=comm.Get_size()

# %%
nEle=20
#msh=create_unit_cube(comm,nEle,nEle,nEle)
msh = create_rectangle(comm,points=((0.0, 0.0), (1.0, 2.0)), n=(16, 32))
x = ufl.SpatialCoordinate(msh)
P1 = element("Lagrange", msh.basix_cell(), 1)
ME = functionspace(msh, mixed_element([P1, P1]))

q, v = ufl.TestFunctions(ME)
u=Function(ME) # solution
uc=Function(ME) #exact solution
u_re,u_im=ufl.split(u)
uc_re,uc_im=ufl.split(uc)
uc.sub(0).interpolate(lambda x: (x[1])) # real
uc.sub(1).interpolate(lambda x: (x[0]**2 + x[1])) # imag
uc.x.scatter_forward()
V_re, map_re = ME.sub(0).collapse()
V_im, map_im = ME.sub(1).collapse()
u.x.array[:]=uc.x.array
# %%

# f = Function(ME)
# f_re,f_im =ufl.split(u)
# f.sub(0).interpolate(lambda x: (5+(x[1]**2+(x[0]**2+x[1])**2)*(x[0]**2+x[1]))) # real
# f.sub(1).interpolate(lambda x:(1-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]))) # imag
# f.x.scatter_forward()
f_re=(5+(x[1]**2+(x[0]**2+x[1])**2)*(x[0]**2+x[1]))
f_im=(1-(x[1]**2+(x[0]**2+x[1])**2)*(x[1]))
# %%
F_re=(u_re.dx(1)*q)*x[0]*ufl.dx-(u_im.dx(0)*q.dx(0))*x[0]* ufl.dx-(f_re*q)*x[0]* ufl.dx \
    + (u_re**2+u_im**2)*u_im*q*x[0]*ufl.dx
F_im=(u_im.dx(1)*v)*x[0]*ufl.dx+(u_re.dx(0)*v.dx(0))*x[0]* ufl.dx-(f_im*v)*x[0]* ufl.dx \
    - (u_re**2+u_im**2)*u_re*v*x[0]*ufl.dx
F=F_re+F_im
Jac = ufl.derivative(F, u) # jacobian
# %%
sdim=msh.topology.dim
msh.topology.create_connectivity(sdim-1, sdim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(msh.topology)
boundary_dofs = dolfinx.fem.locate_dofs_topological(ME, msh.topology.dim-1, boundary_facets)
bc = dolfinx.fem.dirichletbc(uc, boundary_dofs)

problem = dolfinx.fem.petsc.NonlinearProblem(F, u,[bc])
dolfinx.log.set_log_level(dolfinx.log.LogLevel.INFO)
solver = NewtonSolver(comm, problem)
# Set Newton solver options
#solver.convergence_criterion = "incremental"
solver.error_on_nonconvergence = False
solver.convergence_criterion = "residual" #"residual" #"incremental"
solver.atol = 1e-6
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
r = solver.solve(u)
print("Number of iterations: %d, converged %d"%(r[0],r[1]))


# %%


fig=plt.figure(1)
ax=plt.subplot(1,1,1)
ax.plot(uc.x.array[map_im], 'o', label='exact')
ax.plot(u.x.array[map_im], '+', label='Fenicsx')
ax.legend()


fig=plt.figure(2)
ax=plt.subplot(1,1,1)
ax.plot(uc.x.array[map_re], 'o', label='exact')
ax.plot(u.x.array[map_re], '+', label='Fenicsx')
ax.legend()