Newton Solver Does not converge in mixed formulation

Hello everyone,

I am trying to solve a coupled problem and running into the problem of ‘Newton Solver did not converge because of maximum number of iterations reached’. I have tried to increase the number of iterations and checked the code multiple times not but can not find the problem. I can not solve even one time step. Is there a problem the way I am dealing with mixed spaces with 4 unknowns? Could anyone please help with it? I will try to provide as many details as possible. I am using dolfinx version 0.7.2.

Problem Statement:
Following is the geometry which is fixed on the left side (This edge number is 5 for bc). A load of -4 N is applied at the bottom on the right (on 0.05m length. This edge number is 2 for bc).

Weak Formulation:

The following weak formulation is chosen. To keep it simple and to keep the code shorter, I am avoiding gravity (b = 0) and growth (rho_hat_S and rho_hat_N = 0).

Stress:

Code:

It is a mixed formulation with four unknowns: u - displacement, p - pressure, nS - solid volume fraction and nN - nutrient volume fraction. Here is the code.

from dolfinx import log, default_scalar_type, fem
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl
import dolfinx
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers
from pathlib import Path
import meshio

#------------------------------Reading Geometry-------------------------------#

proc = MPI.COMM_WORLD.rank

def create_mesh(mesh, cell_type, prune_z=False):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:, :2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
    return out_mesh

if proc == 0:
    # Read in mesh
    msh = meshio.read(Path("geometry")/"split_ver.msh")

    # Create and save one file for the mesh, and one file for the facets
    triangle_mesh = create_mesh(msh, "quad", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    domain = xdmf.read_mesh(name="Grid") # cell tags
    ct = xdmf.read_meshtags(domain, name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    facet_tag = xdmf.read_meshtags(domain, name="Grid") # facet tags

fdim = domain.topology.dim - 1
dim = domain.topology.dim    
#-----------------------Function Spaces and Mixed Space-------------------#

u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) # u quad element
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # p linear element
nS_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # nS linear element
nN_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # nN linear element
element = ufl.MixedElement([u_el, p_el, nS_el, nN_el]) # Mixed Element

V = fem.FunctionSpace(domain, element) # Mixed Function space for all dofs
up = fem.Function(V) # total dofs in the problem/ func for all dofs

u_sp = fem.FunctionSpace(domain, u_el) # u Function space
p_sp = fem.FunctionSpace(domain, p_el) # p Function space
nS_sp = fem.FunctionSpace(domain, nS_el) # nS Function space
nN_sp = fem.FunctionSpace(domain, nN_el) # nN Function space

Vu, up_to_u = V.sub(0).collapse() # dofs in u func subspace
u_n = fem.Function(Vu) # dofs in u func

Vp, up_to_p = V.sub(1).collapse() # dofs in u func subspace
p_n = fem.Function(Vp) # dofs in u func

VnS, up_to_nS = V.sub(2).collapse() # dofs in u func subspace
nS_n = fem.Function(VnS) # dofs in u func

VnN, up_to_nN = V.sub(3).collapse() # dofs in u func subspace
nN_n = fem.Function(VnN) # dofs in u func

u, p, nS, nN = ufl.split(up)
(del_u, del_p, del_nS, del_nN) = ufl.TestFunctions(V)

#---------------------------temporal parameters---------------------------#

# Define temporal parameters
t = 0  # Start time
T = 2000.0  # Final time
num_steps = 1000
dt = T / num_steps  # time step size

#---------------------------Initial Conditions---------------------------#

# Initialize history values
u_n.x.array[:] = 0.0
p_n.x.array[:] = 0.0
nS_n.x.array[:] = 0.0
nN_n.x.array[:] = 0.0

def nS_init(x):
    values = np.zeros((1, x.shape[1]))
    values[0] = 0.5
    return values

def nN_init(x):
    values = np.zeros((1, x.shape[1]))
    values[0] = 0.4
    return values

nN_n.interpolate(nN_init)
nS_n.interpolate(nS_init)

#---------------------------Boundary Conditions---------------------------#

def bc_vec_zero(x):
    return np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)

