3D Patient Specific

[mhadidy]

I have this geometry which is a patient specific artery with one inlet and two outlets.

I’m trying to run 3D steady Navier stocks to get velocity, I tried to thing one of them is the following nonlinear forum, and also tried the IPC method.

But none of them are converging, the nonlinear method kind of stuck at

and the error almost flattened,

Newton iteration 24: r (abs) = 7.332e-05 (tol = 1.000e-05) r (rel) = 3.905e-06 (tol = 1.000e-08)

this is my script for this one:

from __future__ import print_function
import os
import time
import numpy as np
from fenics import *

# ==============================================================================
# FEniCS SETTINGS
# ==============================================================================
parameters["ghost_mode"] = "shared_facet"
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2
set_log_level(LogLevel.PROGRESS)

comm = MPI.comm_world
rank = MPI.rank(comm)

def log(msg):
    if rank == 0:
        print(msg, flush=True)

# ==============================================================================
# USER PARAMETERS
# ==============================================================================
MU  = 0.04
RHO = 1.06
NU  = MU / RHO

INLET_ID   = 1
OUTLET_IDS = [2, 3]
WALL_ID    = 4

INLET_U_TARGET = 57.0
NUM_STEPS = 500

MESH_FILE   = "Mesh/mesh.xdmf"
FACET_FILE  = "Mesh/facet_tags.xdmf"
OUTPUT_FILE = "results/velocity_steady.xdmf"

os.makedirs("results", exist_ok=True)

# ==============================================================================
# HELPER FUNCTIONS
# ==============================================================================
def compute_inlet_normal(coords, domain_center):
    log("[NORMAL] Computing inlet normal via PCA")

    inlet_center = coords.mean(axis=0)
    X = coords - inlet_center
    C = np.dot(X.T, X) / float(X.shape[0])
    eigvals, eigvecs = np.linalg.eigh(C)

    n = eigvecs[:, np.argmin(eigvals)]
    n /= np.linalg.norm(n)

    if np.dot(n, domain_center - inlet_center) < 0:
        n = -n
        log("[NORMAL] Flipped normal to point INTO domain")

    log(f"[NORMAL] Normal vector = {n}")
    return n, inlet_center


class ParabolicInlet(UserExpression):
    def __init__(self, center, radius, normal, Umax, **kwargs):
        self.c = np.array(center, dtype=float)
        self.R = float(radius)
        self.n = np.array(normal, dtype=float)
        self.Umax = float(Umax)
        super().__init__(**kwargs)

    def value_shape(self):
        return (3,)

    def eval(self, values, x):
        d = np.array(x) - self.c
        rvec = d - np.dot(d, self.n) * self.n
        r = np.linalg.norm(rvec)

        u_mag = 0.0
        if r <= self.R:
            u_mag = self.Umax * (1.0 - (r / self.R) ** 2)

        values[:] = u_mag * self.n

# ==============================================================================
# LOAD MESH
# ==============================================================================
log("====================================================")
log("[STEP] Loading mesh")

mesh = Mesh()
with XDMFFile(MESH_FILE) as f:
    f.read(mesh)

log(f"[MESH] Number of cells    : {mesh.num_cells()}")
log(f"[MESH] Number of vertices : {mesh.num_vertices()}")

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim() - 1)
with XDMFFile(FACET_FILE) as f:
    f.read(mvc, "f")
facets = MeshFunction("size_t", mesh, mvc)

log(f"[FACETS] Unique facet IDs found: {np.unique(facets.array())}")

# ==============================================================================
# EXTRACT INLET GEOMETRY
# ==============================================================================
log("====================================================")
log("[STEP] Extracting inlet geometry")

inlet_facets = np.where(facets.array() == INLET_ID)[0]
log(f"[INLET] Number of inlet facets: {len(inlet_facets)}")

vertex_ids = []
for f in inlet_facets:
    vertex_ids.extend(Facet(mesh, int(f)).entities(0))

vertex_ids = np.unique(vertex_ids)
coords = mesh.coordinates()[vertex_ids]

if coords.shape[0] < 3:
    raise RuntimeError("[ERROR] Not enough inlet vertices")

domain_center = mesh.coordinates().mean(axis=0)
normal, inlet_center = compute_inlet_normal(coords, domain_center)

