Enhancing Arterial Blood Flow Simulations through Physics-Informed Neural Networks
Shivam Bhargava Nagaiah Chamakuri




1]Shivam Bhargava
shivam19@iisertvm.ac.in
1,2]Nagaiah Chamakuri
nagaiah.chamakuri@iisertvm.ac.in
[1]School of Mathematics, IISER Thiruvananthapuram, 695551 Kerala, India.
[2]Center for High-Performance Computing, IISER Thiruvananthapuram, 695551 Kerala, India


This study introduces a computational approach leveraging physics-informed neural networks (PINNs) for the efficient computation of arterial blood flows, particularly focusing on solving the incompressible Navier-Stokes equations by using the domain decomposition technique. Unlike conventional computational fluid dynamics methods, PINNs offer advantages by eliminating the need for discretized meshes and enabling the direct solution of partial differential equations (PDEs). In this paper, we propose the weighted Extended Physics-Informed Neural Networks (WXPINNs) and weighted Conservative Physics-Informed Neural Networks (WCPINNs), tailored for detailed hemodynamic simulations based on generalized space-time domain decomposition techniques. The inclusion of multiple neural networks enhances the representation capacity of the weighted PINN methods. Furthermore, the weighted PINNs can be efficiently trained in parallel computing frameworks by employing separate neural networks for each sub-domain. We show that PINNs simulation results circumvent backflow instabilities, underscoring a notable advantage of employing PINNs over traditional numerical methods to solve such complex blood flow models. They naturally address such challenges within their formulations. The presented numerical results demonstrate that the proposed weighted PINNs outperform traditional PINNs settings, where sub-PINNs are applied to each subdomain separately. This study contributes to the integration of deep learning methodologies with fluid mechanics, paving the way for accurate and efficient high-fidelity simulations in biomedical applications, particularly in modeling arterial blood flow.



Physics-Informed Neural Networks, Hemodynamic Simulations, Incompressible Navier-Stokes, Domain decomposition, weighted XPINNs and weighted CPINNs. 





§ INTRODUCTION

The study of arterial blood flow occupies a pivotal junction in cardiovascular research, bridging advanced theoretical fluid dynamics with pivotal clinical applications. This interdisciplinary research is fundamental for deciphering the complex dynamics governing blood movement through arterial pathways, contributing significantly to our comprehension of the cardiovascular system. The precise simulation of arterial blood flow is not only critical for advancing theoretical understanding but also instrumental in improving diagnostic and therapeutic strategies in cardiovascular medicine <cit.>.

Central to the simulation of arterial blood flow is the incompressible Navier-Stokes equations, which present formidable challenges due to their nonlinear nature and the complexity of arterial geometries. In recent decades, researchers have extensively utilized various computational techniques to simulate and analyze blood flow, aiming to comprehend the correlation between vascular diseases and hemodynamics <cit.>. Tokuda et al. <cit.> employed the finite element method to numerically simulate blood flow, with a focus on understanding stroke mechanisms during cardio-pulmonary bypass. However, such numerical simulations of hemodynamics pose significant computational challenges in terms of time and memory requirements. To address this issue, numerous approaches have been explored to expedite the resolution of these problems <cit.>. Such numerical methods for solving these equations fall short for numerous medical applications, like swiftly calculating hemodynamics for emergency cases and providing real-time surgical guidance. Moreover, in clinical practice, employing computational fluid dynamics necessitates repetitive simulations for varying patients, adding to the physician's workload. Given these challenges, there's a growing need to seek non-invasive, precise, cost-effective, and computationally efficient methods for capturing cardiovascular flow dynamics. Data-driven deep learning algorithms offer a promising avenue to address these issues. This necessitates the exploration of novel computational methodologies that can offer both accuracy and efficiency.

In this context, Physics-Informed Neural Networks (PINNs) emerge as a groundbreaking approach, melding the predictive power of machine learning with the foundational principles of physics <cit.>. PINNs leverage the universal approximation capabilities of neural networks, constrained by the governing physical laws, to solve differential equations efficiently <cit.>. This research innovates by applying PINNs to the domain of arterial blood flow, tackling the Navier-Stokes equations across varied geometries with unprecedented computational efficiency and model robustness. Many investigations have investigated employing deep learning for predicting hemodynamics. These studies encompass forecasting hemodynamic parameters or clinical metrics, local blood flow velocity, vessel cross-sectional area, and blood pressure <cit.>. Until now, most of the literature concerning DL-driven flow modeling predominantly utilizes fully-connected neural network architectures <cit.>. Only a handful of publications explore the application of more intricate network structures, such as convolutional neural networks <cit.>. In computational fluid dynamics (CFD) simulations, encountering backflow at open boundaries poses a frequent challenge, often resulting in unphysical oscillations and instability issues, even at moderate Reynolds numbers, see <cit.> and referenced therein. Traditional numerical methods require stabilization terms to address these backflow instabilities.
Conversely, PINNs offer a different approach. They don't rely on stabilization terms to mitigate backflow instabilities and inherently tackle such challenges within their formulations.

