# Custom Newton Solver NOT converging due to singularity at corners

I am able to reproduce the error I am getting in the following minimal example: Consider

\begin{align*} &\sigma \boldsymbol{u} - div(\nu \nabla \boldsymbol{u}) = \boldsymbol{f}(T) \;\; \mbox{in} \;\; \Omega := [0,1] \times [0,1], \\ &-div(D \nabla T) + (\boldsymbol{u} \cdot \nabla) T = 0 \;\; \mbox{in} \;\; \Omega, \\ &\frac{\partial T}{\partial \boldsymbol{n}} = 0 \;\; \mbox{on} \;\; \Gamma_{U}, \Gamma_{D}, \; \boldsymbol{u} = \boldsymbol{0} \;\; \mbox{on} \;\; \Gamma, \; T = T_0 \;\; \mbox{on} \;\; \Gamma_L, \; T = T_1 \;\; \mbox{on} \;\; \Gamma_R. \end{align*}
Here, \boldsymbol{f}(T) := (G_r T)g is a linear function of T and G_r is a positive constant and g is a vector. I tried solving this problem using the built in Newton solver (with parameter choices specified in the code):

from dolfin import *
import matplotlib.pyplot as plt

# Parameters
Ra = Constant(100);  Rk = Constant(1); Pr = Constant(10); Da = Constant(10E-7); Le = Constant(0.8);
Sr = Constant(0); Du = Constant(.1); N = Constant(0);
GrT = Ra/(Pr*Da); GrC = N*GrT
Sc = Le*Pr
g = Constant((1,0))
nu  = Constant(1)
sigma = 1/Da       # Reaction Domination

#D
DT = as_matrix([[Rk/Pr, Du],[Sr, 1/Sc]])

# T0 and T1
TD0 = Constant(-1); TD1 = Constant(1);

# Mollified at the ends
# TD0 = Expression('(x[1]>0 && x[1] < 1) ? -1/(exp(1/((x[1])*(1-x[1])))) : 0', degree = 4)
# TD1 = Expression('(x[1]>0 && x[1] < 1) ?  1/(exp(1/((x[1])*(1-x[1])))) : 0', degree = 4)

uD = Constant((0,0))

mesh = UnitSquareMesh(30,30)

# Finite element spaces
cell = triangle;
k = int(1)
P1v = VectorElement('CR', cell, k)
P1s = FiniteElement('CR', cell, k)
elem = MixedElement([P1v, P1s])
W = FunctionSpace(mesh, elem)
v, s = TestFunctions(W)
delta = Function(W)
u, T  = split(delta)

# Boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

bc_uT = DirichletBC(W.sub(0), uD, boundaries, 2)
bc_uB = DirichletBC(W.sub(0), uD, boundaries, 4)
bc_uL = DirichletBC(W.sub(0), uD, boundaries, 1)
bc_uR = DirichletBC(W.sub(0), uD, boundaries, 3)

bc_TL = DirichletBC(W.sub(1), TD1, boundaries, 1)
bc_TR = DirichletBC(W.sub(1), TD0, boundaries, 3)

bc = [bc_uT, bc_uB, bc_uL, bc_uR,
bc_TL, bc_TR]

he = FacetArea(mesh)
n  = FacetNormal(mesh)

StabC = sqrt(sigma)*10
# Penalty to handle reaction domination
PenaltyState = (StabC)/avg(he) * nu * dot((jump(u)),jump(v)) * dS + (StabC)/he*dot((nu*u),v) * ds
Eq1s = (inner(nu*grad(u), grad(v))*dx) + sigma*dot(u,v)*dx + PenaltyState

FY = (GrT*T)*g
Eq8s = inner(v,FY)*dx

# unwinding
upwindT = (abs(dot(avg(u),n('+')))*(jump(s)*jump(T)))*dS

Eq1Ts = dot((DT*grad(T)), grad(s))*dx +  dot(u,grad(T))*s*dx + upwindT

F = Eq1s + Eq1Ts - Eq8s

# Solve
solve(F == 0, delta, bc, solver_parameters={"newton_solver":{"relative_tolerance": 1e-6}})
uh, Th = delta.split()

plt.figure(1)
plt.title('u_h')
pu = plot(uh)
plt.colorbar(pu)

plt.figure(2)
plt.title('T_h')
pT = plot(Th)
plt.colorbar(pT)

plt.show()


Observe that due to the choice of parameters, the boundary conditions T_0 = -1 and T_1 = 1 create singularities at the left and right corners of the domain, which the built in Newton solver is able to handle. Following is the output:

