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
, densityrho = 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!