The primary objective of this work is to demonstrate the efficacy of PINNs in simulating arterial blood flow, overcoming existing computational barriers while ensuring high fidelity in the simulations. For extensive simulations of Partial Differential Equations (PDEs), the high training cost associated with Physics-Informed Neural Networks (PINNs) makes tackling large-scale PDEs inherently expensive, impacting their efficiency compared to traditional numerical methods. Thus, there's a critical need to expedite the convergence of these models without compromising performance. Domain decomposition techniques, widely utilized in conventional numerical methods like finite difference, finite volume, and finite element methods, offer a promising solution. Here, the computational domain is subdivided into multiple subdomains, with interactions occurring solely at shared boundaries where continuity conditions are enforced. Within the realm of deep learning frameworks, the application of domain decomposition approaches in PINN is explored, notably in the conservative PINN (cPINN) method <cit.> for conservation laws. Beyond cPINN, other strategies within the Scientific Machine Learning (SciML) domain include employing local neural networks on partitioned subdomains, as demonstrated by Li et al.<cit.>, who leveraged the variational principle. Similarly, Kharazmi et al. <cit.> proposed a variational PINN framework. Moreover, eXtended PINNs (XPINNs) <cit.> offer a comprehensive solution for solving various PDEs, exhibiting advantages akin to cPINN, such as employing separate neural networks in each subdomain, efficient hyper-parameter adjustment, facile parallelization, and substantial representation capacity. Further refining of the simulation accuracy and efficiency was studied in <cit.>.
In this study, we introduce weighted XPINN and CPINN methodologies to tackle hemodynamics model equations, aiming to enhance computational efficiency through domain decomposition strategies. By leveraging a finite number of subdomains, massively parallel computation becomes feasible, empowering effective handling of large-scale problems through domain decomposition. These weighted PINN approaches are more generalizations of the cPINNs <cit.> and XPINNs<cit.>. Through a meticulous application of ADAM <cit.> and L-BFGS <cit.> optimization techniques and domain decomposition strategies, this research enhances the scalability and robustness of PINNs. These contributions mark significant advancements in integrating machine learning with fluid mechanics, offering a new paradigm for high-fidelity simulations in biomedical engineering and beyond.

The subsequent sections of this paper are structured as follows: Section 2 presents mathematical models concerning blood dynamics, encompassing the Navier-Stokes equations. Section 3 elaborates on computational methodologies, with particular emphasis on Physics-Informed Neural Networks (PINNs) and their diverse iterations. The application and refinement of these models in simulating blood flow dynamics are presented in Section 4. Following that, Section 5 unveils the outcomes and evaluations of the simulation performance. Lastly, Section 6 encapsulates the findings, addresses limitations, and delineates potential avenues for future research, highlighting the fusion of machine learning with fluid mechanics.



§ MATHEMATICAL FORMULATION

Blood displays viscoelastic properties and non-Newtonian behavior primarily at low shear rates due to the aggregation and alignment of erythrocytes, which enhance apparent viscosity and induce shear-thinning properties. Despite these complex behaviors, blood can be effectively modeled as a Newtonian fluid within arterial flow regimes, where shear rates are substantially higher. This model simplifies due to the uniform alignment of erythrocytes under high shear conditions, thereby simplifying the computational model without significantly sacrificing the accuracy of flow predictions in arterial geometries.


The dynamics of blood flow within the arterial system are governed by the incompressible Navier-Stokes equations, formulated within the computational domain Ω⊂ℝ^2 and expressed over the simulation interval (0, T]. These equations capture the physiological nuances of arterial blood flow <cit.>:

    ρ( u_t + 𝐮·∇𝐮) - ∇ p + μ∇^2 𝐮   = 0   in   (0, T] ×Ω, 
    ∇·𝐮   = 0   in   (0, T] ×Ω,

where 𝐮 is the velocity field, p is the pressure field, ρ = 1060  kg/m^3 is the blood density, and μ = 3.5 × 10^-3 Pa·s is the dynamic viscosity, assuming a Newtonian fluid approximation. T denotes the terminal simulation time.


Incorporating the Cauchy stress tensor σ, the equations can be alternatively formulated as:

    u_t + 𝐮·∇𝐮 = 1/ρ∇·σ,


with σ = -p𝐈 + μ (∇𝐮 + (∇𝐮)^T), and the pressure p derived from:

    p = -1/2tr(σ).


For a two-dimensional flow, where velocity components are u and v, the stress tensor σ is given by:

    σ = 
        [   -p + 2μ u_x μ (u_y + v_x); μ (v_x + u_y)   -p + 2μ v_y ],


Initial and boundary conditions:
To initiate the simulations, the velocity field throughout the computational domain is set to a quiescent state:

    𝐮(x, y, 0) = [ 0; 0 ],   ∀ (x, y) ∈Ω.

This simplifies the analysis of fluid dynamics as the simulation progresses from a state of rest.

Boundary conditions are designed to mimic physiological conditions realistically and are adapted to the specific geometry and expected flow characteristics of each test case:

Inlet conditions:
At all inlet boundaries, the velocity is defined by a function that incorporates both the geometric specifics of the domain and the time-dependent aspects of physiological flow, often modulating according to a function that mimics cardiac pulsatility:

    𝐮(x_in, y, t) = U_max f(y, t) [ 1; 0 ],

where x_in represents the inlet boundary coordinate, U_max is a scaling factor for maximum velocity, and f(y, t) is a profile function tailored to the specific domain geometry and flow conditions.

