# Newton's Solver does not converge

**URL:** <https://fenicsproject.discourse.group/t/newtons-solver-does-not-converge/8105>\
**Category:** Uncategorized\
**Created:** [April 20, 2022, 9:53pm UTC](https://fenicsproject.discourse.group/t/newtons-solver-does-not-converge/8105 "2022-04-20T21:53:52Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![pm0212](https://avatars.discourse-cdn.com/v4/letter/p/74df32/32.png) [@pm0212](https://fenicsproject.discourse.group/u/pm0212)\
**Post date:** [April 20, 2022, 9:53pm UTC](https://fenicsproject.discourse.group/t/newtons-solver-does-not-converge/8105/1 "2022-04-20T21:53:52Z")

</div>

Hello! me and my group mates are trying to run a nonlinear structural problem to solve for the variable ho0 and are encountering the error that Newton’s solver was unable to solve the Non-Linear system and the reason being that maximum number of iterations are reached and the solver did NOT converge.  
Following is the minimal code for our example-

```auto
from dolfin import *
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

tol=1E-12

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 2) # upped dimension

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
plt.show()

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

top = Top()

boundarie = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundarie.set_all(0)
top.mark(boundarie, 1)

ds = ds(subdomain_data=boundarie) 

ho0=TrialFunction(V)
s=TestFunction(V)
t = TestFunction(V)

ho0=np.random.rand(1) # initial guess

ho0= Function(V)

F0 = dot((-1*(ho0**2)/0.001),s)* ds(1) + dot((3-u),t)*ds(1) + Constant(0) * ho0 * s * dx

ho0J=derivative(F0,ho0)

Ho0problem = NonlinearVariationalProblem(F0, ho0, None , ho0J)

Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0prm = Ho0solver.parameters

Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0prm = Ho0solver.parameters
# makeing the error more visible through less iterations and lowering the tolerances to
# increase the chance of convergence
Ho0prm['newton_solver']['absolute_tolerance'] = 1E-5
Ho0prm['newton_solver']['relative_tolerance'] = 1E-5
Ho0prm['newton_solver']['maximum_iterations'] = 5 
Ho0prm['newton_solver']['relaxation_parameter'] = 1.0

Ho0solver.solve()

ho0plot=plot(ho0, title='show ho0')
plt.colorbar(ho0plot)
plt.show()

print ('rho0 calculated')

```

Some additional background-  
Previously we were facing the PETSc function error and we were suggested that the error occurs because of having no integral over the domain in our definition of F0, so we have made changes to that. However, solver did NOT converge. Few more changes we implemented include-

1. Scaling up the dimension of the function space to 2.
2. Increasing the tolerance value
3. Trying out SNES solver

Any help is much appreciated. Thanks!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [April 20, 2022, 10:20pm UTC](https://fenicsproject.discourse.group/t/newtons-solver-does-not-converge/8105/2 "2022-04-20T22:20:57Z")

</div>

Integrating something that is Zero over the rest of the domain is not going to improve the iterative solver.

> [@pm0212](#):
>
> `Constant(0) * ho0 * s * dx`

Do you actually only want to solve a PDE on the boundary of the mesh? If so, have you considered @cdaversin’s meshview?  
Ref

> **[Abstractions and Automated Algorithms for Mixed Domain Finite Element Methods...](https://dl.acm.org/doi/abs/10.1145/3471138)**
>
> Mixed dimensional partial differential equations (PDEs) are equations coupling unknown
> fields defined over domains of differing topological dimension. Such equations naturally
> arise in a wide range of scientific fields including geology, physiology,...

> **[Mixed-dimensional-examples v2019.1](https://zenodo.org/records/3525001)**
>
> Supplementary material (code) for the numerical experiments presented in "Abstractions and automated algorithms for mixed-dimensional finite element methods" by C. Daversin-Catty, C. N. Richardson, Ada J. Ellingsrud and M. E. Rognes (2019)....

---

<div class="post-metadata">

**Author:** ![Mister\_Worst](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mister_worst/32/3444_2.png) [@Mister\_Worst](https://fenicsproject.discourse.group/u/Mister_Worst)\
**Post date:** [April 29, 2022, 12:42am UTC](https://fenicsproject.discourse.group/t/newtons-solver-does-not-converge/8105/3 "2022-04-29T00:42:04Z")

</div>

Hi,  
Thank you for your help,  
I’m in the same group as pm0212 and we have been trying to implement the method you suggested but we keep running into two problems depending on how we modify the code.

Here is a code for our current code

```auto
from dolfin import *
from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as plt

tol=1E-12

mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'Lagrange', 2) # upped dimension

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

Laplace=plot(u, title='show u')
plt.colorbar(Laplace)
#plt.show()

###########################################################
marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
    for k in facets(mesh):
        marker[k] = 0.5 - DOLFIN_EPS < k.midpoint().x() < 0.5 + DOLFIN_EPS

## Create submesh ##
submesh = MeshView.create(marker, 1)

Q = FunctionSpace(mesh, "Lagrange", 1)
LM = FunctionSpace(submesh, "Lagrange", 1)
W = MixedFunctionSpace(Q,LM)

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

top = Top()

boundarie = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundarie.set_all(0)
top.mark(boundarie, 1)

dBound = ds(subdomain_data=boundarie) 
dSurf = Measure("ds", domain=W.sub_space(1).mesh())

ho0,NotUsed0=TrialFunctions(W)
(s,s_lm)=TestFunctions(W)

#F0 = dot((-1*(ho0**2)/0.001),s)* dBound(1) + dot((3-u),s_lm)*dBound(1)
F0 = dot((-1*(ho0**2)/0.001),s)* dSurf + dot((3-u),s_lm)*dSurf

ho0J=derivative(F0,ho0)
print("ho0J = ", ho0J)
  
ho0= Function(W)

Ho0problem = NonlinearVariationalProblem(F0, ho0, None , ho0J)

Ho0solver = NonlinearVariationalSolver(Ho0problem)

Ho0solver.solve()

ho0plot=plot(ho0, title='show ho0')
plt.colorbar(ho0plot)
plt.show()

print ('ho0 calculated')

```

The problems we keep running into are either

```auto
Computing derivative of form w.r.t. 'v_1^0'. Supply Function as a Coefficient
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/formmanipulations.py", line 85, in derivative
    raise RuntimeError("Computing derivative of form w.r.t. '{}'. Supply Function as a Coefficient".format(u))
  File "/home/ucarsten/Desktop/FEniCS/Example-ho0-poisson_mod.py", line 87, in <module>
    ho0J=derivative(F0,ho0)

```

Which seems to be a problem with the derivative code probably because we now got the two FunctionSpaces and the two TrialFunctions.

or when we use the other F0 function

> F0 = dot((-1\*(ho0\*\*2)/0.001),s)\* dBound(1) + dot((3-u),s\_lm)\*dBound(1)  
> We get

```auto
'<' not supported between instances of 'Mesh' and 'Mesh'
  File "/home/ucarsten/Desktop/FEniCS/Example-ho0-poisson_mod.py", line 84, in <module>
    F0 = dot((-1*(ho0**2)/0.001),s)* dBound(1) + dot((3-u),s_lm)*dBound(1)

```

Sadly we have not been able to resolve either problem yet.  
Again thank you for your help any further suggestions would be appreciated.