Custom Newton Solver NOT converging due to singularity at corners

jt509 · July 19, 2023, 6:04pm

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:

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?
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?
Is there anyway to resolve this? Any hints or ideas in this direction will help.

Thanks!

violetus · July 21, 2023, 12:53am

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.

jt509 · July 24, 2023, 6:46am

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!

Topic		Replies	Views
Newton Solver did not converge mathematics	1	1187	August 13, 2020
Newton solver, nonlinear problem convergence	4	902	November 28, 2020
Why does Solver not converge in specific cases?	1	683	April 27, 2021
The solver did not converge when solving the N-S eqation with large inlet velocity Errors	4	586	November 23, 2021
Newton solver (manually computed) issue in a time loop variational formulation	1	502	March 30, 2022

Custom Newton Solver NOT converging due to singularity at corners

Related Topics