Outlet and wall conditions:
At the outlet, a condition is set to simulate a specific pressure condition depending on the study's need:

    p(x_out, y) = p_out,

where x_out is the outlet boundary coordinate, and p_out is typically set to zero or another value reflecting ambient or downstream pressure conditions.


No-slip boundary conditions are enforced along the walls of the domain to reflect the physical reality that the fluid velocity at a solid boundary is zero due to viscosity:

    𝐮(x, y_wall) = [ 0; 0 ],

where y_wall denotes the coordinates of the wall boundary, applicable to both straight and curved segments in any domain shape.

These boundary conditions are crafted to ensure that the flow dynamics within the simulation accurately reflect the complex interactions expected in real-world fluid flow scenarios, especially within physiological contexts. They provide a robust framework that can adapt to both the specific details in the test cases developed in Section 4. 



§ COMPUTATIONAL METHODOLOGY


Physics-Informed Neural Networks (PINNs) <cit.> are designed to predict solutions 𝐮(𝐱, t) for a set of partial differential equations (PDEs) across a specified domain Ω and time interval [0, T], where 𝐱 represents spatial coordinates and t total simulation time. The neural network, denoted by 𝒩_θ(𝐱, t), with parameters θ, is tasked with approximating the solution to the PDEs subject to initial and boundary conditions.

Given a differential operator 𝒟, a function f representing source terms, boundary conditions described by a function g, and initial conditions defined by h, the problem can be formally expressed as:

    𝒟[𝐮](𝐱, t)    = f(𝐱, t),    𝐱   ∈Ω,   t ∈ [0, T], 
    ℬ[𝐮](𝐱, t)    = g(𝐱, t),    𝐱   ∈∂Ω,   t ∈ [0, T], 
    𝐮(𝐱, 0)    = h(𝐱),    𝐱   ∈Ω.


The core of the PINN methodology lies in the formulation of a physics-informed loss function, which quantifies the network's deviation not only from empirical data but also from the physical laws encapsulated by the governing differential equations. This loss function is constructed by evaluating the residuals associated with the governing PDEs (R_g), boundary conditions (R_bc), and initial conditions (R_ic) across the domain and its boundaries. These residuals are defined as follows:

    R_g(𝐱, t; θ)    = 𝒟[𝒩_θ](𝐱, t) - f(𝐱, t), 
    
    R_bc(𝐱, t; θ)    = ℬ[𝒩_θ](𝐱, t) - g(𝐱, t), 
    
    R_ic(𝐱; θ)    = 𝒩_θ(𝐱, 0) - h(𝐱).


The aggregated loss function for the weighted PINN, combining these residuals, is given by:

    ℒ_WPINN(θ) = 1/N_g∑_j=1^N_g R_g(𝐱_j, t_j; θ) ^2, + β (1/N_bc∑_k=1^N_bc R_bc(𝐱_k, t_k; θ) ^2 + 1/N_ic∑_l=1^N_ic R_ic(𝐱_l; θ) ^2)

where N_g, N_bc, and N_ic denote the number of collocation points for the governing equations, boundary conditions, and initial conditions, respectively. Here β>0 is taken as a user-defined weighting coefficient that balances ℒ_g (governing equations loss) and ℒ_bc/ic (boundary and initial conditions loss) and accelerates convergence, enabling prioritization based on their physical relevance to enhance model accuracy and stability. 
Through the minimization of the above physics-informed loss, weighted PINNs effectively learn the dynamics of the system directly from the governing equations, offering a powerful tool for solving complex physical problems where traditional numerical methods may face limitations.

The techniques XPINNs<cit.> and CPINNs<cit.> advance the capabilities of PINNs through domain decomposition. The main idea behind this methodology is partitioning the computational domain into multiple subdomains, Ω_i, where each with a dedicated sub-network, 𝒩_θ_i, that facilitates localized solutions of complex PDEs.

For each subdomain Ω_i, a specific sub-PINN 𝒩_θ_i is trained to approximate the local solution 𝐮_i. This process adheres to the governing differential equations, boundary conditions, and initial conditions particular to Ω_i:

    𝒟_i[𝐮_i](𝐱, t; θ_i)    = f_i(𝐱, t),    𝐱∈Ω_i,   t ∈ [0, T], 
    ℬ_i[𝐮_i](𝐱, t; θ_i)    = g_i(𝐱, t),    𝐱∈∂Ω_i,   t ∈ [0, T], 
    𝐮_i(𝐱, 0; θ_i)    = h_i(𝐱),    𝐱∈Ω_i,   t = 0.

These equations ensure the model captures the local physics accurately within each subdomain.

For XPINNs, ensuring smooth transitions between subdomains involves calculating an averaged solution at the interface points, Γ_ij, between adjacent subdomains Ω_i and Ω_j. This averaging is represented as:


    𝐮_avg,ij(𝐱, t) = 1/2( 𝒩_θi(𝐱, t) + 𝒩_θ_j(𝐱, t) ).


The interface residual, R_interface,i, measures the difference at these points between the sub-PINN's prediction for Ω_i and the averaged solution, ensuring interface continuity:


    R_interface,i(𝐱, t; θ_i, θ_j) = 𝒩_θ_i(𝐱, t) - 𝐮_avg,ij(𝐱, t).


Minimizing R_interface,i across interfaces achieves a cohesive global solution across the computational domain.