radius = np.max(np.linalg.norm(coords - inlet_center, axis=1))

log(f"[INLET] Center = {inlet_center}")
log(f"[INLET] Radius = {radius}")

# ==============================================================================
# FUNCTION SPACES
# ==============================================================================
log("====================================================")
log("[STEP] Creating function spaces")

V_el = VectorElement("CG", mesh.ufl_cell(), 2)
P_el = FiniteElement("CG", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V_el, P_el]))

(v, q) = TestFunctions(W)
w = Function(W)
(u, p) = split(w)

log(f"[DoF] Velocity DoFs : {W.sub(0).dim()}")
log(f"[DoF] Pressure DoFs : {W.sub(1).dim()}")

# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
log("====================================================")
log("[STEP] Defining boundary conditions")

inlet_expr = ParabolicInlet(
    center=inlet_center,
    radius=radius,
    normal=normal,
    Umax=INLET_U_TARGET,
    degree=2
)

bcs = [
    DirichletBC(W.sub(0), inlet_expr, facets, INLET_ID),
    DirichletBC(W.sub(0), Constant((0, 0, 0)), facets, WALL_ID),
    DirichletBC(W.sub(1), Constant(0.0), facets, OUTLET_IDS[0]),
    DirichletBC(W.sub(1), Constant(0.0), facets, OUTLET_IDS[1]),
]

# ==============================================================================
# VARIATIONAL FORM
# ==============================================================================
log("====================================================")
log("[STEP] Building variational form")

nu = Constant(NU)

F = (
    dot(dot(u, nabla_grad(u)), v) * dx
    + nu * inner(grad(u), grad(v)) * dx
    - div(v) * p * dx
    - q * div(u) * dx
)

dw = TrialFunction(W)
dF = derivative(F, w, dw)

problem = NonlinearVariationalProblem(F, w, bcs, dF)
solver  = NonlinearVariationalSolver(problem)

prm = solver.parameters
prm["nonlinear_solver"] = "newton"
prm["newton_solver"]["absolute_tolerance"] = 1e-5
prm["newton_solver"]["relative_tolerance"] = 1e-8
prm["newton_solver"]["maximum_iterations"] = 30
prm["newton_solver"]["linear_solver"] = "gmres"
prm["newton_solver"]["preconditioner"] = "hypre_amg"

# "gmres"      # General nonsymmetric (most robust)
# "cg"         # Conjugate Gradient (SPD only)
# "bicgstab"   # Nonsymmetric, faster but less robust
# "tfqmr"      # Similar to BiCGStab
# "minres"     # Symmetric indefinite
# "richardson" # Simple iteration (rarely useful alone)

# prm["newton_solver"]["preconditioner"] = "jacobi"

# "jacobi"     # Diagonal scaling
# "sor"        # Successive over-relaxation
# "ilu"        # Incomplete LU (default good choice)

# "icc"        # Incomplete Cholesky (SPD only)
# "asm"        # Additive Schwarz (domain decomposition)
# "hypre_amg"  # Algebraic multigrid (very strong)
# "petsc_amg"  # PETSc AMG

# Scability
# prm["newton_solver"]["linear_solver"] = "gmres"
# prm["newton_solver"]["preconditioner"] = "hypre_amg"



prm["newton_solver"]["krylov_solver"]["relative_tolerance"] = 1e-7

# ==============================================================================
# SOLVER LOOP (PER-VELOCITY RAMP TIMING)
# ==============================================================================
log("====================================================")
log("[STEP] Starting Newton continuation solves")

xdmf = XDMFFile(OUTPUT_FILE)
xdmf.parameters["rewrite_function_mesh"] = False
xdmf.parameters["flush_output"] = True

U_vals = np.linspace(
    (1.0 / NUM_STEPS) * INLET_U_TARGET,
    INLET_U_TARGET,
    NUM_STEPS
)

global_start_time = time.time()

for i, U in enumerate(U_vals):
    inlet_expr.Umax = float(U)

    ramp_start_time = time.time()

    log("----------------------------------------------------")
    log(f"[RAMP] {i+1}/{NUM_STEPS} | Inlet Umax = {U:.6f}")
    log("[RAMP] Newton solve starting...")

    solver.solve()

    ramp_time = time.time() - ramp_start_time
    total_time = time.time() - global_start_time

    log(f"[TIME] Ramp {i+1} time           : {ramp_time:.2f} s")
    log(f"[TIME] Total elapsed so far     : {total_time:.2f} s")

    if i == NUM_STEPS - 1:
        u_sol, p_sol = w.split(deepcopy=True)
        u_sol.rename("velocity", "velocity")
        xdmf.write(u_sol, i)

