PaperlayerView exact revision on arXiv

Source-linked AI summary

A Bifurcation-Based Domain Decomposition Method with Neural Operators for Blood Flow Simulation

Yuzhou Zhao, Han Zhang, J. Matias Di Martino, Jean-Michel Morel, Guillermo Sapiro

arXiv:2610.02238v2cs.CEphysics.comp-ph

TL;DR

Complex vascular simulation is costly, and existing learned approaches can require extensive data and topology-specific adaptation. This paper combines bifurcation-unit operators, Windkessel parameter lumping, and waveform relaxation; it reports 13×–17× speedups with low errors and topology generalization.

  • Problem

    Existing supervised approaches require expensive high-fidelity training data and may need retraining or substantial adaptation when vascular geometry or topology changes.

    Source: §1 · p. 1

  • Method

    The method decomposes vascular trees into bifurcation units, lumps downstream Windkessel parameters, predicts local fields with an operator, and couples units through Schwarz waveform relaxation.

    Source: §4.2 · pp. 9, 11, §7 · p. 25

  • Results

    13× and 17× speedups over Nektar1D on the 7- and 55-segment networks, respectively, with 97.8–99.7% agreement and PWV biomarkers within 1–2% of the reference.

    Source: §6.8 · p. 25

  • Takeaways & Limitations

    The same checkpoint generalizes to pruned arterial trees without retraining, while full-network solutions preserve PWV biomarkers within 1–2% of the Nektar1D reference.

    Source: §6.8 · p. 25, §7 · p. 25

    A rigorous convergence proof for multiplicative SWR on coupled vascular networks with bifurcations remains open.

    Source: §7 · p. 26

Abstract

from arXiv · show

Fast and accurate simulation of hemodynamic behavior within vascular networks is essential for numerous clinical applications. However, obtaining high-quality and computationally efficient flow measurements across complex vascular networks remains challenging. To address this, we first decompose the vascular network into a set of bifurcation units and then develop an operator network capable of mapping unit-specific parameters to the local solution fields of each bifurcation unit. By lumping the Windkessel-model outlet parameters and incorporating inlet boundary conditions from the solution of parent units, the flow and pressure fields can be rapidly approximated. Subsequently, operator-network-driven Schwarz waveform relaxation is applied across bifurcation units to correct discontinuities and improve numerical accuracy. On 7-segment and 55-segment arterial tree models, the proposed method achieves $13\times$ to $17\times$ wall-clock speedups over conventional 1D numerical simulation, with relative $L^2$ errors of 1% in both pressure and velocity. The resulting pulse wave velocity biomarkers agree with the conventional reference to within 1--2%, and the same trained operator generalizes to different tree-like 1D vascular network topologies.

1. Introduction.

The paper addresses costly, topology-sensitive vascular simulation with reusable bifurcation-unit operators and operator-driven Schwarz waveform relaxation. It reports faster simulations with low errors and topology generalization.

1. Introduction.

  • Existing supervised models require costly high-fidelity data and often need retraining or substantial adaptation when vascular geometry or topology changes.Accurate simulation data for one complex network may take hours to days to generate.Source: §1 · p. 1
  • The method decomposes networks into reusable bifurcation units and learns operators mapping unit-specific parameters to local solution fields.Predicted unit solutions are assembled into a global approximation, avoiding a separate global model for each network configuration.Source: §1 · pp. 1–2
  • Lumped Windkessel outlet parameters and parent-unit inlet waveforms provide unit boundary conditions; operator-driven Schwarz waveform relaxation then corrects interface discontinuities.The operator network approximates flow and pressure fields for each unit.Source: §1 · p. 2
  • On a 7-segment network and the 55-segment full arterial tree, operator-driven SWR achieved 13× to 17× wall-clock speedups over conventional 1D simulation, with mean relative L2 errors of the order of 1% in pressure and velocity.Downstream biomarkers were consistent with numerical reference solutions, and the method generalized to incomplete vascular topologies.Source: §1 · pp. 2–3

2. Related Work.

Prior work includes reduced-order blood-flow models, operator learning, and Schwarz waveform relaxation, but the paper targets their composition for topology-varying vascular networks. It uses learned bifurcation-unit solvers within waveform relaxation.

2.1. Blood Simulation.

  • Existing blood-flow approaches include 1D models and nodal predictors, while the paper targets full local pressure and velocity fields with topology-generalizable bifurcation units.The related-work discussion identifies limited scalability to unseen topologies and nodal-only outputs as gaps in existing methods.Source: §2.1 · pp. 3–4