In weighted XPINNs (WXPINNs), the physics-informed loss for each subdomain includes weighted terms for governing equations, boundary conditions, initial conditions, and interfaces. The overall loss function is given by:

    ℒ_WXPINN(θ_i) =     1/N_g,i∑_j=1^N_g,i R_g,i^2 + β( 1/N_bc,i∑_k=1^N_bc,i R_bc, i^2 + 1/N_ic,i∑_l=1^N_ic,i R_ic,i^2 ) 
        + γ∑_j ∈𝒩(i)1/N_ij∑_𝐱∈Γ_ij R_interface,i^2.


where N_ij is the number of interface points between subdomains Ω_i and Ω_j, and 𝒩(i) represents the neighboring subdomains of Ω_i. Here, β>0 serves as a user-defined weighting coefficient that tackles the balance of boundary and initial conditions loss, and γ>0 serves as a user-defined weighting coefficient that targets the enhancement of interface continuity. γ controls the emphasis on the interface loss component, ℒ_interface, which quantifies discrepancies and ensures smooth transitions at the inter-domain boundaries between adjacent sub-networks. Adjusting γ allows for fine-tuning the model's capacity to seamlessly integrate solutions across these interfaces, which is pivotal for achieving a coherent and accurate global solution. 

Through the minimization of the above physics-informed loss, weighted XPINNs solve complex PDEs across large and geometrically intricate domains by addressing them piecemeal yet in a manner that preserves the integrity and continuity of the overall solution.

 

In CPINNs, conservation across subdomain interfaces is ensured by aligning flux vectors, 𝐅_i and 𝐅_j, from adjacent subdomains Ω_i and Ω_j. The conservation condition at these interfaces, Γ_ij, is expressed as:


    ( 𝐅_i(𝐱, t; θ_i) - 𝐅_j(𝐱, t; θ_j) ) ·𝐧_ij = 0  ,

where 𝐧_ij denotes the unit normal vector at Γ_ij, pointing from Ω_i to Ω_j. This equation ensures the net flux of conserved quantities across the interface is zero, upholding the essential physical conservation laws.
The flux residual, quantifying the adherence to conservation laws at the interface Γ_ij, is defined as:

    R_flux,ij(𝐱, t; θ_i, θ_j) = ( 𝐅_i(𝐱, t; θ_i) - 𝐅_j(𝐱, t; θj) ) ·𝐧_ij ,


In weighted CPINNs, the physics-informed loss for each subdomain incorporates residuals from the governing equations, boundary conditions, initial conditions, interfaces, and conservation of flux:

    ℒ_WCPINN(θ_i) =     1/N_g,i∑_j=1^N_g,i R_g,i^2 + β (1/N_bc,i∑_k=1^N_bc,i R_bc,i^2  + 1/N_ic,i∑_l=1^N_ic,i R_ic,i^2) 
        + γ∑_j ∈𝒩(i)1/N_ij∑_𝐱∈Γ_ij R_interface,i^2 + δ∑_j ∈𝒩(i)1/N_ij∑_𝐱∈Γ_ij R_flux,i^2


Here, β>0 and γ>0 serve as user-defined weighting coefficients that tackle the balance of boundary and initial conditions loss and tackle interface continuity, respectively, and δ>0 acts as a user-defined weighting coefficient that emphasizes the conservation of flux at the interfaces between subdomains. This coefficient is critical for ensuring that key conservation properties, such as mass and momentum in the case of the incompressible Navier-Stokes equations, are maintained across subdomain boundaries. By tuning δ, the model can more effectively handle fluid dynamics problems where the accurate representation of flow and pressure continuity is essential for achieving realistic and physically accurate simulations. Adjusting δ allows the model to precisely manage how these conservation laws influence the overall solution, ensuring that the continuity of physical quantities is preserved.

Through the minimization of the above physics-informed loss, weighted CPINNs solve complex PDEs across large and geometrically intricate domains by addressing them piecemeal yet in a manner that preserves the integrity and continuity of the overall solution.





§ MODEL IMPLEMENTATION


To apply weighted Physics-Informed Neural Networks (WPINNs) to the incompressible Navier-Stokes equations, a neural network that maps the spatiotemporal variables {t, 𝐱}^T to the mixed-variable solution {ψ, p, σ} is constructed. The input 𝐱 has 2 components (x, y). The output σ has three components σ=(σ_11, σ_12,σ_22)^T.



The stream function ψ is used rather than velocity 𝐮 directly to ensure that the divergence-free condition of the flow is maintained. The velocity 𝐮 = (u', v') where u' and v' are predicted velocities can be computed as:


    u' = ψ_y,    v' = -ψ_x


Thus the residuals obtained from u and v can be computed as:


    [ R_u; R_v ] = 
    [ ρ u'_t + ρ (u' u'_x + v' u'_y) - (-p_x + 2μ u'_xx) - (μ (u'_y + v'_x)_y); ρ v'_t + ρ (u' v'_x + v' v'_y) - (μ (v'_x + u'_y)_x) - (-p_y + 2μ v'_yy) ]


The residual obtained from p for predicted pressure p' can be computed as:


    R_p = p' + 1/2(σ_11 + σ_22)


