Newton's Solver does not converge

1

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-

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!

2

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

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

1 Like
3

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

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

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

'<' 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.