(fenicsproject) jaitushar@Jais-MacBook-Pro ~ %  cd /Users/jaitushar ; /usr/bin/env /opt/anaconda3/envs/fenicsproject/bin/python /Users/jaitushar/.vscode/extens
ions/ms-python.python-2023.10.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher 64109 -- /Users/jaitushar/Documents/myFEniCS/myExamples/MinimalEx
_II_Newton.py
No Jacobian form specified for nonlinear variational problem.
Differentiating residual form F to obtain Jacobian J = F'.
Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 2.329e+01 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = 2.135e-01 (tol = 1.000e-10) r (rel) = 9.167e-03 (tol = 1.000e-06)
Newton iteration 2: r (abs) = 4.348e+00 (tol = 1.000e-10) r (rel) = 1.867e-01 (tol = 1.000e-06)
Newton iteration 3: r (abs) = 1.019e+00 (tol = 1.000e-10) r (rel) = 4.374e-02 (tol = 1.000e-06)
Newton iteration 4: r (abs) = 1.775e-01 (tol = 1.000e-10) r (rel) = 7.623e-03 (tol = 1.000e-06)
Newton iteration 5: r (abs) = 8.842e-03 (tol = 1.000e-10) r (rel) = 3.797e-04 (tol = 1.000e-06)
Newton iteration 6: r (abs) = 2.188e-05 (tol = 1.000e-10) r (rel) = 9.393e-07 (tol = 1.000e-06)
Newton solver finished in 6 iterations and 6 linear solver iterations.


However, when I code the same problem using a custom newton solver by providing Jacobian and residual manually and using a linear solver (mumps), the solver is not able to resolve the singularity and fluctuates. Following is the code:

# Import Libraries
from dolfin import *
import matplotlib.pyplot as plt

Ra = Constant(100);  Rk = Constant(1); Pr = Constant(10); Da = Constant(10E-7); Le = Constant(0.8);
Sr = Constant(0); Du = Constant(.1); N = Constant(0);
GrT = Ra/(Pr*Da); GrC = N*GrT
Sc = Le*Pr
g = Constant((1,0))
nu  = Constant(1)
sigma = 1/Da

DT = as_matrix([[Rk/Pr, Du],[Sr, 1/Sc]])
TD0 = Constant(-1); TD1 = Constant(1);
# TD0 = Expression('(x[1]>0 && x[1] < 1) ? -1/(exp(1/((x[1])*(1-x[1])))) : 0', degree = 4)
# TD1 = Expression('(x[1]>0 && x[1] < 1) ?  1/(exp(1/((x[1])*(1-x[1])))) : 0', degree = 4)
uD = Constant((0,0))

mesh = UnitSquareMesh(30,30)

# Finite element spaces
cell = triangle;
k = int(1)
P1v = VectorElement('CR', cell, k)
P1s = FiniteElement('CR', cell, k)
elem = MixedElement([P1v, P1s])
W = FunctionSpace(mesh, elem)
v, s = TestFunctions(W)
u, T = TrialFunctions(W)

# T0D = Expression('(x[1]>0.2+eps && x[1] < 0.8-eps) ? -1/(exp(1/((x[1]-0.2)*(0.8-x[1])))) : 0', eps = DOLFIN_EPS, degree = 2)
# T1D = Expression('(x[1]>0.2+eps && x[1] < 0.8-eps) ?  1/(exp(1/((x[1]-0.2)*(0.8-x[1])))) : 0', eps = DOLFIN_EPS, degree = 2)

# Boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)

class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 1.0)

class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0.0)

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()

# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

bc_uT = DirichletBC(W.sub(0), uD, boundaries, 2)
bc_uB = DirichletBC(W.sub(0), uD, boundaries, 4)
bc_uL = DirichletBC(W.sub(0), uD, boundaries, 1)
bc_uR = DirichletBC(W.sub(0), uD, boundaries, 3)

bc_TL = DirichletBC(W.sub(1), TD1, boundaries, 1)
bc_TR = DirichletBC(W.sub(1), TD0, boundaries, 3)

bc = [bc_uT, bc_uB, bc_uL, bc_uR,
bc_TL, bc_TR]

# Initial guesses satisfying the boundary conditions
u0 = Function(W)

T0 = Function(W)
bc_TL.apply(T0.sub(1).vector()); bc_TR.apply(T0.sub(1).vector())

for bcs in bc:
bcs.homogenize()

# Mesh info
he = FacetArea(mesh)
n  = FacetNormal(mesh)

# Assemble

iter = 0; Er = 1; a_tol = 1E-8

while (Er > a_tol) and (iter < 20):