# ==============================================================================
# DONE
# ==============================================================================
log("====================================================")
log("[DONE] Steady Navier–Stokes solve completed")
log("====================================================")

while the IPC method it’s giving me this error

*** -------------------------------------------------------------------------
*** 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: 9e1777c0df7c762c7230189c87c13fd5ad33a4f0
*** -------------------------------------------------------------------------

srun: error: notch482: task 0: Exited with exit code 1

Below is my IPC Script

#!/usr/bin/env python3
# coding: utf-8

from __future__ import print_function
import os
import time
import numpy as np
from fenics import *

# ==============================================================================
# GLOBAL LOGGING HELPERS
# ==============================================================================
comm = MPI.comm_world
rank = MPI.rank(comm)

def log(msg):
    if rank == 0:
        print(msg, flush=True)

def banner(msg):
    if rank == 0:
        print("\n" + "="*80)
        print(msg)
        print("="*80 + "\n", flush=True)

# ==============================================================================
# FEniCS SETTINGS
# ==============================================================================
banner("FEniCS INITIALIZATION")

parameters["ghost_mode"] = "shared_facet"
parameters["form_compiler"]["optimize"] = True
parameters["form_compiler"]["quadrature_degree"] = 2
set_log_level(LogLevel.PROGRESS)

log(f"[MPI] Rank {rank} / {MPI.size(comm)} active")

# ==============================================================================
# USER PARAMETERS
# ==============================================================================
banner("USER PARAMETERS")

MU  = 0.04
RHO = 1.06
NU  = MU / RHO

INLET_ID   = 1
OUTLET_IDS = [2, 3]
WALL_ID    = 4

INLET_U_TARGET = 57.0
N_RAMP = 100

DT = 1e-4
T_MAX = 10.0
MAX_STEPS = int(np.ceil(T_MAX / DT))

STEADY_TOL = 1e-6
STEADY_STREAK = 300

SAVE_EVERY  = 500
PRINT_EVERY = 50

log(f"MU={MU}, RHO={RHO}, NU={NU}")
log(f"INLET_U_TARGET={INLET_U_TARGET}")
log(f"N_RAMP={N_RAMP}")
log(f"DT={DT}, MAX_STEPS={MAX_STEPS}")
log(f"STEADY_TOL={STEADY_TOL}, STEADY_STREAK={STEADY_STREAK}")

# ==============================================================================
# FILES
# ==============================================================================
banner("FILES & DIRECTORIES")

MESH_FILE   = "Mesh/mesh.xdmf"
FACET_FILE  = "Mesh/facet_tags.xdmf"

OUT_DIR     = "results_unsteady_ramp_VERBOSE"
U_TRANSIENT = os.path.join(OUT_DIR, "velocity_unsteady.xdmf")
U_STEADY    = os.path.join(OUT_DIR, "velocity_steady.xdmf")

os.makedirs(OUT_DIR, exist_ok=True)

log(f"MESH_FILE   = {MESH_FILE}")
log(f"FACET_FILE  = {FACET_FILE}")
log(f"OUT_DIR     = {OUT_DIR}")

# ==============================================================================
# LOAD MESH
# ==============================================================================
banner("LOADING MESH")

mesh = Mesh()
with XDMFFile(comm, MESH_FILE) as f:
    f.read(mesh)

log(f"[MESH] Cells    : {mesh.num_cells()}")
log(f"[MESH] Vertices : {mesh.num_vertices()}")
log(f"[MESH] Dim      : {mesh.geometry().dim()}")

mvc = MeshValueCollection("size_t", mesh, mesh.topology().dim()-1)
with XDMFFile(comm, FACET_FILE) as f:
    f.read(mvc, "f")
facets = MeshFunction("size_t", mesh, mvc)

if rank == 0:
    log(f"[FACETS] Unique IDs: {np.unique(facets.array())}")

# ==============================================================================
# INLET GEOMETRY (PCA NORMAL)
# ==============================================================================
banner("INLET GEOMETRY EXTRACTION")