2.2. Operator Learning.

  • Prior operator-learning methods address fixed or single domains; this work maps varying bifurcation-unit parameters to local solution fields and composes units through Schwarz iterations.The stated gap includes varying geometry, lumped boundary conditions, and network topology.Source: §2.2 · p. 4

2.3. Domain Decomposition Methods and Schwarz Waveform Relaxation.

  • Schwarz waveform relaxation exchanges time-dependent interface waveforms, and the paper applies it to vascular 1D Navier-Stokes equations using surrogate rather than conventional subdomain solvers.The related-work passage states that prior cited SWR applications had not addressed these vascular equations.Source: §2.3 · pp. 4–5

3. Problem and Preliminaries.

The paper models tree-like vascular networks with 1D Navier-Stokes equations and bifurcation junctions. It aggregates downstream Windkessel conditions to support unit-wise decomposition.

3. Problem and Preliminaries.

  • The network is represented as directed vessel segments connected at nodes, with binary bifurcations used as the common unit; higher-branch junctions can be reformulated as sequential bifurcations.The stated scope is tree-like 1D topologies.Source: §3 · p. 5

3.1. Blood Flow via 1D Navier-Stokes.

  • Under impermeable-wall and incompressible Newtonian-fluid assumptions, the 1D Navier-Stokes model uses area, flow rate, and pressure as its primary unknowns.The model includes physical, geometric, and boundary-condition parameters; inlet flow is prescribed and the initial state is quiescent.Source: §3.1 · pp. 5–6

3.2. Parameter Lumping.

  • Child-vessel Windkessel models are reduced and combined into an equivalent parent-segment model, recursively aggregating downstream resistance and compliance.The equivalent proximal resistance is chosen as characteristic impedance to minimize wave reflections.Source: §3.2 · pp. 7–8

4. Bifurcation-Based Domain Decomposition.

Bifurcation units make vascular-network decomposition tractable by inheriting child inlets from parent solutions and representing downstream effects with lumped Windkessel outlets. A top-down initialization is then refined through iterative waveform exchange across units.

4.1. Rationale for Using Bifurcation Units.

  • Bifurcation units avoid the unavailable flow-splitting relations required by segment-wise decomposition, while letting child inlets inherit states from parent solutions and downstream effects use lumped Windkessel models.The entire network otherwise must be solved simultaneously to satisfy junction and boundary conditions.Source: §4.1 · p. 8, §4 · p. 8

4.2. Multi-Level Schwarz Waveform Relaxation by BU with Parameter Lumping.

  • The three-phase workflow aggregates Windkessel outlet parameters bottom-up, initializes BU solutions top-down, then applies Schwarz waveform relaxation to update interfaces.The root uses a prescribed-flow inlet; downstream BUs inherit full-state inlet conditions from parent solutions, while terminal BUs retain Windkessel outlets.Source: Figure 4 · p. 9, §4.2 · pp. 10–11
  • SWR exchanges full Q,P waveforms across BU interfaces in sequential top-down iterations; transferring only one variable generally led to instability and divergence.Terminal BUs retain their prescribed Windkessel outlet models during the iterations.Source: §4.2 · pp. 9, 11, 12
  • Iterations reduce interface mismatches until a convergence criterion is met, and trained operators replace expensive high-fidelity numerical solves for each BU.The authors report monotonically decreasing residuals empirically and use a fixed iteration budget in practice.Source: §4.2 · pp. 11–12

5. Operator Network for Bifurcation Unit.

The paper introduces a local BU operator that maps unit parameters and boundary conditions to solution fields without using global network topology. BUFormer encodes these inputs into tokens, decodes temporal waveforms, and is trained against reference fields using relative L2 error.

5. Operator Network for Bifurcation Unit.

  • BUFormer serves as a reusable local solver: one forward pass predicts converged fields without global-topology inputs, enabling reuse across vascular trees decomposable into BUs.The same trained operator is described as applicable without retraining or architectural modification.Source: §5 · p. 12, §5.1 · p. 13

5.1. BU Input-Output Map.

  • The operator maps BU geometry and physical parameters plus inlet and outlet conditions to solution fields on its three segments, covering six possible boundary-condition configurations.Inputs include root inflow or inherited full-state traces, and Windkessel parameters or downstream Dirichlet traces.Source: §5.1 · pp. 12–13, Table 2 · p. 13

5.2. Network Architecture.

  • BUFormer encodes scalar parameters and boundary waveforms into four tokens; spatial query points cross-attend through Transformer blocks, and a Fourier Decoder produces temporal solution waveforms.Mask-based switches route the different inlet and outlet boundary-condition types through the appropriate encoding branch.Source: §5.2 · pp. 13–14, Figure 5 · p. 14

