Modeling Double Descending Thoracic Aortic Aneurysms Using CFD and ResNet-PINNs

A Physics-Informed Neural Network (PINN) framework with ResNet architecture for predicting hemodynamic parameters in double descending thoracic aortic aneurysms (DTAA).

The model learns velocity and pressure fields from CFD-derived wall shear stress (WSS) data while enforcing steady incompressible Navier-Stokes equations in dimensionless form.

Key Highlights

• Clinical Application: Predicts rupture risk in aortic aneurysms
• Performance: R² = 0.9916, RMSE = 0.040 Pa (>99% accuracy)
• Speed: Instant inference after training (no mesh generation required)
• Architecture: ResNet-PINN with self-adaptive loss weighting
• Physics-Informed: Enforces Navier-Stokes equations via automatic differentiation
• Dataset: 202,877 wall surface points from 3 aneurysm geometries

Authors: Noor Fatima, Muhammad Hamza, M. Asif Farooq, Michael Ajao-Olarinoye

Affiliations:

• Department of Mathematics, School of Natural Sciences (SNS), National University of Sciences and Technology (NUST), Islamabad, Pakistan
• Department of Physics, School of Natural Sciences (SNS), National University of Sciences and Technology (NUST), Islamabad, Pakistan
• Computational Science and Mathematical Modelling, Coventry University, Coventry, United Kingdom

Research Overview

This work investigates the effect of aneurysm enlargement on hemodynamic parameters across double descending thoracic aortic aneurysms (DTAA) with variable diameters (4cm, 5cm, 6cm). The PINN framework learns from 202,877 CFD-derived wall surface points with WSS measurements, achieving R² = 0.9916 and RMSE = 0.040 Pa. The model uses ResNet architecture with self-adaptive loss weighting, Random Fourier Features, and gradient-based WSS computation from learned velocity fields.

Quick Start

Paper Configuration

The following configuration reproduces the results reported in the paper:

python main.py --mode both --data-dir data --epochs 6000 --lr 1e-3 \
  --hidden-dim 128 --num-layers 10 \
  --scheduler-step-size 1000 --scheduler-gamma 0.998

Evaluation Only

To evaluate a pre-trained model:

python main.py --mode evaluate --model-path results/models/aneurysm_flow_pinn.pth --data-dir data --hidden-dim 128 --num-layers 10

Methodology

Governing Equations

The model solves the steady incompressible Navier-Stokes equations in dimensionless form:

Continuity:

∇·u* = 0

Momentum:

u*·∇u* + ∇p* - (1/Re)∇²u* = 0

where u* = (u*, v*, w*) is the dimensionless velocity, p* is dimensionless pressure, and Re is the Reynolds number.

Non-Dimensionalization

All quantities are transformed to dimensionless form using reference scales computed exclusively from training data (not from arbitrary choices like inlet velocities). This physics-based approach ensures:

• Dimensional consistency across all terms in governing equations
• Improved numerical conditioning (variables maintain O(1) magnitudes)
• Scale-invariant learning relationships

Reference Scales:

Scale	Definition	Physical Meaning
L_ref	mean([max(X)-min(X), max(Y)-min(Y), max(Z)-min(Z)])	Characteristic length: average spatial extent
τ̄	mean(WSS_training)	Mean wall shear stress from training data
U_ref	(τ̄ × L_ref) / μ	Characteristic velocity derived from WSS
WSS_ref	μ × U_ref / L_ref = τ̄	Reference WSS (equals mean by construction)
D_ref	mean(diameters_training)	Mean aneurysm diameter
Re	ρ × U_ref × L_ref / μ	Reynolds number (typically 100-1000)
P_scale	WSS_ref × L_ref	Pressure scale (for momentum balance)

Transformation Formulas:

# Spatial coordinates (centered by centroid)
x* = (x - x_mean) / L_ref
y* = (y - y_mean) / L_ref
z* = (z - z_mean) / L_ref
# Result: x*, y*, z* ∈ approximately [-0.5, 0.5]

# Velocity and pressure
u* = u / U_ref
v* = v / U_ref
w* = w / U_ref
p* = p / (ρ × U_ref²)

# Wall shear stress
τ_w* = τ_w / WSS_ref

# Diameter (conditioning parameter)
d* = d / D_ref

Physical Constants (Blood Properties):

• Density (ρ): 1060 kg/m³
• Dynamic viscosity (μ): 0.0035 Pa·s (Newtonian approximation)
• Gravitational acceleration (g): 9.81 m/s²