inlet_facet_ids = np.where(facets.array() == INLET_ID)[0]
log(f"[INLET] Number of inlet facets: {len(inlet_facet_ids)}")

vertex_ids = []
for fid in inlet_facet_ids:
    vertex_ids.extend(Facet(mesh, int(fid)).entities(0))
vertex_ids = np.unique(vertex_ids)

coords = mesh.coordinates()[vertex_ids]
log(f"[INLET] Unique inlet vertices: {coords.shape[0]}")

domain_center = mesh.coordinates().mean(axis=0)

X = coords - coords.mean(axis=0)
C = np.dot(X.T, X) / X.shape[0]
eigvals, eigvecs = np.linalg.eigh(C)
normal = eigvecs[:, np.argmin(eigvals)]
normal /= np.linalg.norm(normal)

if np.dot(normal, domain_center - coords.mean(axis=0)) < 0:
    normal *= -1
    log("[INLET] Normal flipped to point INTO domain")

radius = np.max(np.linalg.norm(coords - coords.mean(axis=0), axis=1))

log(f"[INLET] Center = {coords.mean(axis=0)}")
log(f"[INLET] Radius = {radius}")
log(f"[INLET] Normal = {normal}")

# ==============================================================================
# FUNCTION SPACES
# ==============================================================================
banner("FUNCTION SPACES")

V = VectorFunctionSpace(mesh, "CG", 2)
Q = FunctionSpace(mesh, "CG", 1)

log(f"[DoF] Velocity DoFs = {V.dim()}")
log(f"[DoF] Pressure DoFs = {Q.dim()}")

# Trial/Test
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Solution fields
u_n   = Function(V, name="u_n")
u_t   = Function(V, name="u_tent")
u_new = Function(V, name="u")

p_n = Function(Q, name="p_n")
p_  = Function(Q, name="p")

# ==============================================================================
# INLET PROFILE
# ==============================================================================
banner("INLET PROFILE DEFINITION")

class ParabolicInlet(UserExpression):
    def __init__(self, center, radius, normal, Umax, **kwargs):
        self.c = np.array(center)
        self.R = radius
        self.n = normal
        self.Umax = Umax
        super().__init__(**kwargs)

    def value_shape(self):
        return (3,)

    def eval(self, values, x):
        d = np.array(x) - self.c
        rvec = d - np.dot(d, self.n) * self.n
        r = np.linalg.norm(rvec)
        if r <= self.R:
            values[:] = self.Umax * (1 - (r/self.R)**2) * self.n
        else:
            values[:] = (0,0,0)

inlet_expr = ParabolicInlet(
    center=coords.mean(axis=0),
    radius=radius,
    normal=normal,
    Umax=0.0,
    degree=2
)

log("[INLET] Parabolic profile initialized (Umax=0)")

# ==============================================================================
# BOUNDARY CONDITIONS
# ==============================================================================
banner("BOUNDARY CONDITIONS")

bcu = [
    DirichletBC(V, inlet_expr, facets, INLET_ID),
    DirichletBC(V, Constant((0,0,0)), facets, WALL_ID),
]

bcp = [DirichletBC(Q, Constant(0.0), facets, i) for i in OUTLET_IDS]

log(f"[BC] Inlet ID  = {INLET_ID}")
log(f"[BC] Wall ID   = {WALL_ID}")
log(f"[BC] Outlet IDs = {OUTLET_IDS}")

# ==============================================================================
# VARIATIONAL FORMS (IPCS)
# ==============================================================================
banner("VARIATIONAL FORMS (IPCS)")

dt = Constant(DT)
mu = Constant(MU)
rho = Constant(RHO)
nrm = FacetNormal(mesh)

def epsilon(u): return sym(nabla_grad(u))
def sigma(u,p): return 2*mu*epsilon(u) - p*Identity(3)

U_mid = 0.5*(u_n + u)

F1 = (
    rho*dot((u-u_n)/dt, v)*dx
    + rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
    + inner(sigma(U_mid, p_n), epsilon(v))*dx
    - dot(mu*nabla_grad(U_mid)*nrm, v)*ds
)