5.3. Optimization Objective.

  • Training minimizes relative L2 error between predicted and reference solution fields concatenated across the three segments, spatial points, fields, and time steps.The same relative L2 metric is used to evaluate accuracy in the experiments.Source: §5.3 · p. 14

6. Experimental Results.

Across the tested arterial networks, operator-driven Schwarz waveform relaxation (SWR) speeds simulation while maintaining agreement with the numerical reference. The results also assess local operator accuracy, pulse wave velocity (PWV) agreement, and robustness to changed topologies and physiological conditions.

6.1. Experimental Setup.

  • The evaluation used 7- and 55-segment networks, virtual-subject datasets, and Nektar1D as the ground-truth reference; the 55-segment dataset required about 35 HPC-cluster hours to generate and 20 hours to train.The solver comparisons included numerical SWR and operator SWR configurations.Source: §6.1 · pp. 15–17

6.2. Warm-Start Ablation.

  • After 20 iterations, warm-start SWR reached residual, pressure-error, and velocity-error values of 2.66 × 10−4, 1.32 × 10−4, and 3.07 × 10−4, respectively, versus 1.05 × 10−2, 7.48 × 10−2, and 4.32 × 10−2 with zero initialization.The warm start began with a residual roughly one order of magnitude lower and maintained its advantage.Source: §6.2 · p. 17, Figure 7 · p. 18

6.3. BU-Level Validation.

  • Median BU relative L2 errors were 0.52% on the 7-segment model and 0.20% on the 55-segment model, with 95th-percentile errors of 0.74% and 0.41%, respectively.The operator also maintained low errors across boundary-condition configurations.Source: §6.3 · pp. 19–20

6.4. Network-Level Performance.

  • 13.2× speedup over Nektar1D with 99.68% pressure and 99.62% velocity agreement on 7 segments; 17.3× speedup with above 97.7% agreement in both fields on 55 segments.Timings are per-case comparisons against Nektar1D’s single-process CPU benchmark.Source: §6.4 · p. 20, Table 6 · p. 20

6.5. Clinical Relevance Through PWV Robustness.

  • PWV errors ranged from 0.7% for cfPWV and baPWV to 1.7% for faPWV; aPWV error was 1.1%.The operator and Nektar1D population-level PWV distributions were nearly indistinguishable.Source: Table 8 · p. 22, §6.5 · p. 23

6.6. Generalization to Different Topologies.

  • At iteration 30, mean relative L2 errors increased from 1.15% pressure and 2.21% velocity on the full 55-segment tree to 3.04% and 8.68% on the pruned 35-segment network.The same full-tree-trained checkpoint was applied to pruned networks without retraining.Source: Figure 11 · p. 23, Table 9 · p. 23, §6.6 · p. 24

6.7. OOD with Elastic Vascular Aging and Hypertension.

  • In tested out-of-distribution aging and hypertension cases, mean pressure/velocity errors rose from 6.2%/8.9% for the milder 40-year phenotype to 9.0%/15.0% for the stiffer 50-year phenotype; all cases remained stable.Errors could accumulate over Schwarz iterations and along the vascular tree, making the most extrapolated conditions and distal segments more challenging.Source: §6.7 · p. 24, Table 10 · p. 24

6.8. Discussion.

  • Domain decomposition makes the problem amenable to surrogate modeling: local BU operator evaluations can run in parallel at each SWR iteration.The framework is not tied to the transformer backbone; the paper also reports testing a Fourier Neural Operator.Source: §6.8 · p. 25

7. Conclusion.

The framework combines bifurcation-unit operators with Windkessel parameter lumping and Schwarz waveform relaxation to assemble globally consistent vascular solutions. Experiments report accurate, faster simulations and topology generalization, while validation and theoretical analysis remain bounded by stated assumptions and open questions.

7. Conclusion.

  • The three-phase method aggregates downstream parameters, warm-starts local predictions, then applies operator-driven multiplicative Schwarz waveform relaxation to assemble a globally consistent solution.The BUFormer replaces the numerical solver at each bifurcation unit.Source: §7 · p. 25
  • Mean relative L2 errors are of order 1% in pressure and velocity, and pulse wave velocity biomarkers agree with the Nektar1D reference within 1–2%.The BU operator's median error is below 0.5%; the conclusion also reports generalization to different tree-like topologies, including incomplete networks, without retraining.Source: §7 · p. 25
  • Validation uses healthy-subject physiological assumptions; extending the training data to pathological regimes such as severe stenoses and aneurysms is identified as future work.The authors also note possible extensions to trifurcations and side-branch junctions.Source: §7 · p. 25
  • A rigorous convergence analysis remains open for the coupled bifurcation network under full-state (Q, P) coupling.The authors identify establishing well-posedness and convergence for the coupled network as a key step.Source: §7 · p. 26