Key Advantage: By anchoring the velocity scale to WSS measurements, the dimensionless velocity field maintains direct physical relevance to the target predictions, improving training stability and convergence.

Wall Shear Stress Computation

WSS is derived from velocity gradients projected onto the tangent plane at wall boundaries:

τ_w = μ ||(I - nn^T) S n||

where S = ½(∇u + ∇u^T) is the strain rate tensor and n is the inward-pointing surface normal. This approach ties WSS predictions directly to the learned velocity field rather than using a separate output head.

Network Architecture

Four independent subnetworks for (u, v, w, p) with:

Random Fourier Features (Tancik et al. 2020):

γ(x*, y*, z*) = [x*, y*, z*, sin(2πB·[x*,y*,z*]), cos(2πB·[x*,y*,z*])]
where B ∈ ℝ^(K×3) ~ N(0, σ²), K=12, σ=1.0

Learnable Swish Activation (Ramachandran et al. 2017):

φ_β(x) = x · σ(β·x),  where σ(x) = 1/(1 + e^(-x))

Layer Normalization (Ba et al. 2016):

LayerNorm(h) = γ ⊙ (h - μ)/√(σ² + ε) + β,  ε = 10⁻⁵

Architecture: Input(28) → Linear(128) → 10×ResidualBlocks(128) → OutputHeads(1 each for u*, v*, w*, p*). Total parameters: ~1.7M.

Loss Function

The total loss combines five components with self-adaptive multiplicative weighting (McClenny & Braga-Neto, 2023):

L_total = λ_WSS × L_WSS + λ_boundary × L_boundary
        + λ_cont × L_continuity + λ_mom × L_momentum
        + γ_p × Var(p*)

where: λ_i = exp(log λ_i) are learnable multiplicative weights

Loss Components

Component	Formula	Purpose
L_WSS	(1/N_w) Σ (τ_w,pred - τ_w,CFD)²	Data fidelity: match CFD wall shear stress
L_boundary	(1/N_w) Σ (u*² + v*² + w*²)	No-slip condition at walls
L_continuity	(1/N_c) Σ (∂u*/∂x* + ∂v*/∂y* + ∂w*/∂z*)²	Incompressibility (mass conservation)
L_momentum	(1/N_c) Σ (R_x² + R_y² + R_z²)	Navier-Stokes momentum balance
Var(p)*	γ_p × Var(p*)	Pressure gauge stability (γ_p = 10⁻⁴)

Momentum Residuals (Dimensionless Form):

R_x = u*·∂u*/∂x* + v*·∂u*/∂y* + w*·∂u*/∂z* + ∂p*/∂x* - (1/Re)∇²u*
R_y = u*·∂v*/∂x* + v*·∂v*/∂y* + w*·∂v*/∂z* + ∂p*/∂y* - (1/Re)∇²v*
R_z = u*·∂w*/∂x* + v*·∂w*/∂y* + w*·∂w*/∂z* + ∂p*/∂z* - (1/Re)∇²w*

Where:

• First term: convective acceleration (inertial effects)
• Second term: pressure gradient (driving force)
• Third term: viscous diffusion (molecular transport)

Self-Adaptive Weighting Mechanism

Key Innovation: Both network parameters θ and loss weights {log λ_i} are jointly optimized via AdamW:

{θ, log λ_i}^(t+1) = AdamW({θ, log λ_i}^(t), ∇L_total)

Advantages of Multiplicative Weighting:

1. Automatic Balance: Network discovers optimal trade-offs without manual tuning
2. Positive Weights: exp(·) transformation guarantees λ_i > 0
3. Dynamic Adaptation: If physics residuals grow large, corresponding weights decrease to restore balance
4. Floor Protection: Minimum weight floor (default 0.1) prevents physics terms from decaying to zero

Initial Values (Learned During Training):

Weight	Initial log λ	Initial λ = exp(log λ)	Rationale
λ_WSS	0.8	2.23	Emphasize data fidelity from start
λ_boundary	-0.4	0.67	Moderate no-slip enforcement
λ_continuity	-0.3	0.74	Balanced mass conservation
λ_momentum	-0.4	0.67	Allow focus on WSS while maintaining physics

These weights automatically evolve during training to achieve optimal convergence.

Collocation Point Generation

Physics constraints are enforced at strategically sampled collocation points generated by offsetting wall points along inward surface normals. This approach focuses physics enforcement where it matters most while avoiding arbitrary interior sampling.

Sampling Strategy

Collocation Point: x_coll = x_wall - d × n_inward

Where:

• x_wall: wall boundary point with known WSS
• n_inward: inward-pointing unit normal (estimated via Open3D or PCA)
• d: offset distance defining proximity to wall

Band Type	Distance from Wall	Points Ratio	Purpose
Near-wall	2-6 mm	40%	Capture steep velocity gradients in boundary layer
Interior	10-50 mm	60%	Enforce PDE in bulk flow region and recirculation zones

Rationale:

• Near-wall focus: WSS computation requires accurate velocity gradients close to walls
• Interior coverage: Momentum balance must hold throughout the domain
• Adaptive density: More points near complex flow regions (aneurysm necks, domes)

Surface Normal Estimation

Two methods available for computing inward-pointing unit normals:

1. Open3D Method (Recommended):

# Radius-based neighborhood with adaptive density
r_i = 3.0 × median(nearest_neighbor_distances)
neighbors = points within r_i of x_i (max 64 points)

# PCA on local neighborhood
C_i = (1/|neighbors|) Σ (x_j - μ_i)(x_j - μ_i)^T
n_i = eigenvector corresponding to smallest eigenvalue

# Orient inward toward global centroid
if (c_global - x_i)·n_i < 0:
    n_i = -n_i

2. PCA Method (Fallback):

# k-nearest neighbors (default k=10)
neighbors = k_nearest_neighbors(x_i)

# Same PCA + orientation as above

Why Open3D?

• Handles irregular sampling patterns robustly
• Adapts to local point density variations
• More stable normals in high-curvature regions (aneurysm domes)

Batch Sampling During Training

Each training epoch samples:

• Wall batch: 1024 points with WSS supervision + surface normals
• Collocation batch: 1024 points for physics constraints
  • 410 near-wall points (2-6mm from wall)
  • 614 interior points (10-50mm from wall)

This ensures balanced supervision between data fidelity (wall) and physics consistency (interior).

Command-Line Interface

Core Arguments

Argument	Default	Description
--mode	train	Execution mode: train, evaluate, or both
--data-dir	data	Path to CFD data directory
--output-dir	results	Output directory for results
--model-path	-	Path to saved model (for evaluation)
--config	-	Path to YAML configuration file
--device	auto	Computing device: auto, cpu, or cuda
--seed	-	Random seed for reproducibility
--log-level	INFO	Logging verbosity
--no-visualizations	False	Disable plot generation

Network Architecture

Argument	Default	Description
--hidden-dim	64	Hidden layer width
--num-layers	5	Number of residual blocks
--activation	swish	Activation function: swish or tanh
--weight-init	glorot	Weight initialization: glorot, kaiming, or default

Fourier Features

Argument	Default	Description
--use-fourier	enabled	Enable Fourier feature embeddings
--no-use-fourier	-	Disable Fourier feature embeddings
--fourier-num-frequencies	12	Number of random frequencies (K)
--fourier-sigma	1.0	Std for random frequency matrix

Training

Argument	Default	Description
--epochs	2000	Maximum training epochs
--lr	1e-3	Learning rate
--batch-size	1024	Wall points batch size
--scheduler-step-size	500	StepLR step size (epochs)
--scheduler-gamma	0.97	StepLR decay factor
--eval-batch-size	1024	Evaluation batch size

Surface Normals

Argument	Default	Description
--normals-method	open3d	Normal estimation: open3d or pca
--normals-k	10	k for PCA k-NN
--o3d-radius-multiplier	3.0	Open3D radius multiplier
--o3d-max-nn	64	Open3D maximum neighbors
--o3d-orient-k	50	Open3D orientation k

Collocation

Argument	Default	Description
--near-bands	"0.002,0.006"	Near-wall bands (dimensionless)
--deep-bands	"0.01,0.02;0.03,0.05"	Interior bands
--coll-ratio-near	0.4	Fraction for near-wall points
--coll-ratio-deep	0.6	Fraction for interior points

Citation

If you use this code in your research, please cite:

@article{fatima2025pinn,
  title={Modeling of Double Descending Thoracic Aortic Aneurysms Using Computational Fluid Dynamics (CFD) and Residual-Based Physics-Informed Neural Networks (ResNets-PINNs)},
  author={Fatima, Noor and Hamza, Muhammad and Farooq, M. Asif and Ajao-Olarinoye, Michael},
  journal={Pramana - Journal of Physics},
  year={2026 },
  publisher={Indian Academy of Sciences}
}

Contact

For questions or collaborations, please contact:

• Noor Fatima: noorfatima.nf801@gmail.com
• Michael Ajao-Olarinoye: olarinoyem@uni.coventry.ac.uk