PETSc solver fails with DIVERGED_NANORINF for Navier-Stokes in left ventricle geometry

General
1

Hi everyone,

I’m simulating blood flow in a 2D geometry resembling the left ventricle of the heart using the IPCS method in FEniCS 2019.2.0. The simulation runs fine for a while but eventually fails with:

RuntimeError: Unable to solve linear system using PETSc Krylov solver.
Reason: Solution failed to converge in 0 iterations (PETSc reason DIVERGED_NANORINF, residual norm ||r|| = -1.000000e+00).

I understand this usually means there’s a NaN or Inf in the matrix or RHS vector, but I’m not sure how it’s introduced in my case. Below is the full code for context:

from fenics import *
from mshr import *
import numpy as np

output_dir = 'results/'
T = 3  #Three cardiac cycles
num_steps = 50000
dt = T / num_steps  # time step size
mu =0.035 # 0.035          # dynamic viscosity
rho = 1.06            # density


# Definir geometría ventrículo
a = 8.0  # Semieje mayor
b = 2  # Semieje menor

#Geometria 
rectangle = Rectangle(Point(-a, -b), Point(0, b))
ellipse = Ellipse(Point(0, 0), a, b)
semi_ellipse = ellipse - rectangle

# Rectángulos para los tubos (válvula aórtica y válvula mitral)
aorta_rect = Rectangle(Point(-2, -1.8), Point(0, -1.8 + 2.4))
mitral_rect = Rectangle(Point(-2, 1), Point(0, 1.8))

# Combinar geometría
semi_ellipse_with_tubes = semi_ellipse + aorta_rect + mitral_rect

 #Generar malla
mesh_resolution = 64 #Resolución
mesh = generate_mesh(semi_ellipse_with_tubes, mesh_resolution)


basis_order = 1 
V = VectorFunctionSpace(mesh, 'P', basis_order)
Q = FunctionSpace(mesh, 'P', 1)


boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)  # 

# Recorres solo las facetas de frontera (exterior)
for facet in facets(mesh):
    if facet.exterior():
        boundaries[facet] = 3  # Paredes por defecto

# Contorno (definimos boundaries): funciones devuelven true si es un punto del contorno considerado
# Contorno de entrada (valvula mitral)
def inflow_mitral_boundary(x, on_boundary):
    return near(x[0], -2, 1e-5) and (-1.8 - 1e-5) <= x[1] <= (0.6 + 1e-5)

def outflow_aorta_boundary(x, on_boundary):
    return near(x[0], -2, 1e-5) and (1 - 1e-5) <= x[1] <= (1.8 + 1e-5)

inflow = AutoSubDomain(inflow_mitral_boundary)
outflow = AutoSubDomain(outflow_aorta_boundary)

# Sobrescribimos las marcas para inflow y outflow
inflow.mark(boundaries, 1)  # 1 para inflow mitral
outflow.mark(boundaries, 2) # 2 para outflow aorta

# Exporta el resultado
File('boundaries.pvd') << boundaries

# class Inflow(SubDomain):
#   def inside(self, x, on_boundary):
#     return near(x[0], 0)
# class Outflow(SubDomain):
#   def inside(self, x, on_boundary):
#     return near(x[0], 3.8)
# class Walls(SubDomain):
#   def inside(self, x, on_boundary):
#     return on_boundary and x[1] >= 0.4 or near(x[1], 0)
# 
# inflow = Inflow()
# outflow = Outflow()
# walls = Walls()
#         
# boundaries = MeshFunction('size_t', mesh, 1)
# boundaries.set_all(0)
# inflow.mark(boundaries, 1)
# outflow.mark(boundaries, 2)
# walls.mark(boundaries, 3)

# File(output_dir + 'boundaries.pvd') << boundaries         

# Define inflow profile
# inflow_profile = ('4.0*1.5*x[1]*(0.81 - x[1]) / pow(0.81, 2)', '0')

#inflow_profile = ('32', '0')

# Define boundary conditions
#bcu_inflow = DirichletBC(V, Expression(inflow_profile, degree = 2), inflow)
bcu_walls = DirichletBC(V, Constant((0., 0.)),boundaries,3)
#bcu = [bcu_inflow, bcu_walls]
bcp_outflow = DirichletBC(Q, Constant(0.), boundaries, 2)
bcp = [bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)


u_n = Function(V) 
u_  = Function(V)
p_n = Function(Q) 
p_  = Function(Q)

# Define expresion
U  = 0.5*(u_n + u)
n  = FacetNormal(mesh)
f  = Constant((0, 0))
k  = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)