Appendix A. Hyperparameters.

  • The 55-segment BUFormer uses a larger embedding dimension and more Transformer layers than the 7-segment setting to accommodate greater configuration diversity.Experiments were conducted on a single NVIDIA A100 GPU.Source: Appendix A · p. 30, Table 11 · p. 30

Appendix B. Baseline Arterial Parameters.

  • Baseline simulations assume purely elastic arterial segments and standard physical values; the authors state that minor deviations are unnecessary for this study's scope.The framework can handle different values in principle.Source: Appendix B · p. 30

Appendix C. Training Data Normalization.

Training inputs and targets are normalized using statistics fixed from the training split, with predictions returned to physical units for coupling and error evaluation. The supplied passages also report SWR convergence across tested factors and identify the source of rare extreme-range breakdowns.

Appendix C. Training Data Normalization.

  • All operator inputs and targets are normalized, with training-split statistics frozen for validation, testing, and inference; predictions are denormalized for SWR coupling and error reporting.The text notes that fields span several orders of magnitude and have different units.Source: Appendix C · p. 30
  • Waveforms are rescaled to normalized time, and pressure and velocity use global training-set means and standard deviations; inlet flow uses its own statistics.Relative L2 accuracy is evaluated per field after denormalization.Source: Appendix C · pp. 30–31
  • Scalar operator inputs are normalized using training-set statistics, with quantities spanning many orders of magnitude log-transformed first.Inputs include geometry, wall stiffness, reference pressure, cardiac period, and lumped Windkessel boundary conditions.Source: Appendix C · p. 31
  • Multiplicative SWR converges across the tested network-size, impedance/compliance, stiffness, heart-rate, and inflow-waveform ranges; rare non-physiological-extreme breakdowns stem from the numerical solver.These experiments use Nektar1D at each bifurcation unit, so the reported convergence concerns the SWR coupling itself.Source: Appendix C · p. 32

Appendix D. SWR Convergence Analysis.

The figures depict SWR convergence through residuals and pressure and velocity errors across network sizes and, separately, across iterations for selected segments of the 55-segment model.

Appendix D. SWR Convergence Analysis.

  • Fig. 12 plots residual and mean relative L2 pressure and velocity errors against iteration for arterial networks ranging from 11 to 55 segments.All networks use warm-start initialization.Source: Figure 12 · p. 32
  • Fig. 13 plots per-segment relative L2 pressure and velocity errors against Schwarz iteration for 12 anatomical segments in the 55-segment model.Curves are colored by cumulative centerline distance from the aortic root to each segment midpoint.Source: Figure 13 · p. 32

D.1. Convergence vs. network size.

SWR converges monotonically on the tested anatomical networks, with convergence slowing as network size increases. Segment-level errors converge fastest near the root, while distal segments lag as boundary information propagates through the top-down sweep.

D.1. Convergence vs. network size.

  • All six networks converge monotonically in the mean, with smaller networks converging faster and the gap across three metrics remaining roughly one order of magnitude.Fig. 12 shows the residual and mean relative L2 errors over 50 warm-started Schwarz iterations for networks from 11 to 55 segments.Source: §D.1 · pp. 32–33
  • By iteration 50, the 11-segment network reaches a residual of 2.81 × 10−4 and mean pressure and velocity errors of 1.12×10−4 and 3.08×10−4, respectively.Source: §D.1 · p. 33
  • Proximal segments near the aortic root converge fastest, while distal segments lag as information propagates outward through successive top-down Schwarz iterations.The reported distal examples include the left femoral and left anterior tibial segments.Source: §D.1 · p. 33
  • Fig. 14 compares self-convergence residual and mean pressure and velocity errors over iterations for balanced and strongly unbalanced synthetic 15-segment trees.In the unbalanced tree, one first-generation subtree has segments four times longer than the other.Source: Figure 14 · p. 33

D.2. Network topology.

The convergence sweeps vary tree topology, terminal loads, vessel stiffness, heart rate, and inlet waveform on a synthetic vascular tree. SWR converges across the tested settings, with accuracy affected by some parameter changes.

