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