Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09) Newton solver finished in 0 iterations and 0 linear solver iterations

General
1

Hi, I’m trying to implement this tutorial on Hertzian contact using penalty approach by @bleyerj in 2D on an unit square domain but unfortunately ran into the error stated above in title. I tried using Custom NewtonSolver made by @nate in this post Set krylov linear solver paramters in newton solver instead of inbuilt NonlinearVariationalProblem but ran into the same error again.

Here is MWE :-

from dolfin import *
from mshr import *
import numpy as np

square = Rectangle(Point(0, 0), Point(1., 1.))
domain = square 
mesh = generate_mesh(domain, 10)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1., 1e-8) and on_boundary
def symmetry_x(x, on_boundary):
        return near(x[0], 0., 1e-8) and on_boundary
def symmetry_y(x, on_boundary):
        return near(x[1], 1., 1e-8) and on_boundary
def bottom(x, on_boundary):
        return near(x[1], 0., 1e-8) and on_boundary

facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
Top().mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d+(pow(x[0],2))/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom),  
     DirichletBC(V.sub(0), Constant(0.), symmetry_x),
     DirichletBC(V.sub(1), Constant(0.), symmetry_y)]

E = Constant(10.)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u, du)

class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)

class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("ksp_monitor")
        PETScOptions.set("pc_type", "ilu")

        self.linear_solver().set_from_options()

problem = Problem(J, form, bc)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

Any help would be highly appreciated.
Thanks.

2

I don’t have time to check myself but consider the tips here: Default absolute tolerance and relative tolerance - #4 by nate

Also for the system to be solvable you may need to ensure both bodies are initially in contact.

1 Like
3

Thanks @nate for suggesting the thread, the error got resolved when I defined Dirichlet boundary condition using obstacle expression like this -

bc =[DirichletBC(V, Constant((0., 0.)), bottom), 
     DirichletBC(V.sub(1), obstacle, bc_top), 
     DirichletBC(V.sub(0), Constant(0.), symmetry_x), 
     DirichletBC(V.sub(0), Constant(0.), symmetry_xx)]

But the problem now is I’m getting way off values compared to analytical solution that is -

print("Maximum pressure FE: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FE: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))

gives

Maximum pressure FE: 24.09987 Hertz:  1.39916
Applied force    FE:  0.00000 Hertz:  0.02930

where the simulation values must be close to the analytical but they are not, I don’t know what mistake I am making. Kindly suggest how to solve this.
Thanks.

PS: MWE :-

from dolfin import *
from mshr import *
import numpy as np

square = Rectangle(Point(0, 0), Point(1., 1.))
domain = square 
mesh = generate_mesh(domain, 10)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1., 1e-8) and on_boundary
def symmetry_x(x, on_boundary):
        return near(x[0], 0., 1e-8) and on_boundary
def symmetry_xx(x, on_boundary):
        return near(x[0], 1., 1e-8) and on_boundary
def bottom(x, on_boundary):
        return near(x[1], 0., 1e-8) and on_boundary

facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
bc_top = Top()
bc_top.mark(facets, 1)
ds = Measure('ds', subdomain_data=facets)

R = 0.5
d = 0.02
obstacle = Expression("-d + (pow(x[0] - 0.5,2) + 1)/2/R", d=d, R=R, degree=2)

V = VectorFunctionSpace(mesh, "CG", 1)
V2 = FunctionSpace(mesh, "CG", 1)
V0 = FunctionSpace(mesh, "DG", 0)

u = Function(V, name="Displacement")
du = TrialFunction(V)
u_ = TestFunction(V)
gap = Function(V2, name="Gap")
p = Function(V0, name="Contact pressure")

bc =[DirichletBC(V, Constant((0., 0.)), bottom), 
     DirichletBC(V.sub(1), obstacle, bc_top), 
     DirichletBC(V.sub(0), Constant(0.), symmetry_x), 
     DirichletBC(V.sub(0), Constant(0.), symmetry_xx)]

E = Constant(10.)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def eps(v):
    return sym(grad(v))
def sigma(v):
    return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
def ppos(x):
    return (x+abs(x))/2.

pen = Constant(1e4)
form = inner(sigma(u), eps(u_))*dx + pen*dot(u_[1], ppos(u[1]-obstacle))*ds(1)
J = derivative(form, u, du)

class Problem(NonlinearProblem):
    def __init__(self, J, F, bcs):
        self.bilinear_form = J
        self.linear_form = F
        self.bcs = bcs
        NonlinearProblem.__init__(self)

    def F(self, b, x):
        assemble(self.linear_form, tensor=b)
        for bc in self.bcs:
            bc.apply(b, x)

    def J(self, A, x):
        assemble(self.bilinear_form, tensor=A)
        for bc in self.bcs:
            bc.apply(A)

class CustomSolver(NewtonSolver):
    def __init__(self):
        NewtonSolver.__init__(self, mesh.mpi_comm(),
                              PETScKrylovSolver(), PETScFactory.instance())

    def solver_setup(self, A, P, problem, iteration):
        self.linear_solver().set_operator(A)

        PETScOptions.set("ksp_type", "gmres")
        PETScOptions.set("ksp_monitor")
        PETScOptions.set("pc_type", "ilu")

        self.linear_solver().set_from_options()

problem = Problem(J, form, bc)
custom_solver = CustomSolver()
custom_solver.solve(problem, u.vector())

def local_project(v, V, u=None):
    dv = TrialFunction(V)
    v_ = TestFunction(V)
    a_proj = inner(dv, v_)*dx
    b_proj = inner(v, v_)*dx
    solver = LocalSolver(a_proj, b_proj)
    solver.factorize()
    if u is None:
        u = Function(V)
        solver.solve_local_rhs(u)
        return u
    else:
        solver.solve_local_rhs(u)
        return

p.assign(-project(sigma(u)[1, 1], V0))
gap.assign(project(obstacle-u[1], V2))

a = sqrt(R*d)
F = 4/3.*float(E)/(1-float(nu)**2)*a*d
p0 = 3*F/(2*pi*a**2)
print("Maximum pressure FE: {0:8.5f} Hertz: {1:8.5f}".format(max(np.abs(p.vector().get_local())), p0))
print("Applied force    FE: {0:8.5f} Hertz: {1:8.5f}".format(4*assemble(p*ds(1)), F))