Divegence problem in nonlinear solver

Errors
1

Hello everybody,
I am trying to solve a strongly coupled, nonlinear problem (6 PDEs) with monolithic time discretization. I encountered a convergence problem as below:

*** Error: Unable to solve nonlinear system with NewtonSolver.
*** Reason: Newton solver did not converge because maximum number of iterations reached.
*** Where: This error was encountered inside NewtonSolver.cpp.
*** Process: 0
so I tried to reduce my problem but the same error existed. So I tried it on a linear problem that I have already solved. However, the nonlinear iterative solver still doesn’t converge.
According to https://fenicsproject.discourse.group/t/default-absolute-tolerance-and-relative-tolerance/3829/3 the newton solver should converge in one iteration for linear problems.
I am sharing my code here. It is a bit long but it is the smallest one that I can share here.
I appreciate any help and hint to solve the problem. My original problem is nonlinear and way complicated than this one

2

This is the linear version which works well

import fenics as fe
import matplotlib.pyplot as plt
import logging
import numpy as np

# geometry
b1 = 1 # width of the domain (m)
h1 = 7 # height of the domain (m)
N_POINTS_X_AXIS = 5
N_POINTS_Y_AXIS = 29

# Initial condition and time increment
t_0 = 0.0
delta_t = 1e-3
t_final = 1
n_t = int(t_final/delta_t) #1000

# parameters

# Density
rho = fe.Constant(0.0)

# Lame parameters
LAME_MU =  5.0e6
LAME_LAMBDA = 8.0e6

DENSITY = 0.0 # b
ACCELERATION_DUE_TO_GRAVITY = 0.016

K_s = 1.0e-5 # intrinsic permeability (m2)
MU_FR = 1.0 # effective dynamic fluid viscosity [Ns/m2]

mesh = fe.RectangleMesh(fe.Point(0,0), fe.Point(b1,h1), N_POINTS_X_AXIS, N_POINTS_Y_AXIS)

u_vector_element = fe.VectorElement("CG", mesh.ufl_cell(), 2)
p_scalar_element = fe.FiniteElement("CG", mesh.ufl_cell(), 1)
mixed_function_space = fe.FunctionSpace(mesh, u_vector_element * p_scalar_element)
'''u_p_init = Expression(('0', '0', 'p_0'), p_0=0, degree=1)  # degree=1

u_p_trial = interpolate(u_p_init, mixed_function_space)
u_p_old.assign(u_solution)'''

# Define trial and test functions

u_p_test = fe.TestFunction(mixed_function_space)
(u_test, p_test) = fe.split(u_p_test)
u_p_trial = fe.TrialFunction(mixed_function_space)
(u_trial, p_trial) = fe.split(u_p_trial)

u_p   = fe.Function(mixed_function_space)
u0_p0  = fe.Function(mixed_function_space)

# initialization

u_p_init = fe.Expression(('0', '0', 'p_0'), p_0=0, degree=1)  # degree=1
u_p = fe.interpolate(u_p_init, mixed_function_space)
u0_p0.assign(u_p)
u, p = fe.split(u_p)
u0, p0 = fe.split(u0_p0)

# Boundary functions

tol = 1e-14

# Top boundary
def Top_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[1], h1, tol))

# Bottom boundary
def Bottom_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[1], 0.0, tol))

# Left boundary
def Left_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[0], 0.0, tol))

# Right boundary
def Right_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[0], b1, tol))

boundary_subdomains = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = fe.AutoSubDomain(Top_boundary)
force_boundary.mark(boundary_subdomains, 3)

# Define measure for boundary condition integral
dss = fe.ds(subdomain_data=boundary_subdomains)

bc_u_bottom = fe.DirichletBC(mixed_function_space.sub(0).sub(1), fe.Constant(0.0), Bottom_boundary)
bc_u_left = fe.DirichletBC(mixed_function_space.sub(0).sub(0), fe.Constant(0.0), Left_boundary)
bc_u_right = fe.DirichletBC(mixed_function_space.sub(0).sub(0), fe.Constant(0.0), Right_boundary)