def bc_scalar_zero(x):
    return np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)

list_bc = []

# Displacement Boundary
u_bc = fem.Function(u_sp)
u_bc.interpolate(bc_vec_zero)

facets = facet_tag.indices[facet_tag.values == 5] 
left_dofs = fem.locate_dofs_topological((V.sub(0), Vu), fdim, facets) 
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))

q_bottom = -4.0
load = fem.Constant(domain, ScalarType(q_bottom))
load_bc = load * ufl.FacetNormal(domain)

#------------------------------Parameters-----------------------------------#

# Material Parameters
mat_k_FOS = 8.3e5
mat_gamma_FR = 1e4
mat_rho_SR = 2.0
mat_rho_LR = 1.0
mat_rho_NR = 2.0
mat_nS_OS = 0.5
mat_nN_OS = 0.4
mat_mu = default_scalar_type(1e4)
mat_lmbda = default_scalar_type(2e3)

# Spatial dimension
d = len(u)

# Identity tensor
I = ufl.variable(ufl.Identity(d))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
b = ufl.variable(F * F.T)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

# Time Derivatives
vel_s = (u - u_n) / dt # Solid Velocity
dnN_dt = (nN - nN_n) / dt # nN's
dnS_dt = (nS - nS_n) / dt # nS's

#Volume Fractions
nF = 1 - nS
nF_OS = 1 - mat_nS_OS

# Darcy
k_F = 8.3e1
w_FS = - (k_F/nF) * (ufl.grad(p)) 
dar = - (k_F) * (ufl.grad(p))

# Stress
fac_1 = (nS/mat_nS_OS)**3
iso1 = (mat_mu/2) * (Ic-3)
iso2 = mat_mu * ufl.ln(J)
iso3 = (mat_lmbda/2) * (ufl.ln(J))**2
T_ef_1 = (iso1 - iso2 + iso3) * I #Ic - 3
T_ef_2 = (mat_mu * (b - I)) + (mat_lmbda * ufl.ln(J) * I)
P = (fac_1 * T_ef_2) - (2.0 * fac_1 * T_ef_1) - (nS * p * I)

# Integration measures
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)
dx = ufl.Measure("dx", domain=domain)

# Weak Form
Mom_B = ufl.inner(P, ufl.grad(del_u)) * dx # Stress
Mom_B += - ufl.inner(del_u, load_bc) * ds(2) # Load

Vol_M = (ufl.div(vel_s) * del_p) * dx # div velocity
Vol_M += - (ufl.inner(dar, ufl.grad(del_p))) * dx # darcy

Vol_S = ((dnS_dt * del_nS) + (nS * ufl.div(vel_s) * del_nS)) * dx # time derivative

Vol_N = (dnN_dt)* del_nN * dx + (nN * ufl.div(vel_s)) * del_nN * dx # time derivative
Vol_N += ufl.inner(ufl.grad(nN), ufl.grad(del_nN)) * dx # artificial
Vol_N += - nN * ufl.inner(w_FS, ufl.grad(del_nN)) * dx # darcy

weak_form = Mom_B + Vol_M + Vol_S + Vol_N

#-------------------------------Solver----------------------------#

# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up, list_bc)

# Initialise newton solver
solver = NewtonSolver(domain.comm, problem)

# Set Newton solver options
solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"
solver.max_it = 10
solver.report = True

# Configure mumps 
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0

for n in range(num_steps):
    t += dt

    # Solve newton-steps
    duration_solve -= MPI.Wtime()
    niter, converged = solver.solve(up)
    duration_solve += MPI.Wtime()

    PETSc.Sys.Print(
        "Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
            t, duration_solve, niter
        )
    )

    u_n.x.array[:] = up.x.array[up_to_u]
    nS_n.x.array[:] = up.x.array[up_to_nS]
    nN_n.x.array[:] = up.x.array[up_to_nN]