a1, L1 = lhs(F1), rhs(F1)
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/dt)*div(u_t)*q*dx
a3 = dot(u, v)*dx
L3 = dot(u_t, v)*dx - dt*dot(nabla_grad(p_ - p_n), v)*dx

log("[FORMS] Assembling matrices...")

A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

for bc in bcu: bc.apply(A1)
for bc in bcp: bc.apply(A2)

log("[FORMS] Matrices assembled and BCs applied")

# ==============================================================================
# OUTPUT
# ==============================================================================
banner("OUTPUT FILES")

xdmf = XDMFFile(comm, U_TRANSIENT)
xdmf.parameters["flush_output"] = True
xdmf.parameters["rewrite_function_mesh"] = False

xdmf_steady = XDMFFile(comm, U_STEADY)
xdmf_steady.parameters["flush_output"] = True
xdmf_steady.parameters["rewrite_function_mesh"] = False

# ==============================================================================
# TIME LOOP
# ==============================================================================
banner("STARTING TIME LOOP (RAMP → STEADY)")

u_prev = u_n.vector().copy()
steady_count = 0
t = 0.0
start_time = time.time()

for step in range(MAX_STEPS):

    inlet_expr.Umax = (
        INLET_U_TARGET*(step+1)/N_RAMP if step < N_RAMP else INLET_U_TARGET
    )

    # Step 1
    b1 = assemble(L1)
    for bc in bcu: bc.apply(b1)
    solve(A1, u_t.vector(), b1, "gmres", "hypre_amg")

    # Step 2
    b2 = assemble(L2)
    for bc in bcp: bc.apply(b2)
    solve(A2, p_.vector(), b2, "gmres", "hypre_amg")

    # Step 3
    b3 = assemble(L3)
    solve(A3, u_new.vector(), b3, "cg", "sor")

    # Steady detection
    du = u_new.vector().copy()
    du.axpy(-1.0, u_prev)

    rel = du.norm("l2") / max(u_new.vector().norm("l2"), 1e-30)

    if step >= N_RAMP and rel < STEADY_TOL:
        steady_count += 1
    else:
        steady_count = 0

    if step % PRINT_EVERY == 0:
        log(f"[STEP {step:6d}] "
            f"Uin={inlet_expr.Umax:7.3f} | "
            f"rel={rel:.2e} | "
            f"streak={steady_count:4d} | "
            f"u_max={u_new.vector().get_local().max():.3e}")

    if step % SAVE_EVERY == 0:
        xdmf.write(u_new, t)

    if steady_count >= STEADY_STREAK:
        banner("STEADY STATE REACHED")
        log(f"Step = {step}")
        log(f"Relative change = {rel}")
        break

    u_prev = u_new.vector().copy()
    u_n.assign(u_new)
    p_n.assign(p_)
    t += DT

# ==============================================================================
# FINAL OUTPUT
# ==============================================================================
banner("WRITING FINAL STEADY FIELD")

u_new.rename("velocity", "velocity")
xdmf_steady.write(u_new, t)

elapsed = time.time() - start_time
log(f"[DONE] Total wall time = {elapsed:.2f} s")
banner("SIMULATION COMPLETE")

I also got the mesh from Simvascular and converted it into .xdmf

[dokken]

@mhadidy

parameters["form_compiler"]["quadrature_degree"] = 2
set_log_level(LogLevel.PROGRESS)

You are severely underintegrating most of the terms in your variational forms. I would remove the pinning of the quadrature degree.

@mhadidy

prm["newton_solver"]["linear_solver"] = "gmres"
prm["newton_solver"]["preconditioner"] = "hypre_amg

To first ensure that you are solving the problem, use an lu solver/direct rather than an iterative solver.

Furthermore, you are not specifying which solver that fails in the splitting scheme approach, and without being able to run the code, as you haven’t shared the mesh, it’s hard to give a more detailed review of the code.

[mhadidy]

Hi @dokken
Thanks for your fast replay.

I tried unpinning the quadarature degree and still diverging.

For the usage of solving direct solver I keep getting Out of memory even I used the cluster with 300 gb of ram.

for the IPC method this is the script the one that which falling is the tentaive velocity step “the first one” but I tried printing the velocity magnitude and it’s start exploding after the first iteration till it fails at the 3rd iteration.

I will share the mesh, it’s not very fine but I’m using it to diagnose.

Mesh.zip

Thanks!