bc_p_top = fe.DirichletBC(mixed_function_space.sub(1), fe.Constant(0.0), Top_boundary)
bc_vector = [bc_u_bottom, bc_u_left, bc_u_right, bc_p_top]

# Neumann boundary condition (loading)
g = fe.Expression(('0.0*x[0]' ,'-5e6'), element = mixed_function_space.sub(0).ufl_element())

# define strain and stress

def epsilon(u):
        engineering_strain = 0.5 * (fe.nabla_grad(u) + fe.nabla_grad(u).T)
        return engineering_strain

def sigma(u, p):
        cauchy_stress = (
            LAME_LAMBDA * fe.tr(epsilon(u)) * fe.Identity(2)
            +
            2 * LAME_MU * epsilon(u)
            -
            p * fe.Identity(2)
        )
        return cauchy_stress

# Weak forms

from ufl.operators import nabla_grad , nabla_div

weak_form_momentum_balance = (fe.inner(sigma(u_trial, p_trial), nabla_grad(u_test)) * fe.dx -
                      fe.dot(g, u_test) * dss(3))

weak_form_volume_balance = (nabla_div((u_trial - u0)/delta_t) * p_test * fe.dx +
                            K_s / MU_FR * fe.dot(nabla_grad(p_trial), nabla_grad(p_test)) * fe.dx)

weak_form = weak_form_momentum_balance + weak_form_volume_balance

# We have a linear PDE that is separable into a lhs and rhs
weak_form_lhs = fe.lhs(weak_form)
weak_form_rhs = fe.rhs(weak_form)

# Save results in a XDMF file
results = fe.XDMFFile("My_results.xdmf")
results.parameters["flush_output"] = True
results.parameters["functions_share_mesh"] = True

# Solve for each point in time
time_current = t_0
for i in range(int(n_t)):
    time_current += delta_t

    # Finite Element Assembly, BC imprint & solving the linear system
    fe.solve(
        weak_form_lhs == weak_form_rhs,
        u_p,
        bc_vector)
    'ME'
    results.write(u_p.split(deepcopy=True)[0], time_current)
    results.write(u_p.split(deepcopy=True)[1], time_current)
    u0_p0.assign(u_p)
3

Here is the same problem which I tried to solve by nonlinear solver

import fenics as fe
import matplotlib.pyplot as plt
import logging
import numpy as np

# geometry
b1 = 1 # width of the domain (m)
h1 = 7 # height of the domain (m)
N_POINTS_X_AXIS = 5
N_POINTS_Y_AXIS = 29

# Initial condition and time increment
t_0 = 0.0
delta_t = 1e-3
t_final = 1
n_t = int(t_final/delta_t) #1000

# parameters

# Density
rho = fe.Constant(0.0)

# Lame parameters
LAME_MU =  5.0e6
LAME_LAMBDA = 8.0e6

DENSITY = 0.0 # b
ACCELERATION_DUE_TO_GRAVITY = 0.016

K_s = 1.0e-5 # intrinsic permeability (m2)
MU_FR = 1.0 # effective dynamic fluid viscosity [Ns/m2] 

mesh = fe.RectangleMesh(fe.Point(0,0), fe.Point(b1,h1), N_POINTS_X_AXIS, N_POINTS_Y_AXIS)

u_vector_element = fe.VectorElement("CG", mesh.ufl_cell(), 2)
p_scalar_element = fe.FiniteElement("CG", mesh.ufl_cell(), 1)
mixed_function_space = fe.FunctionSpace(mesh, u_vector_element * p_scalar_element)

'''u_p_init = Expression(('0', '0', 'p_0'), p_0=0, degree=1)  # degree=1

u_p_trial = interpolate(u_p_init, mixed_function_space)
u_p_old.assign(u_solution)'''

# Define trial and test functions

u_p_test = fe.TestFunction(mixed_function_space)
(u_test, p_test) = fe.split(u_p_test)
u_p_trial = fe.TrialFunction(mixed_function_space)
(u_trial, p_trial) = fe.split(u_p_trial)

u_p   = fe.Function(mixed_function_space)
u0_p0  = fe.Function(mixed_function_space)