Geometry msh file is attached here: split_ver.msh - Google Drive

If you add log.set_log_level(log.LogLevel.INFO) to your code, and change the convergence type to “residual”, you will see:

2023-12-19 07:51:39.484 (   0.379s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 0.103003 (tol = 1e-08) r (rel) = inf(tol = 1e-08)
2023-12-19 07:51:39.773 (   0.668s) [main            ]              petsc.cpp:698   INFO| PETSc Krylov solver starting to solve system.
2023-12-19 07:51:40.109 (   1.004s) [main            ]         TimeLogger.cpp:28    INFO| Elapsed wall, usr, sys time: 0.335990, 0.290000, 0.050000 (PETSc Krylov solver)
2023-12-19 07:51:40.147 (   1.042s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = -nan (tol = 1e-08) r (rel) = -nan(tol = 1e-08)

which means that the first solution of the problem becomes infinite.
There is quite alot going on in your code, and as I’m not an expert when it comes to the set of equations you are solving, I cannot spot a mistake in your code (as it is quite lengthy).

What I can spot is that:

is not used anywhere.

Thank you for the response. Yes, I am not using using scalar zero function. Sorry I should have removed it. One question is that is the way of introducing mixed elements with four solution variables is correct in the code above?

Also, setting nS = 0.2 as initial condition like this makes sense?

u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) # u quad element
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # p linear element
nS_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # nS linear element
element = ufl.MixedElement([u_el, p_el, nS_el]) # Mixed Element

nS_sp = fem.FunctionSpace(domain, nS_el) # nS Function space

VnS, up_to_nS = V.sub(2).collapse() # dofs in u func subspace
nS_n = fem.Function(VnS) # dofs in u func

nS_n.x.array[:] = 0.0

def nS_init(x):
    values = np.zeros((1, x.shape[1]))
    values[0] = 0.2
    return values

nS_n.interpolate(nS_init)

# Time derivative of nS
dnS_dt = (nS - nS_n) / dt # nS's

# Use of time derivative in weak formulation
Vol_S = (dnS_dt + (nS * ufl.div(vel_s))) * del_nS * dx 

I see that after nS_n.interpolate(nS_init) and writing it into a file, nS value is set to 0.2. But on reducing tolerance and solving the system, nS is still 0. There seems to be a problem.

Ho Dokken,

I have narrowed down the problem to the stress equation. If I multiply the factor (nS/0.2) to the stress equation, this problem appears, eg: in the first line below. What would be right way to do this?

P_ef_2 = (nS/0.2) *( (mu_1 * (b - I)) + lmbda_1 * ufl.ln(J) * I)
P = P_ef_2 - nS*p*I

Is it possible to print test function values? For example, I would like to check if the value of nS is going negative or not.

P_ef_2 is not in your original post.

I am assuming it is related to:

fac_1 = (nS/mat_nS_OS)**3
iso1 = (mat_mu/2) * (Ic-3)
iso2 = mat_mu * ufl.ln(J)
iso3 = (mat_lmbda/2) * (ufl.ln(J))**2
T_ef_1 = (iso1 - iso2 + iso3) * I  # Ic - 3
T_ef_2 = (mat_mu * (b - I)) + (mat_lmbda * ufl.ln(J) * I)
P = (fac_1 * T_ef_2) - (2.0 * fac_1 * T_ef_1) - (nS * p * I)

However I cannot see how these map one to one.

Yes it was not there. I was just testing a simpler equation and this was more of an example because I saw that variable nS is a problematic one. That is why I was wondering if it is possible to print the value of nS in this case while it is being solved by the Newton Solver. Until now, I could not figure out how to print the value of nS or any other variable like iso1, T_ef_1 or P_ef_2. Could you please guide me in that direction?

The easiest is to create a custom newton solver, as shown in:
https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html