Unable to solve linear system using PETSc Krylov solver for incompressible flow simulation

Hi folks,

I am simulating the incompressible flow in a cylinder. Consider the axisymmetric, I created the 2D model. In my model, the geometry is a rectangle which is the axisymmetric part of a cylinder. There are left, top, bottom, and right boundaries. The left boundary is the axisymmetric boundary where dot(u, n) and shear stress are zero. On all other boundaries u=(0,0). So sigma * n * v = 0 on ds and there is no ds term in the variational form. The x-axis in the model is actually r-axis and y-axis is the z-axis. I convert the nabla_grad and div operators to cylindrical coordinates and dv should be 2pirdrdz=x[0]dx as the 2*pi are eliminated in the variational form.

The model cannot converge after a few time steps as the solved velocity is super high (~10e103). I tried a lot of different initial conditions but they were not helpful. Any suggestions would be appreciated.

My DOLFIN version is 2018.1.0 and the error message says:

*** -------------------------------------------------------------------------
*** 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 linear system using PETSc Krylov solver.
*** Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = 0.000000e+00).
*** Where: This error was encountered inside PETScKrylovSolver.cpp.
*** Process: 0


*** DOLFIN version: 2018.1.0
*** Git changeset:
*** -------------------------------------------------------------------------

My code is:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Apr 21 10:50:16 2023

@author: xin
"""

from fenics import *
from mshr import *

import matplotlib.pyplot as plt
import numpy as np

tt = 0.01                 # Simulation time
num_steps = 10      # Number of time steps
dt = tt/num_steps      # Time step size
k = Constant(dt)

Lv = 0.4                # Length of the vaporizer 0.4m
La = 1.2                # Length of the adiabatic section 1.2m
Lc = 0.4                # Length of the condenser 0.4m
tot_l = Lv+La+Lc        # Total length
h = 0#0.02                # Thickness of the tube
r = 0.1                 # Radius of heat pipe

# Define material properties

# viscosity
mu = 1.34*10**(-5)

rho = 0.5542              # Density of vapor

# Create domain
whole_domain = Rectangle(Point(0,0),Point(h+r,tot_l+2*h))             # Whole domain

tol = 1E-14                   # Tolorence used to distinguish different domain and define boundaries

# Define boundaries
class top(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], tot_l + 2*h)
class bot(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 0)
class left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0],0)

# Right boundary of the condenser
class right_con(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], r+h) and x[1]>= tot_l +h - Lc

# Right boundary of the adiabatic section
class right_adi(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], r+h) and x[1]<= tot_l -Lc + h and x[1]>=Lv + h

# Right boundary of the vaporizer
class right_vap(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], r+h) and x[1]<=Lv + h

top = top()
bot = bot()
left = left()
right_vap = right_vap()
right_adi = right_adi()
right_con = right_con()

# Generate Mesh
mesh = generate_mesh(whole_domain, 150)

x = SpatialCoordinate(mesh)

# normal vector
n = FacetNormal(mesh)

# Mark the tube domain as 0 and the fluid domain as 1
# material = MeshFunction('size_t', mesh, 2, mesh.domains())

# vtkfile = File('meshf.pvd')
# vtkfile<<material

# Mark the boundary
bf = MeshFunction('size_t', mesh,1)
bf.set_all(0)
top.mark(bf, 1)
bot.mark(bf, 2)
left.mark(bf, 3)
right_con.mark(bf, 4)
right_adi.mark(bf, 5)
right_vap.mark(bf, 6)

vtkfile = File('meshf_f.pvd')
vtkfile<<bf

# Create function space
V = VectorFunctionSpace(mesh, 'P', 2)    # for velosity      
Q = FunctionSpace(mesh, 'CG', 1)          # for pressure, phase field, and temp

# Define test and trial functions

# For velocity
u = TrialFunction(V)                 
v = TestFunction(V)

# For Pressure
P = TrialFunction(Q)
q = TestFunction(Q)

#Define functions
u_n = Function(V)    # save the result at n time step
u_ = Function(V)     # save the tentative result


P_n = Function(Q)
P_ = Function(Q)


# define operations in cylindrical coordinates

# define gradient in cylindrical coordintes
def nabla_grad_vector_innertensor(v):
    return as_tensor([[v[0].dx(0),0, v[0].dx(1)],
                      [0,v[0]/x[0] ,0],
                          [v[1].dx(0),0, v[1].dx(1)]])
# nabala grad for inner dot of tensor and vector
def nabla_grad_vector_innervector(v):
    return as_tensor([[v[0].dx(0), v[0].dx(1)],
                          [v[1].dx(0), v[1].dx(1)]])

def nabla_grad_scalar(v):
    return as_vector((v.dx(0),v.dx(1)))

def div_vector(v):
    return (1/x[0])*(x[0]*v[0]).dx(0) + v[1].dx(1)

# create XDMF files for visualization output
file_u = File('heatpipe/velocity.pvd')
file_p = File('heatpipe/pressure.pvd')

# initialize the velocity field
u_n = interpolate(Expression(('0','6.0*x[0]*(0.41-x[0])/pow(0.41,2)'), degree=2), V)

# define the variational form for incompressible flow (splitting method)

#define body force
# gravity density
f = Constant((0,-9.8))

U = 1/2*(u_n+u)

# define bc for velosity

bcu_top = DirichletBC(V, Constant((0,0)),top)
bcu_bot = DirichletBC(V, Constant((0,0)),bot)
bcu_left = DirichletBC(V.sub(0), Constant(0),left)
bcu_right_con = DirichletBC(V, Constant((0,0)),right_con)
bcu_right_adi = DirichletBC(V, Constant((0,0)),right_adi)
bcu_right_vap = DirichletBC(V, Constant((0,0)),right_vap)

bcu = [bcu_top,bcu_left,bcu_bot,bcu_right_con,bcu_right_adi,bcu_right_vap]


# Define strain tensor
def epsilon(u):
    return sym(nabla_grad_vector_innertensor(u))
# Define Stress tensor
def sigma(u,P):
    return 2*mu*epsilon(u) - P*Identity(3)

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*x[0]*dx \
+ rho*dot(dot(u_n, nabla_grad_vector_innervector(u_n)), v)*x[0]*dx \
+ inner(sigma(U, P_n), epsilon(v))*x[0]*dx \
- dot(f, v)*x[0]*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad_scalar(P), nabla_grad_scalar(q))*x[0]*dx
L2 = dot(nabla_grad_scalar(P_n), nabla_grad_scalar(q))*x[0]*dx - (1/k)*div_vector(u_)*q*x[0]*dx

# Define variational problem for step 3
a3 = dot(u, v)*x[0]*dx
L3 = dot(u_, v)*x[0]*dx - k*dot(nabla_grad_scalar(P_ - P_n), v)*x[0]*dx

# Assemble matrices 
# since there is no parameters or function need to be updated at the left hand side, so we assemble it once
    # step 1: solve the tentative velocity
A1 = assemble(a1)  # pho need to be updated after each time step
[bc.apply(A1) for bc in bcu]
A2 = assemble(a2)
A3 = assemble(a3)

# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"

# time steping
t = 0
for n in range(num_steps):
# for n in range(1):
    print(n)
    t = t+dt
    
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(),b1,'gmres','default')
    
    # step 2: pressure correction step
    b2 = assemble(L2)
    solve(A2, P_.vector(),b2,'cg',prec)
    
    # step 3: velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(),b3,'gmres','default')
    
    # update previous solution
    u_n.assign(u_)
    P_n.assign(P_)
    
    # save solution to files
    file_u<<(u_,t)
    file_p<<(P_,t)