# initialization
u_p_init = fe.Expression(('0', '0', 'p_0'), p_0=0, degree=1)  # degree=1
u_p = fe.interpolate(u_p_init, mixed_function_space)
u0_p0.assign(u_p)

u, p = fe.split(u_p)
u0, p0 = fe.split(u0_p0)

tol = 1e-14

# Top boundary
def Top_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[1], h1, tol))

# Bottom boundary
def Bottom_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[1], 0.0, tol))

# Left boundary
def Left_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[0], 0.0, tol))

# Right boundary
def Right_boundary(x, on_boundary):
    return (on_boundary and fe.near(x[0], b1, tol))

boundary_subdomains = fe.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundary_subdomains.set_all(0)
force_boundary = fe.AutoSubDomain(Top_boundary)
force_boundary.mark(boundary_subdomains, 3)

# Define measure for boundary condition integral
dss = fe.ds(subdomain_data=boundary_subdomains)

bc_u_bottom = fe.DirichletBC(mixed_function_space.sub(0).sub(1), fe.Constant(0.0), Bottom_boundary)
bc_u_left = fe.DirichletBC(mixed_function_space.sub(0).sub(0), fe.Constant(0.0), Left_boundary)
bc_u_right = fe.DirichletBC(mixed_function_space.sub(0).sub(0), fe.Constant(0.0), Right_boundary)

bc_p_top = fe.DirichletBC(mixed_function_space.sub(1), fe.Constant(0.0), Top_boundary)

bc_vector = [bc_u_bottom, bc_u_left, bc_u_right, bc_p_top]

g = fe.Expression(('0.0*x[0]' ,'-5e6'), element = mixed_function_space.sub(0).ufl_element())

# define strain and stress

def epsilon(u):
        engineering_strain = 0.5 * (fe.nabla_grad(u) + fe.nabla_grad(u).T)
        return engineering_strain

def sigma(u, p):
        cauchy_stress = (
            LAME_LAMBDA * fe.tr(epsilon(u)) * fe.Identity(2)
            +
            2 * LAME_MU * epsilon(u)
            -
            p * fe.Identity(2)
        )
        return cauchy_stress

# Weak forms

from ufl.operators import nabla_grad , nabla_div

weak_form_momentum_balance = (fe.inner(sigma(u, p), nabla_grad(u_test)) * fe.dx -
                      fe.dot(g, u_test) * dss(3))

weak_form_volume_balance = (nabla_div((u - u0)/delta_t) * p_test * fe.dx +
                            K_s / MU_FR * fe.dot(nabla_grad(p), nabla_grad(p_test)) * fe.dx)

weak_form = weak_form_momentum_balance + weak_form_volume_balance

results = fe.XDMFFile("My_results.xdmf")
results.parameters["flush_output"] = True
results.parameters["functions_share_mesh"] = True

J  = fe.derivative(weak_form, u_p, u_p_trial)   # Gateaux derivative in dir. of u

time_current = t_0
for i in range(int(n_t)):
    time_current += delta_t

    # Finite Element Assembly, BC imprint & solving the linear system
    fe.solve(
        weak_form==0,
        u_p,
        bc_vector)
    'ME'
    results.write(u_p.split(deepcopy=True)[0], time_current)
    results.write(u_p.split(deepcopy=True)[1], time_current)
    u0_p0.assign(u_p)

And the error appears as:

RuntimeError                              Traceback (most recent call last)
<ipython-input-2-e5c71400c768> in <module>
    138         weak_form==0,
    139         u_p,
--> 140         bc_vector)
    141     'ME'
    142     results.write(u_p.split(deepcopy=True)[0], time_current)

1 frames
/usr/local/lib/python3.7/dist-packages/dolfin/fem/solving.py in _solve_varproblem(*args, **kwargs)
    312             solver = NonlinearVariationalSolver(problem)
    313             solver.parameters.update(solver_parameters)
--> 314             solver.solve()
    315 
    316 

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to solve nonlinear system with NewtonSolver.
*** Reason:  Newton solver did not converge because maximum number of iterations reached.
*** Where:   This error was encountered inside NewtonSolver.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  74a6ac2095059ea6582fe4ea11686c1c723e14f5
*** -------------------------------------------------------------------------