The residual obtained from σ_11, σ_12, σ_22 can be computed as:


    [ R_σ_11; R_σ_12; R_σ_22 ] = [   (-p' + 2μ u'_x) - σ_11; (μ u'_y + μ v'_x) - σ_12;   (-p' + 2μ v'_y) - σ_22 ]


The residual from the governing equations is given by:

    R_g(x, y, t)^2 = R_u^2 + R_v^2 + R_p^2 + R_σ_11^2 + R_σ_12^2 + R_σ_22^2


The residuals from the boundary and initial conditions are given by:

    R_bc(x, y, t)^2    = u(x, y, t) - u'(x, y, t)^2 + v(x, y, t) - v'(x, y, t)^2 
       + p(x, y, t) - p'(x, y, t)^2 
    R_ic(x, y, 0)^2    = u(x, y, 0) - u'(x, y, 0)^2 + v(x, y, 0) - v'(x, y, 0)^2 
       + p(x, y, 0) - p'(x, y, 0)^2


The physics-informed loss and initial/boundary condition loss can now be given by:

    ℒ_g   = 1/N_g∑_i=1^N_g R_g(x_i, y_i, t_i) ^2 
    ℒ_bc/ic   = 1/N_bc∑_i=1^N_bc R_bc(x_i, y_i, t_i) ^2 + 1/N_ic∑_i=1^N_ic R_ic(x_i, y_i, 0) ^2

where N_(·) represents the number of collocation points.

The loss function for the Physics-Informed Neural Network (PINN) model is formulated as follows:

    ℒ_WPINN = ℒ_g + βℒ_bc/ic,

where β serves as a weighting factor to balance the components of the loss function.

To apply the weighted XPINN to the incompressible Navier-Stokes equations, the domain Ω is decomposed into M non-overlapping subdomains {Ω_i}_i=1^M, each with its own sub-PINN 𝒩_i. The interface between adjacent subdomains Ω_i and Ω_j is denoted as Γ_i,j where j ∈{i-1, i+1}. For each subdomain Ω_i, interface conditions are enforced at the shared boundaries with adjacent subdomains. The predicted velocity fields u_i' and v_i' and the pressure field p_i' from the sub-PINN for Ω_i, and u_j', v_j', and p_j' from Ω_j are used to define the interface residuals.

The interface residual for velocity and pressure for subdomain Ω_i at the interface Γ_i,j is given by:


    R_interface,u(x, y, t)    = u_i'(x, y, t) - u_i'(x, y, t) + u_j'(x, y, t)/2, 
    
    R_interface,v(x, y, t)    = v_i'(x, y, t) - v_i'(x, y, t) + v_j'(x, y, t)/2, 
    
    R_interface,p(x, y, t)    = p_i'(x, y, t) - p_i'(x, y, t) + p_j'(x, y, t)/2,


The interface loss ℒ_interface,ij for subdomain Ω_i for the interface Γ_i,j is calculated based on the interface residuals:


    0.9ℒ_interface,ij = 1/N_ij∑_k=1^N_ij(  R_interface,u(x_k, y_k, t_k) ^2 +  R_interface,v(x_k, y_k, t_k) ^2 +  R_interface,p(x_k, y_k, t_k) ^2 )


Thus the interface loss ℒ_interface,i for subdomain Ω_i becomes:

    ℒ_interface,i = ℒ_interface,i(i+1) + ℒ_interface,i(i-1) .


The total loss function for each sub-PINN 𝒩_i in subdomain Ω_i is a combination of the physics-informed loss, boundary, and initial condition loss, and the interface loss:


    ℒ_WXPINN,i =  ℒ_g,i + βℒ_bc/ic,i + γℒ_interface,i

where γ serves as a weighting factor to balance the components of the loss function.


To apply the weighted CPINN to the incompressible Navier-Stokes equations, the domain Ω is decomposed into M non-overlapping subdomains {Ω_i}_i=1^M, each with its own sub-PINN 𝒩_i. The interface between adjacent subdomains Ω_i and Ω_j is denoted as Γ_i,j where j ∈{i-1, i+1}. For each subdomain Ω_i, interface conditions are enforced at the shared boundaries with adjacent subdomains. The predicted velocity fields u_i' and v_i' and the pressure field p_i' from the sub-PINN for Ω_i, and u_j', v_j', and p_j' from Ω_j are used to define the flux residuals for conservation of mass and conservation of momentum.

The flux residuals for mass and momentum for subdomain Ω_i at the interface Γ_i,j are given by:


    R_flux,m(x, y, t)    = ρ( u_i'(x, y, t) - u_j'(x, y, t) ), 
    
    R_flux,μ(x, y, t)    = ρ( u_i'(x, y, t)^2 + p_i'(x, y, t) ) - ρ( u_j'(x, y, t)^2 + p_j'(x, y, t) ),


The flux loss ℒ_f,ij for subdomain Ω_i for the interface Γ_i,j is calculated based on the flux residuals:


    ℒ_flux,ij = 1/N_ij∑_k=1^N_ij(  R_flux,m(x_k, y_k, t_k) ^2 +  R_flux,μ(x, y, t)(x_k, y_k, t_k) ^2 )


Thus the flux loss ℒ_flux,i for subdomain Ω_i becomes:

    ℒ_flux,i = ℒ_flux,i(i+1) + ℒ_flux,i(i-1)