iter += 1

############### Assemble non-linearities ###############
b = u0.sub(0)
StabC = sqrt(sigma)*10
PenaltyState = (StabC)/avg(he) * nu * dot((jump(u)),jump(v)) * dS + (StabC)/he*dot((nu*u),v) * ds
Eq1s = (inner(nu*grad(u), grad(v))*dx) + sigma*dot(u,v)*dx + PenaltyState

FY = (GrT*T0.sub(1))*g
Eq8s = inner(v,FY)*dx

upwindT = (abs(dot(avg(b),n('+')))*(jump(s)*jump(T)))*dS
Eq1Ts = dot((DT*grad(T)), grad(s))*dx + dot(b,grad(T))*s*dx + upwindT

# System
J = Eq1s + Eq1Ts

res = - action(Eq1s, u0)  + Eq8s - action(Eq1Ts, T0)

# Solve
delta = Function(W)
solve(J == res, delta, bc, solver_parameters={'linear_solver':'mumps'})
un, Tn= delta.split()

# Termination criterion
error = (delta)**2*dx
Er = sqrt(abs(assemble(error)))

# Update the solution
u0.vector()[:] += un.vector()[:]
T0.vector()[:] += Tn.vector()[:]

print('|delta| = %e, iteration %d \n' % (Er, iter))

## Visualize
plt.figure(1)
plt.title('u_h')
pu = plot(u0.sub(0))
plt.colorbar(pu)

plt.figure(2)
plt.title('T_h')
pT = plot(T0.sub(1))
plt.colorbar(pT)

plt.show()


Output:

(fenicsproject) jaitushar@Jais-MacBook-Pro ~ %  cd /Users/jaitushar ; /usr/bin/env /opt/anaconda3/envs/fenicsproject/bin/python /Users/jaitushar/.vscode/extens
ions/ms-python.python-2023.10.1/pythonFiles/lib/python/debugpy/adapter/../../debugpy/launcher 64168 -- /Users/jaitushar/Documents/myFEniCS/myExamples/MinimalEx
_II.py
Solving linear variational problem.
|delta| = 5.676462e-01, iteration 1

Solving linear variational problem.
|delta| = 5.546556e+00, iteration 2

Solving linear variational problem.
|delta| = 4.399264e-01, iteration 3

Solving linear variational problem.
|delta| = 4.358093e+00, iteration 4

Solving linear variational problem.
|delta| = 1.333091e-01, iteration 5

Solving linear variational problem.
|delta| = 1.318313e+00, iteration 6

Solving linear variational problem.
|delta| = 4.795220e-02, iteration 7

Solving linear variational problem.
|delta| = 4.542295e-01, iteration 8

Solving linear variational problem.
|delta| = 1.668294e-02, iteration 9

Solving linear variational problem.
|delta| = 1.470277e-01, iteration 10

Solving linear variational problem.
|delta| = 5.031541e-03, iteration 11

Solving linear variational problem.
|delta| = 4.315923e-02, iteration 12

Solving linear variational problem.
|delta| = 1.448626e-03, iteration 13

Solving linear variational problem.
|delta| = 1.234712e-02, iteration 14

Solving linear variational problem.
|delta| = 4.584585e-04, iteration 15

Solving linear variational problem.
|delta| = 3.990136e-03, iteration 16

Solving linear variational problem.
|delta| = 2.415253e-04, iteration 17

Solving linear variational problem.
|delta| = 2.231978e-03, iteration 18

Solving linear variational problem.
|delta| = 2.039454e-04, iteration 19

Solving linear variational problem.
|delta| = 1.926179e-03, iteration 20


Questions:

1. I suspect ill conditioning is being introduced due to the boundary conditions and the problem. How is the Newton solver able to handle this but not the custom newton solve using the linear solver that I wrote?

2. Does the Newton solver introduce some kind of preconditioning in the system by default? If so how can I use it in my linear solve step?

3. Is there anyway to resolve this? Any hints or ideas in this direction will help.

Thanks!

You might want to have a look at Default absolute tolerance and relative tolerance.

I think it is rather an issue of Jacobian, kindly verify the Jacobian if it is defined correctly. Also consider reading Solution nonconvergence.

Kindly see source code here.

As I said, please verify Jacobian once if it is correctly built because when you were using the built in Newton solver, Jacobian is defined as -

But when the custom solver is used then you are explicitly Jacobian which might be causing issue.
I hope this helps.

Hi @violetus. I can confirm that it was indeed a case of malformed Jacobian alongwith providing a more refined meshes (to resolve the singularity). Thanks!