# Define symmetric gradient
def epsilon(u):
   return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
   return 2*mu*epsilon(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
  + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
  + inner(sigma(U, p_n), epsilon(v))*dx \
  + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
  - dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
#[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Create XDMF files 
xdmffile_u = XDMFFile(output_dir + 'velocity.xdmf')
xdmffile_p = XDMFFile(output_dir  + 'pressure.xdmf')

file_u = File("RESULTADOS/navier_stokes_ventriculo_iteracionparabolicosinusoidal/velocidad.pvd")
file_p = File("RESULTADOS/navier_stokes_ventriculo_iteracionparabolicosinusoidal/presion.pvd")


# Create time series (for use in reaction_system.py)
timeseries_u = TimeSeries(output_dir + 'velocity_series')
timeseries_p = TimeSeries(output_dir  + 'pressure_series')

# Save mesh to file (for use in reaction_system.py)
#File(output_dir + 'aneurysm.xml.gz') << mesh


# Create progress bar
#progress = Progress('Looping', num_steps)
#set_log_level(LogLevel.PROGRESS)

# Time-stepping
t = 0
f=1
A=25
for n in range(num_steps):
  
  v_inlet_BC = A * np.abs(np.sin(np.pi * f * t) )

  inflow_profile = Expression(
    ("v_inlet_BC * (1.0 - pow((x[1]+ 0.6)/1.2, 2))","0.0"),  # ejemplo para flujo en eje x[0]
    degree=2,
    v_inlet_BC= v_inlet_BC,
    R=1.2
)
  bcu_inflow = DirichletBC(V, inflow_profile, boundaries, 1)
  bcu = [bcu_inflow, bcu_walls]

  # Apply boundary conditions to matrices
  [bc.apply(A1) for bc in bcu]

  # Step 1: Tentative velocity step
  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)
  [bc.apply(b2) for bc in bcp]
  solve(A2, p_.vector(), b2, 'gmres', 'default')

  # Step 3: Velocity correction step
  b3 = assemble(L3)
  solve(A3, u_.vector(), b3, 'cg', 'sor')

  # Save solution to file (XDMF/HDF5)
  if n % 500 == 0 and t > 2 :
    file_u << (u_, t)
    file_p << (p_, t)
    
    #   xdmffile_p.write(p_, t)
      # Save nodal values to file
    #timeseries_u.store(u_.vector(), t)
      #timeseries_p.store(p_.vector(), t)

  # Update previous solution
  u_n.assign(u_)
  p_n.assign(p_)


  # Update current time
  t += dt

  #print 'u max:', u_.vector().get_local().max()
  if n % 50 == 0:
    if MPI.rank(MPI.comm_world) == 0:  
      print('u max:', u_.vector().get_local().max(),flush=True)
      print('t', t,flush=True)
      print(f"Paso {n}, t = {t:.3f} s, v_inlet_BC = {v_inlet_BC:.3f} cm/s")

I’m using the IPCS fractional-step method with:

  • dt = 3 / 50000
  • Viscosity mu = 0.035, density rho = 1.06
  • Sinusoidal parabolic inflow
  • Geometry: semi-ellipse + rectangles for mitral and aortic valves

The matrix A1 is assembled only once before the loop and reused (with BCs reapplied). At each time step I update the inflow boundary condition.

I’d like to ask:

  • Any idea what might be causing the solver to diverge?
  • More importantly: is there a better numerical scheme you’d recommend that can handle more turbulent flows or be more stable within FEniCS?
    *In future steps I want to move/deform the ventricle geometry. What’s the recommended way to remesh or adapt the mesh in FEniCS during time-stepping? Any examples or libraries that help with that?
    Thanks so much in advance!
2

In your solution profiles, do you happen to notice any backflow just before the convergence fails? That is, inflow at the outflow boundary. That could cause your issue. The solution would be to add a ‘directional do-nothing’ boundary condition at the outflow.

3

Hi, thanks for the suggestion regarding the directional do-nothing condition at the outflow. I’ve tried to implement it in my FEniCS 2019 setup, but unfortunately I haven’t managed to make it work properly yet.

What I attempted is modifying the variational formulation in the tentative velocity step (F1) to only apply the stress boundary term at the outflow (ds(2)) when the flow is outgoing, using conditional(...) and ge(dot(U, n), 0) like this:

U = 0.5*(u_n + u_)

F1 = rho*dot((u - u_n) / k, v)*dx \
  + rho*dot((3/2)*dot(u_n, nabla_grad(u_n)) - (1/2)*dot(u_n_1, nabla_grad(u_n_1)), v)*dx \
  + inner(sigma(U, p_n), epsilon(v))*dx \
  + conditional(ge(dot(U, n), 0), p_n * dot(n, v), 0) * ds(2) \
  - conditional(ge(dot(U, n), 0), mu * dot(as_vector(nabla_grad(U) * n), v), 0) * ds(2) \
  - dot(f, v)*dx

On the other hand I also changed the variational formulation from the original to see if the instability (which appears when using low viscosity values) could be alleviated, but unfortunately the simulation still breaks at the same time step, when the sinusoidal inflow velocity starts to decrease at the inlet.

I’m solving the unsteady incompressible Navier–Stokes equations using a second-order time-stepping scheme:

F1 = rho*dot((u - u_n) / k, v)*dx \
  + rho*dot((3/2)*dot(u_n, nabla_grad(u_n))-(1/2)*dot(u_n_1,nabla_grad(u_n_1)), v)*dx \
  + inner(sigma(U, p_n), epsilon(v))*dx \
  + dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
  - dot(f, v)*dx
  • Adams–Bashforth for the convective term,
  • Crank–Nicolson for the viscous term,
  • and a standard IPCS fractional step method.

If anyone has a minimal working example or any suggestion to correctly implement a directional do-nothing condition that avoids backflow (especially in FEniCS 2019), I’d really appreciate it.

Also, would this kind of instability be easier to resolve using FEniCSx? I’m considering migrating my code but not sure if it would help in this particular case.

Thanks again for your time and help!