The total loss function for each sub-PINN 𝒩_i in subdomain Ω_i is a combination of the physics-informed loss, boundary, and initial condition loss, interface loss, and the flux loss:


    ℒ_WCPINN,i =  ℒ_g,i + βℒ_bc/ic,i + γℒ_interface,i + δℒ_flux, i


where δ serves as a weighting factor to balance the components of the loss function.


This study encompasses two primary test cases designed to demonstrate the efficacy of the PINN framework in simulating non-Newtonian fluid flow dynamics: a rectangular domain and a semi-circular domain.

Rectangular Domain:
Numerical Simulation of flow in a rectangle using PINNs is presented. The computational domain is a rectangle in two dimensions: a length of 1.1 cm and a height of 0.41 cm. The viscosity ν is defined as 0.01 gs^-1cm^-1. The total simulation time, T, is 0.5 seconds with a time step Δ t of 0.01 seconds.


The initial velocity field at t = 0 seconds is given by:

    𝐮(x, y, 0) = 
    [ u(x, y, 0); v(x, y, 0) ]
    =
    [ 0; 0 ]

for all points (x, y) within the domain.


The boundary conditions are given by:

    
  * Inlet (x = 0):
     A parabolic profile gives the velocity vector at the inlet boundary:
    
    𝐮(0,y) = 
        [ u(0, y); v(0, y) ]
        =
       0.9[ 4U_maxy(H - y)H^2(sin(π tT + 3π2) + 1);                                      0 ]


    where H is the height of the rectangle and U_max=0.5.

    
  * Outlet (x = L):
    At the outlet, a pressure condition:
    
    p(L, y) = 0


    Where L is the length of the rectangle.
   
    
  * Solid Walls (y = 0 and y = H):
    No-slip conditions at the solid walls:
    
    u(x, 0)    = 0,    v(x, 0)    = 0, 
    u(x, H)    = 0,    v(x, H)    = 0,




Our implementation of the weighted PINNs is based on using the TensorFlow framework. The neural network architecture consists of an input layer with three neurons, an output layer with five neurons, and seven hidden layers containing 50 neurons. The tanh function is utilized as the activation function for each hidden layer. A total of 3321 collocation points are generated for network training, which includes 244 boundary points (N_b) and 81 points on the inlet/outlet boundary (N_in/out), using Latin hypercube sampling (LHS). The balancing coefficient β is taken as 1, 2, 5, and 10.

For optimization, the network employs both Adam and L-BFGS optimizers. The network undergoes an initial training phase with the Adam optimizer for 5000 iterations, followed by a subsequent stage with the L-BFGS optimizer until convergence. The Hager-Zhang line search algorithm is utilized with a maximum of 50 iterations to determine the optimal step length α_k.

For prediction, 64561 collocation points are used, which includes 1124 boundary points (N_b) and 161 points on the inlet/outlet boundary (N_in/out).

For the WXPINN model, the computational domain was divided into M= 2, 3, and 4 non-overlapping subdomains. Each subdomain was modeled using a separate sub-PINN, whose network architecture mirrored that of the global model, and the interface conditions were enforced as described in the previous sections. The same number of collocation points as the global model were taken and equally distributed across the subdomains. The balancing coefficient γ is taken as 1, 2, 5, and 10 while keeping β fixed at 1. 

For the WCPINN model, the computational domain was divided into M= 2, 3, and 4 non-overlapping subdomains. Each subdomain was modeled using a separate sub-PINN, whose network architecture mirrored that of the global model, and the interface conditions were enforced as described in the previous sections. The same number of collocation points as the global model were taken and equally distributed across the subdomains. The balancing coefficient δ is taken as 1, 2, 5, and 10 while keeping β fixed at 1 and γ taken as 1 and 5. 

Semi-Circular Domain: 
Numerical Simulation of flow in a semi-circular domain using PINNs is presented. The computational domain is a semi-circular pipe with smooth disturbances in two dimensions, possessing a cross-sectional radius a = 1.6 cm and curvature radius R = 2.9 cm. The kinematic viscosity ν is set to 0.4 gs^-1 cm^-1. The total simulation time, T, is 6 seconds with a time step Δ t of 0.01 seconds.


The initial velocity field at t = 0 seconds is given by:

    𝐮(x, y, 0) = 
    [ u(x, y, 0); v(x, y, 0) ]
    =
    [ 0; 0 ]

for all points (x, y) within the domain.


The boundary conditions are given by:

    
  * Inlet (at the straight segment):
    A parabolic profile gives the velocity vector at the inlet boundary:
    
    𝐮_inlet = 
        [ u_inlet; v_inlet ]
        =
        0.9[ 4U_maxx(D - x)D^2(sin(π tT + 3π2) + 1);                                      0 ]

    where U_max = 0.75 is the maximum velocity at the inlet, and D is the cross-sectional diameter of the inlet segment.

    
  * Outlet (at the straight segment):
    At the outlet, a pressure condition:
    
    p(x_out, y) = 0

    where x_out denotes the position at the outlet.

    
  * Curved Boundary (semi-circular edge):
    No-slip conditions are applied at the curved boundary:
    
    𝐮_curved = 
        [ u_curved; v_curved ]
        =
        [ 0; 0 ]