D.2. Network topology.

  • On the unbalanced tree, SWR converges to a residual of ≈4 × 10−6, while mean relative L2 error settles at ≈1.4 × 10−4, above the balanced tree’s ≈7 × 10−6.The unbalanced tree has one first-generation subtree with segments four times longer than the other.Source: §D.2 · pp. 33–34

D.3. Terminal resistance and compliance.

  • Compliance scaling leaves mean relative L2 error near 5 × 10−6, while high resistance (RT ×2) reaches 1.3 × 10−4 by iteration 20 with error still dropping.At the two resistance extremes, the local numerical solver fails before Schwarz coupling starts.Source: §D.3 · p. 34, Figure 15 · p. 34

D.4. Vessel stiffness.

  • Across tested wall stiffnesses, SWR converges; the soft wall (β×0.5) reaches ≈1.8×10−4, while at β×0.25 the monolithic reference solver crashes before comparison.The solver crash at β×0.25 is attributed to the reference solver, not the SWR coupling.Source: §D.4 · p. 34

D.5. Heart rate.

  • Across cardiac periods from 1.5 to 0.5 s (40 to 120 bpm), multiplicative SWR converges at every rate, with mean relative L2 error of ≈3–12 × 10−6.The half-sine inlet period tracks the cardiac period in this sweep.Source: §D.5 · pp. 34–35, Figure 17 · p. 35

D.6. Inflow waveform.

  • Four structurally distinct inlet waveforms with matched cycle-mean flow all converge, with residuals ≈2–3 × 10−6 and pressure errors ≈1.7 × 10−5 to ≈6.5 × 10−5.The tested shapes are physiological, half-sine, narrow Gaussian, and triangular pulses.Source: §D.6 · p. 35

Appendix E. Parallelized Throughput and Hardware-Usage Disclosure.

The throughput analysis measures the operator on one A100 GPU and the reference solver on a CPU core, and evaluates batched inference on the 55-segment tree.

Appendix E. Parallelized Throughput and Hardware-Usage Disclosure.

  • The reported speedups compare GPU operator deployment with a CPU-core numerical solver, rather than measuring per-FLOP efficiency.Reference simulations run one per CPU core; multi-core throughput scales with core count for independent cases.Source: Appendix E · p. 35
  • At batch size 128, bf16 inference takes 0.54 s per case, a 3.8× speedup over fp32 with only a ∼0.2 percentage-point increase in pressure and velocity relative L2 errors.In fp32, per-case time falls from 6.30 s at batch size 1 to 2.05 s at batch size 128.Source: Appendix E · p. 36
  • At this throughput, the batched operator matches Nektar1D running on roughly 210 CPU cores, and the cited HPC node is 3–4× more expensive than one A100.The cost comparison is stated at equivalent throughput.Source: Appendix E · p. 36

Appendix F. Ablation Studies on Backbone.

The backbone ablation tests whether the bifurcation-unit decomposition and Schwarz coupling can use an operator architecture other than BUFormer.

Appendix F. Ablation Studies on Backbone.

  • Because a plain FNO cannot handle BUFormer’s BC-dependent boundary layouts, the study trains a separate, parameter-matched FNO for each boundary-condition type.BUFormer remains the default because its shared mask-based encoder handles all boundary-condition combinations in one model.Source: Appendix F · p. 37
  • The decomposition and Schwarz coupling also function with an FNO-per-type backbone on both trees, supporting that the framework is not tied to BUFormer.Table 14 reports full-network SWR accuracy on the 7- and 55-segment trees, with BUFormer as reference.Source: Appendix F · p. 37, Table 14 · p. 37

Appendix G. Ablation Studies on Fourier Mode and Balanced Sampling.

The ablations vary retained Fourier modes and balanced sampling over boundary-condition types in the 7-segment setting, changing one factor at a time.

G.1. Fourier modes.

  • At K = 128, velocity error is 5.2 × 10−3, down from 4.7 × 10−2 at K = 16; it reverses at K = 256, while pressure error is smallest at K = 128.The reported pressure error stays below 1 × 10−2; K = 128 is used for the main model.Source: §G.1 · pp. 37–38, Figure 19 · p. 38

G.2. Balanced sampling over BC types.

  • With balanced sampling, every type reaches ≈2–3×10−3, versus ≈1.8 × 10−2 without balancing; the rare root type 1 is worst without balancing.Without balancing, all types stall after ∼60 epochs; with balancing, every type continues improving.Source: §G.2 · pp. 38–39, §G.1 · p. 38, Figure 20 · p. 38

Related papers

Loading 2610.02238v2…