Big error when solving non linear complex PDE

variational formulation
1

Hello,

I am trying to solve the pde -laplacianU + q*U^2 =0 with with Dirichlet complex boundary values, when solving it for constant values for the function q the error is small but with a complex variable function q the error is big even with higher resolution or degree for Lagrange poly, like in this example :


#(FEniCS 2019.1)

from dolfin import *

import numpy as np



# -----------------------------------------

# Mixed real–imag Newton solve

# -----------------------------------------

def newton_solve_complex(mesh, V, q_fun, uD_fun, p=2, quad_deg=10):
    # Mixed space W = V x V  (real, imag)

    W = FunctionSpace(mesh, MixedElement(\[V.ufl_element(), V.ufl_element()\]))
        
        
        
     # Unknown, tests, trial for Jacobian
        
    w  = Function(W)                 # (u_R, u_I)
        
    u_R, u_I = split(w)

    v_R, v_I = TestFunctions(W)
    
    dw = TrialFunction(W)
    
    
    
    # q(x): split to real/imag via UserExpression -> interpolate to V
    
    class QRe(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.real(q_fun(x)))
    
    def value_shape(self): return ()
    
    class QIm(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.imag(q_fun(x)))
    
    def value_shape(self): return ()
    
    q_R = interpolate(QRe(degree=4), V)
    
    q_I = interpolate(QIm(degree=4), V)
    
    
    
    # Dirichlet data u_D = uD_fun(x) split into real/imag
    
    class UDR(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.real(uD_fun(x)))
    
    def value_shape(self): return ()
    
    class UDI(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.imag(uD_fun(x)))
    
    def value_shape(self): return ()
    
    uD_R = interpolate(UDR(degree=4), V)
    
    uD_I = interpolate(UDI(degree=4), V)
    
    
    
    bc_R = DirichletBC(W.sub(0), uD_R, "on_boundary")
    
    bc_I = DirichletBC(W.sub(1), uD_I, "on_boundary")
    
    bcs  = \[bc_R, bc_I\]
    
    
    
    # Nonlinearity: (u_R + i u_I)^2
    
    u2_R = u_R\*u_R - u_I\*u_I
    
    u2_I = 2.0\*u_R\*u_I
    
    
    
    # q u^2
    
    Re_qu2 = q_R\*u2_R - q_I\*u2_I
    
    Im_qu2 = q_R\*u2_I + q_I\*u2_R
    
    
    
    dxm = dx(metadata={"quadrature_degree": quad_deg})
    
    
    
    # Residual (two real equations)
    
    F_R = inner(grad(u_R), grad(v_R))\*dxm + Re_qu2\*v_R\*dxm
    
    F_I = inner(grad(u_I), grad(v_I))\*dxm + Im_qu2\*v_I\*dxm
    
    F   = F_R + F_I
    
    
    
    \# Jacobian
    
    J = derivative(F, w, dw)
    
    
    
    prm = {"newton_solver": {"absolute_tolerance": 1e-10,
    
    "relative_tolerance": 1e-10,
    
    "maximum_iterations": 50,
    
    "linear_solver": "mumps"}}
    
    
    
        solve(F == 0, w, bcs, J=J, solver_parameters=prm,
    
    form_compiler_parameters={"quadrature_degree": quad_deg})
    
    
    
    uR_fun, uI_fun = w.split(deepcopy=True)
    
    return uR_fun, uI_fun



# -----------------------------------------

# Manufactured exact solution & q

# -----------------------------------------

k = 5.0

l = 5.0



def u_exact(x):  # complex

    return np.exp(1j\*(k\*x\[0\] + l\*x\[1\]))



def q_function(x):  # complex

    return -(k\*\*2 + l\*\*2) \* np.exp(-1j\*(k\*x\[0\] + l\*x\[1\]))



# -----------------------------------------

# Run small test

# -----------------------------------------

if __name__ == "__main__":

    # Mesh and space
    
    resolution = 64
    
    mesh = RectangleMesh(Point(-1, -1), Point(1, 1), resolution, resolution)
    
    V = FunctionSpace(mesh, "Lagrange", 2)
    
    
    
    # Solve
    
    u_R, u_I = newton_solve_complex(mesh, V, q_function, u_exact, p=2, quad_deg=12)
    
    
    
    # Exact functions on V (for errors)
    
    class UExactRe(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.real(u_exact(x)))
    
    def value_shape(self): return ()
    
    class UExactIm(UserExpression):
    
    def eval(self, values, x): values\[0\] = float(np.imag(u_exact(x)))
    
    def value_shape(self): return ()
    
    uR_ex = interpolate(UExactRe(degree=5), V)
    
    uI_ex = interpolate(UExactIm(degree=5), V)
    
    
    
    dxm = dx(metadata={"quadrature_degree": 12})
    
    dsm = ds(metadata={"quadrature_degree": 12})
    
    
    
    # Domain L2 error
    
    err_dom = np.sqrt(assemble(((u_R - uR_ex)\*\*2 + (u_I - uI_ex)\*\*2) \* dxm))
    
    print("Domain L2 error (complex):", err_dom)
    
    
    
    # Boundary normal derivative error
    
    n = FacetNormal(mesh)
    
    dn_num_re = dot(grad(u_R), n);  dn_ex_re = dot(grad(uR_ex), n)
    
    dn_num_im = dot(grad(u_I), n);  dn_ex_im = dot(grad(uI_ex), n)
    
    err_bnd = np.sqrt(assemble(((dn_num_re - dn_ex_re)\*\*2 + (dn_num_im - dn_ex_im)\*\*2) \* dsm))
    
    print("Boundary L2 error of normal derivative (complex):", err_bnd)