The weighted PINN model is implemented in TensorFlow. The neural network architecture consists of an input layer with three neurons, an output layer with five neurons, and seven hidden layers containing 50 neurons. The tanh function is utilized as the activation function for each hidden layer. In total, 29760 collocation points are employed, with a batch size of 20000 points used in each training iteration. The balancing coefficient β is set to 1, 2, 5, and 10.

For optimization, the network employs both Adam and L-BFGS optimizers. The network undergoes an initial training phase with the Adam optimizer for the first few iterations, followed by a subsequent stage with the L-BFGS optimizer until convergence. The Hager-Zhang line search algorithm is utilized with a maximum of 50 iterations to determine the optimal step length α_k. For prediction, the exact total of 29760 collocation points is utilized. 
The numerical results showcased in this study were obtained using a Linux cluster located at IISER Thiruvananthapuram. The cluster comprises 88 nodes, with each node housing 28 compute cores running at 2.60 GHz and equipped with 128GB of RAM.



§ RESULTS



The weighted PINN model's capability to simulate fluid dynamics within a rectangular domain was evaluated through the visualization of velocity and pressure fields at various timestamps. These visualizations illustrated the dynamic flow patterns, identifying regions of varying velocity and pressure essential for understanding the fluid's behavior and the impact of forces on domain boundaries.


The model effectively captured the evolution of the flow field over time, highlighting significant variations in velocity across the domain, which is depicted in Figure <ref> for the rectangular domain. The velocity profile is consistent, with the maximum velocity observed in the middle of the inlet decreasing towards the walls due to the no-slip boundary conditions. The velocity in the middle of the channel is observed to continuously increase in the simulation time period due to the parabolic inlet profile.





Similarly, pressure distribution plots underscored the changes in pressure gradients driving the flow, pinpointing areas of high and low pressure crucial for fluid dynamics analysis, which is depicted in Figure <ref> for the rectangular domain. A pressure gradient decreasing from the inlet towards the outlet is observed, reflecting the zero-pressure boundary condition at the outlet, which drives the flow through the domain. The pressure plots exhibit pulsatile behavior, correlating with the sinusoidal inlet velocity. A cyclic increase and decrease in pressure is observed, mimicking the systolic and diastolic phases of blood flow. Comparable findings were obtained in <cit.> through conventional numerical techniques.




The weighted PINN model performance metrics are presented below in the case of the rectangular domain. The study explored the influence of different β values (1, 2, 5, 10) on model accuracy and computational efficiency. The metrics of interest were final loss (accuracy indicator) and computation time, alongside the number of iterations required for convergence, which are shown in Table <ref>.


 
β=1 yielded the most accuracy while β=2 took the least computation time after 26,438 iterations and 24,081 iterations, respectively, optimizing the balance between accuracy and computational effort. 
 
The progression of the loss function highlighted the model’s convergence behavior across different β settings is depicted in Figure <ref>, offering insights into performance optimization.



The weighted PINN model was also applied to a semi-circular domain to assess its performance in capturing fluid dynamics within more complex geometries.


In CFD simulations, it's a common challenge that the presence of incoming flow at open boundaries (backflow) can lead to unphysical oscillations and instability problems, even at moderate Reynolds numbers, see <cit.>. This issue arises from the convective energy entering the domain through the open boundary, which becomes problematic when the boundary velocity is unknown. Various stabilized finite element formulations have been proposed to address this problem in solving the incompressible Navier–Stokes equations. These formulations introduce stabilization terms based on the residual of a weak Stokes problem normal to the open boundary, driven by an approximate boundary pressure gradient.
In contrast, the framework of PINNs doesn't necessitate such stabilization terms to mitigate backflow instabilities. PINNs naturally handle such challenges without additional stabilization. In this study, we utilize a benchmark example problem to showcase the effectiveness of the PINNs approach. The PINNs model effectively captured the evolution of the flow field over time, highlighting significant variations in velocity across the domain, as shown in Figure <ref>. A disturbance in the flow profile is observed near the stenotic region of the artery as the velocity profile has a higher magnitude close to the origin of stenosis.





Similarly, pressure distribution plots underscored the changes in pressure gradients driving the flow, pinpointing areas of high and low pressure crucial for fluid dynamics analysis, which is presented in Figure <ref>. The pressure distribution also shows significant gradients around the stenotic region of the artery, with higher pressure upstream and lower downstream of the constriction. This highlights the impact of the narrowing on flow dynamics and the forces exerted on arterial walls.
Here, it's evident that we don't encounter any backflow instabilities, which stands out as the major advantage of using PINNs for addressing complex problems like this.




The weighted PINN model performance metrics are presented below in the case of the semi-circular domain. 
The study explored the influence of different β values (1, 2, 5, 10) on model accuracy and computational efficiency in Table <ref>. 


A β of 10 achieved the lowest final loss after 17,389 iterations, indicating high model accuracy. Conversely, β of 5 was the most efficient in terms of computation time and convergence speed, with 13,976 total iterations. The progression of the loss function highlighted the model’s convergence behavior across different β settings, offering insights into performance optimization.



The performance of the weighted Extended Physics-Informed Neural Network (WXPINN) model was assessed by varying the balancing coefficient γ and the number of subdomains M while testing with β=1. The choice of β=1 was taken, focusing on the model’s accuracy and computational efficiency, as that gave the best convergence in the weighted PINN case.

The study considered γ values of 1, 2, 5, and 10, alongside subdomain counts of 2, 3, and 4. Metrics of interest included final loss, computation time, and iterations required for convergence, providing a comprehensive view of the model's efficiency and accuracy.

Below, we demonstrate the WXPINNs efficiency w.r.t the number of subdomains:. 


    
  * Case I (M=2): In this scenario, γ=10 recorded the shortest computation time, while γ=5 achieved the lowest final loss, indicating the highest accuracy. See Table <ref> for more information.
    
  * Case II (M=3): In this scenario, γ=1 required the least computation time, and γ=5 presented the lowest final loss, showcasing optimal accuracy. See Table <ref> for more information.
    
  * Case III (M=4): In this scenario, γ=1 was most efficient in terms of computation time, whereas γ=5 maintained the lowest final loss, demonstrating superior accuracy. See Table <ref> for more information.




The evolution of the loss function across different configurations of the WXPINN model, varying by the number of subdomains (M) and balancing coefficient (γ), is visualized below. 



Consequently, a comparative analysis of computation times and loss evolution across various subdomain counts (M=2, 3, 4) elucidated the impact of domain partitioning on model efficiency, facilitating an understanding of accuracy-computational efficiency trade-offs.



The performance of the weighted Conservative Physics-Informed Neural Network (WCPINN) model was evaluated by varying the balancing coefficient δ and the number of subdomains M while maintaining a constant β=1 with γ varied as 1 and 5. The choice of γ=5 was taken, focusing on the model’s accuracy and computational efficiency, as that gave the best convergence in the WXPINN case.

The study considered δ values of 1, 2, 5, and 10, alongside subdomain counts of 2, 3, and 4. Metrics of interest included final loss, computation time, and iterations required for convergence, providing a comprehensive view of the model's efficiency and accuracy.

Below, we demonstrate the WCPINNs efficiency w.r.t the number of subdomains:


    
  * Case I (M=2): In this scenario, for γ=1, δ=10
     recorded the shortest computation time, while δ=2 achieved the lowest final loss, and for γ=5, δ=5 recorded the shortest computation time, while δ=2 achieved the lowest final loss, indicating the highest accuracy. See Table <ref> for more information.
    
  * Case II (M=3): In this scenario, for γ=1, δ=5 required the least computation time, and δ=2 presented the lowest final loss, and for γ=5, δ=1 required the least computation time, and δ=2 presented the lowest final loss, showcasing optimal accuracy. See Table <ref> for more information.
    
  * Case III (M=4): In this scenario, for γ=1, δ=10 was most efficient in terms of computation time, whereas δ=2 maintained the lowest final loss, and for γ=5, δ=5 was most efficient in terms of computation time, whereas δ=2 maintained the lowest final loss, demonstrating superior accuracy. See Table <ref> for more information.




The evolution of the loss function across different configurations of the WCPINN model, varying by the number of subdomains (M) and balancing coefficient (δ), is visualized below. 





A comparative evaluation of computation times further elucidates the impact of domain partitioning on computational efficiency, assisting in navigating the trade-offs between simulation accuracy and computational demand.




§ CONCLUSION

This study rigorously explores the capabilities of Physics-Informed Neural Networks (PINNs), including their advanced variants as weighted Extended PINN (WXPINN) and weighted Conservative PINN (WCPINN), to solve the blood flow model equations within cardiovascular systems. The study demonstrates that strategic subdomain partitioning and the optimal selection of balancing coefficients (β, γ, δ) are essential for enhancing computational efficiency while preserving high model accuracy. The adaptability of these models to diverse geometries underscores their potential for broader applications in simulating real-world fluid dynamics scenarios. Our PINNs simulation results clearly demonstrate the absence of backflow instabilities, highlighting a significant advantage of utilizing PINNs over traditional numerical methods for tackling complex problems like this.

However, the performance of PINN-based models heavily relies on hyperparameter configuration, necessitating extensive tuning for optimal settings. While domain decomposition improves generalization across subdomains, it introduces a tradeoff: simplifying each segment reduces required learning complexity but also decreases available training data per subdomain, potentially increasing overfitting risk and diminishing model generalizability. To address these challenges, we investigated WXPINN and WCPINN models incorporating weighted methodologies to optimize the distribution of training data and computational resources across subdomains for solving Navier-Stokes equations on complex domains. Based on presented numerical results and strategic adaptations, WXPINN and WCPINN provide a general approach to managing inherent tradeoffs in domain-decomposed neural network training, thereby enhancing the accuracy and applicability of PINN simulations in complex fluid dynamics.

Future endeavors should prioritize improving the generalization abilities of weighted PINNs across parallel computing frameworks and GPU accelerators. This enhancement aims to enhance computational efficiency, facilitating the scaling of these models on complex geometries and achieving real-time performance.



§ ACKNOWLEDGEMENT
 
We gratefully acknowledge the financial support under R&D projects leading to HPC Applications (DST/NSM/R&D_HPC_Applications/Extension/2023/33), National Supercomputing Mission, India, and we greatly acknowledge the support for high-performance computing facilities at the Padmanabha cluster, IISER Thiruvananthapuram, India. Additionally, we extend our gratitude to Miss Soundarya Sarathi for her initial code implementation.

ieeetr


