Generate a physical vascular network¶
The vascularization library grows physical vessels inside a prepared 3D
background. It keeps vessel identity separate from numerical 1D cells.
NetworkData stores physical node, vessel and tissue-terminal identifiers,
radii, viscosity, system, orientation, ownership and provenance. A successful
GenerationResult contains the proposed network, its numerical subdivision,
common tissue territories and validation reports. The caller adopts it only
when validation.valid() is true; input data are never modified.
Generation from an empty arterial input is collective and supports one or several MPI processes. The background point owner creates the single real seed root; its branch may cross partitions immediately. Imported arterial networks also support collective startup from one or several selected leaves and collective terminal-count calibration. Arterio-venous generation starts from real venous outlets or selected input venous leaves and uses the same MPI path with one or several processes. Use a first-order mapping on an all-tetrahedral or all-hexahedral physical background. Geometry checks include the vessel wall and clearance. An ambiguous boundary separation at the finite geometry budget fails conservatively. Explicit planar inlet/outlet sections are supported; a warped boundary face is not replaced by a planar aperture.
The algorithmic references are the OpenCCO article and implementation and svVascularize. ImmersX implements the background-owned MPI growth and commit protocol described below using deal.II ownership, spatial search and sparse exchanges.
Length scale and physical data¶
Vascularization/Length scale defaults to local mesh size. The configured
lengths are coefficients of a physical reference \(H\):
Policy |
Physical reference |
|---|---|
|
One metre |
|
Diameter of the containing active background cell |
|
Global minimum active-cell diameter |
|
Global maximum active-cell diameter |
Only owned cells contribute to the global bounds. Empty ranks participate. On interfaces, the smallest containing diameter is used, including across MPI partitions. The prepared background and its reference scale remain fixed during generation and time integration.
The generator retains one background search context across startup, both systems and calibration attempts. It is created after background validation and destroyed when generation returns; it cannot survive mesh or mapping changes. Each calibration attempt starts with independent network and tissue state. Final coverage, subdivision and output remain separate measured phases.
Growth/Seed length defaults to one and gives a physical length
\(cH(x_{attachment})\), for both new seeds and imported attachments. The fixed
seed or attachment radius remains physical. A venous startup uses this seed
stem followed by a straight continuation to the assigned territory endpoint.
The junction may bend; an unbranched junction preserves its radius and carries
the same flow on both sides. Deterministic geometric retries can fail when the
requested stem or continuation does not fit. They do not shorten the seed or
change its radius.
Growth/Minimum vessel length defaults to 0.1 and requires
\(\int_\Gamma ds/H(x)\ge c\) for every new physical branch. The complete exact
cell-crossing profile is used, including interface portions. Optimization
prechecks use the global minimum reference as a necessary lower bound; selected
proposals undergo the exact complete-profile check before adoption. A midpoint or
endpoint does not replace that integral. Physiology/Minimum radius defaults
to zero; a positive value requires \(r(x)/H(x)\ge c\) throughout each generated
branch. The conservative check includes the largest crossed scale. Murray
radii are computed once from demand and fixed physical attachment radii.
Growth/Clearance defaults to zero. Let \(H_{max}(\Gamma)\) be the maximum
reference along a branch. Its wall clearance is \(cH_{max}(\Gamma)\). For two
branches, their additional separation is
\(c\max(H_{max}(\Gamma_a),H_{max}(\Gamma_b))\) beyond their radii. This rule is
symmetric and includes the scale variation along both complete branches.
An imported branch’s local maximum includes only its inside portions. A fully
exterior obstacle has no local scale; the generated branch supplies the defined
scale in their pair maximum. New exterior branches remain ineligible.
Routing boxes and ghost halos use the global maximum reference. Geometry
profiles are cached independently of radii and discarded when their geometry
leaves the current snapshot or trial set.
Venous network/Terminal separation defaults to 0.25 times \(H\) at the common
territory centre, shared by the pair. This separates endpoints by a quarter
cell diameter while leaving each physical radius and all supplied physiology
explicit; it is not an organ-specific calibration. The examples and geometric
gates exercise it with the new seed and minimum-length coefficients.
Objectives/Maximum service distance defaults to one. Coverage uses
\(cH(x)\) at each tissue quadrature point, with physical JxW weights. It remains
independent of perfusion and of nearest-territory demand assignment.
Domain size, coordinates, imported geometry and properties, fixed seed and attachment radii, viscosity, pressures, total flow and supplied bed resistances remain in SI. Generated coordinates, radii, archives and solver transfers also remain in SI; the configured coefficients are not overwritten. Local \(H\) is undefined outside the background: exterior input is preserved, new exterior attachments are rejected, and local normalized diagnostics explicitly report only inside portions. Constant policies are defined outside too.
Arterial and venous trees¶
This synthetic example uses metres, pascals, seconds and cubic metres per
second. Its values are not an organ calibration. The background is a unit cube.
The VascularAV integration tests check paired flow and Murray constraints in
Debug and Release.
#include <deal.II/distributed/tria.h>
#include <deal.II/distributed/fully_distributed_tria.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/grid_generator.h>
#include <immersx/vascularization/generation.h>
using namespace dealii;
using namespace ImmersX::Vascularization;
parallel::distributed::Triangulation<3> background(MPI_COMM_WORLD);
GridGenerator::hyper_cube(background);
background.refine_global(2);
MappingQ1<3> mapping;
Parameters parameters;
parameters.seed_position = Point<3>(.2, .5, .5);
parameters.seed_direction = Point<3>(1., 0., 0.);
parameters.seed_radius = .003;
parameters.seed_length = 1.;
parameters.viscosity = .004;
parameters.terminal_count = 1;
parameters.terminal_count_multiplier_factor = "background cell count";
parameters.reference_flow = 2.e-9;
parameters.tissue_demand = "1+x+2*y";
parameters.arterial_murray_exponent = 2.7;
parameters.venous.enabled = true;
parameters.venous.seed_position = Point<3>(.8, .5, .5);
parameters.venous.seed_radius = .004;
parameters.venous.murray_exponent = 3.2;
parameters.venous.terminal_separation = .25;
parameters.reference_conditions.pressures = {{0, 12000.}, {2, 4000.}};
auto result = generate(background, mapping, {}, parameters, MPI_COMM_WORLD);
AssertThrow(result.validation.valid(), ExcMessage("Vascular generation failed"));
parallel::fullydistributed::Triangulation<1, 3> vessels(MPI_COMM_WORLD);
const auto physical_cells =
create_triangulation(result.network, result.discretization, vessels);
For empty input, the arterial inlet is physical node 0. Arterial seed creation
also creates its first endpoint, node 1. When a new venous outlet is needed,
its physical node is reserved next, node 2, before further insertions. Its
identifier therefore stays fixed when terminal-count calibration changes the
number of territories. These pressure keys are physical node identifiers,
not triangulation vertex numbers or DoF indices. Imported networks retain their
own port identifiers.
The arterial forest defines common territory centres and their demand. The same CCO engine then grows veins from the declared outlet or selected venous attachments toward these territories. The first venous connection joins a fixed-length seed stem to an admissible territory, trying territories in distance/id order and deterministic stem directions. Each later insertion tries six Cartesian endpoint offsets of the configured separation. The territory centre is the arterial endpoint; its venous endpoint is distinct. Each territory receives its integrated tissue demand, shared exactly by both terminals. A vein’s geometry is oriented from outlet toward tissue, and its physical reference flow is negative in that orientation.
The arterial forest remains fixed while the venous forest grows. All venous proposals, including radius changes on their ancestors, are checked against both forests and imported obstacles. This ordered construction can exhaust its finite geometry budget even when another joint arrangement might exist. Such failure is reported; it does not establish geometric impossibility.
Each system has its own Murray exponent. The radius of a generated vessel
satisfies \(r_e=R_{\mathrm{root}}(D_e/D_{\mathrm{root}})^{1/\gamma}\), where
\(D_e\) is the sum of downstream territory demands. Root radii stay fixed. The
minimum-radius check includes every generated branch, including internal
branches crossing coarse background cells, after each joint demand update.
terminal_weight and density_beta control where territories are sampled;
tissue_demand controls supply and radii independently.
Integrated tissue demand¶
Physiology/Tissue demand supplies the relative density \(d(x)\). The tissue is
assigned to its nearest common territory centre; stable physical territory
identifiers break equal-distance ties. This assignment is independent of the
maximum service distance used for coverage. A territory’s demand is
The implementation approximates this integral with physical JxW weights on
owned background cells. Adding centres changes existing territories too: their
absolute demands, dependent Murray radii and growth costs are updated together.
The selected batch is integrated jointly, since its new demands can differ from
isolated candidate previews. Rejection restores tissue assignments and physical
state together. Arteries and veins use the same \(D_i\) for each common territory;
the explicit total flow is distributed as \(Q_i=Q_{total}D_i/\sum_jD_j\).
Refining the proposal density therefore does not multiply the total tissue
integral. Input terminal weights outside the generated forest remain unchanged.
The integration controls are in Physiology:
Parameter |
Default |
Meaning |
|---|---|---|
|
2 |
Gauss order on reference hexahedra or tetrahedra; supported orders are 1 through 8. |
|
64 |
Initial global quadrature-point target per requested territory, after the terminal multiplier. |
|
4000000 |
Maximum global number of active demand quadrature points. |
|
8 |
Additional integration subdivision rounds permitted in each territorial preview. |
Integration subdivisions leave the background mesh and its length scales
unchanged. Dyadic subdivisions can exceed the initial point target; the point
budget remains a hard limit. A territory with no quadrature mass triggers local
integration refinement. If positive demands cannot be resolved within the
budgets, the candidate is rejected explicitly; no demand floor is substituted.
integration_budget_rejections in each growth CSV records these rejections.
demand_quadrature_points is the accepted global point count and
demand_quadrature_peak_points is its global proposed peak, including rejected
previews. Use the maximum across ranks for these cardinalities. Assignment
queries, hits, generated/refined points and reassigned samples are local work
counts and can be summed across ranks. Index counters include persistent
network indices, small MPI owner-box indices and temporary candidate trees.
Read their insertion and hit totals alongside build counts to see the number
of records processed. The work columns include unsuccessful
attempts; the initialization_ and finalization_ prefixes separate their
one-time contributions from ordinary insertion work. Timings use the same
phase separation, with finalization outside the accumulated batch total.
Initialization work also includes startup and cache preparation outside the
forest-growth call. That work contributes to overall generation_seconds;
initialization_seconds times initialization inside forest growth.
Failure to reach the terminal target also reports any encountered integration
budget exhaustion. Invalid or nonrepresentable density/JxW values are input
errors rather than geometric retries.
The density and its weighted samples must be positive, finite and representable. Finite quadrature approximates territory boundaries: obtaining positive mass for every territory is not an accuracy estimate. Increase the quadrature order or samples per terminal and compare the resulting demands when individual territory accuracy matters.
Terminal sampling and numerical resolution¶
Arterial territory proposals use the relative physical density
Here \(|K(x)|\) is the physical volume of the containing active background cell,
computed by deal.II for the first-order mapping. On cell and MPI interfaces,
the smallest containing diameter selects the cell, with cell identity breaking
ties; the lowest containing process owns the point.
terminal_weight supplies the nonnegative dimensionless expression \(w(x)\),
bounded by terminal_weight_bound. The corresponding standalone parameters
are Sampling/Spatial weight and Sampling/Spatial weight bound.
density_beta (Sampling/Density mesh exponent) defaults to three. Zero makes
the sampling density independent of mesh size. Three gives
equal integrated proposal weight per cell when \(w=1\), including cells of
different shape or aspect ratio. Geometric constraints can change the
accepted terminal distribution.
The generator first chooses an owned cell with mass
\(|K|^{1-\beta/3}\), then draws uniform points in its physical bounding box.
The native mapping inverse and ReferenceCell::contains_point select points
inside that same cell. Mapped cell volumes use physical
JxW; reference-uniform points would bias distorted cells. Rank proposal
probabilities are proportional to their total cell mass. The remaining rejection
uses only \(w/\mathrm{terminal\_weight\_bound}\).
The immutable cell table is cached with the background, including across
calibration attempts. Its common logarithmic normalization leaves relative
weights unchanged. This avoids spending the insertion budget on global-box
draws rejected solely because a small refined region has tiny cells. The
sampler reports exhaustion of its 4096 inverse-map attempts
explicitly; it never switches to a different cell or substitutes a density.
Unrepresentable cell-mass ranges are rejected before growth.
sampled_points counts generated physical-cell points before spatial-weight and
geometry acceptance. Rank mass skips produce no point. Inverse-map retries
contribute to sampling time. There is no separate
\(\lambda\) parameter: terminal_count and its multiplier set the initial number
of new tissue terminals, and calibration changes that number to seek the volume
target.
Sampling/Terminal count multiplier factor selects how the initial terminal
count follows the prepared background:
Value |
Initial number of new territories |
|---|---|
|
|
|
|
|
|
The default requests one new territory per active cell when Terminal count
is one. Select none to request an absolute count.
Shared MPI vertices count once, including hanging nodes; ghost and artificial
cells do not increase the cell count. The count is recomputed on each call to
generate() after the caller has prepared the background. A zero or overflowing
product is rejected collectively. The supplied parameters are not modified.
In arterio-venous mode the count is the number of common territories, with one
arterial and one venous terminal per new territory, not their combined count.
For example, to request two new territories per active background cell:
subsection Vascularization
subsection Sampling
set Terminal count = 2
set Terminal count multiplier factor = background cell count
set Density mesh exponent = 3
end
subsection Objectives
set Target volume fraction = 0
set Minimum tissue coverage = 0
end
end
This scales the total; it does not assign an exact quota to each cell or vertex.
Density mesh exponent and Spatial weight still control relative proposal
locations. Uniform refinement increases the total when a background multiplier
is selected, even with exponent zero.
With calibration enabled, multiplication applies once to its initial guess.
Minimum new terminals and a nonzero Maximum new terminals are absolute
counts. A zero maximum (the default) derives an upper bound from the scaled
initial count or the explicit minimum, whichever is larger. Automatic limits
therefore permit thinning from the mesh-driven request; supply a larger maximum
to permit increases. Explicit bounds clamp the initial guess and subsequent
trials search within the chosen limits without multiplying again. Both targets
default to zero, disabling calibration and requesting the scaled count directly.
Existing input terminals are excluded from the count of new territories.
Containment, collision and growth constraints can bias the accepted terminal distribution away from the proposal density. The illustrated density plots use actual accepted physical terminals.
Numerical subdivision has separate controls. By default it targets interval length \(h\), with desired range \(0.5h\)–\(2h\), using the smallest diameter of background cells actually crossed by each segment. Short vessels and crossing fragments can prevent the lower bound; the discretization report records these exceptions. Refining these 1D cells changes neither physical topology nor integrated volume and resistance, and does not change which physical vessels are selected as growth candidates.
Arterial, portal and venous reference networks¶
Three distinct geometric networks can feed and drain the same tissue territory.
Use System::portal for the portal vessels and terminals, and
Port::portal_inlet for their pressure-controlled root. Physical vessel endpoint
order runs from each root toward tissue for all three systems; solved venous
flow has the opposite sign.
Register one terminal from each network with the same territory identifier and
positive tissue demand. A TerminalTriplet identifies those three terminals and
three positive resistances in arterial, portal, venous order:
std::vector<TerminalTriplet> triplets;
// Supply each territory once globally, on a rank of the caller's choice.
if (rank == territory_owner)
triplets.push_back({territory_id, arterial_terminal, portal_terminal,
venous_terminal, {{3.e11, 4.e10, 2.e10}}});
auto reference = solve_reference(network, {}, triplets, conditions,
reference_parameters, communicator);
Each resistance connects its terminal pressure to one common microvascular
pressure p_t. The equations are
Q_a = (p_a - p_t) / R_a
Q_p = (p_p - p_t) / R_p
Q_v = (p_t - p_v) / R_v
Q_a + Q_p = Q_v
The common pressure is eliminated from the distributed graph operator. It does
not create a geometric vessel, contribute vessel volume or split the territory
between the two supplies. The reported territory record contains
arterial_flow, portal_flow, venous_flow and compartment_pressure;
flow remains total venous drainage. Nonpositive supply or drainage invalidates
the result, including portal backflow caused by an incompatible pressure or
resistance choice. The supplied resistances and pressures are model inputs,
not patient-specific physiological estimates.
For an uncoupled demand-flow design solve, use the overload of
demand_flow_conditions() taking a std::map<System, double> with explicit
arterial, portal and venous totals. Choose Q_v = Q_a + Q_p. The scalar-flow
overload retains its arterial/venous behavior and rejects selected portal
terminals rather than silently assigning the full tissue flow to each supply.
The coupled reference solve determines the actual split from resistances and
pressures; prescribing desired totals does not calibrate those resistances.
The triplet overload also accepts ordinary two-terminal pairs, with unique
territory and terminal identifiers across both lists. Each triplet is supplied
once globally and resolved through sparse owner exchanges; terminal ghosts are
not required on the supplying rank. The generator can create the same three-system registration. Enable both
parameters.venous.enabled and parameters.portal.enabled, supply their seeds
or select imported terminal ids, and set a positive terminal separation for
each companion network. Arterial endpoints define the common territories;
portal and venous endpoints remain geometrically distinct. The same CCO engine
checks every generated vessel against all three systems.
parameters.reference_flow is total venous drainage. With portal enabled,
parameters.portal_flow_fraction assigns f Q to portal supply and (1-f) Q
to arterial supply, with 0 < f < 1. This fraction is a declared modeling
input. Each system distributes its rate over the same positive tissue demands,
with its own Murray exponent. Existing two-system configurations retain total
flow Q in each tree.
Leave parameters.compartment_resistances at its three zero defaults to derive
positive arms from a demand-flow design solve. For each territory the generator
chooses p_t = p_v + (min(p_a,p_p)-p_v)/2 and derives the three arms from the
requested flows. This midpoint is an illustrative resistance closure, not a
measured tissue pressure. An empty admissible pressure interval invalidates the
result. Alternatively provide three positive arms in arterial, portal, venous
order; when a reference flow is also requested, the coupled solve must realize
each system’s desired terminal flow within the reference tolerances. Keep the
legacy venous.microvascular_resistance at zero in three-system mode.
The resulting terminal_triplets belong to the arterial territory owner.
Use the save_network() overload accepting both pairs and triplets, or pass
both lists to the metadata-taking generate() overload when completing a saved
network. All three endpoints of an input compartment must be extended together;
untouched compartments keep their resistances. Three-system archives use format
3; two-system writes retain format 2, and both formats remain readable. Portal
and triple service-distance coverage are reported independently of the existing
arterial/venous joint coverage. Three-system calibration uses triple coverage.
The MetricFlowX transfer uses the same eliminated three-terminal resistor graph through native terminal connections. It preserves three distinct 1D networks and reverses only venous cells into the direction of flow. Portal inlet ports receive independent native boundary ids through the boundary-physics callback (with a null terminal pointer), rather than sharing the arterial inlet law. Provide their physical root boundary data explicitly in that callback.
Reference pressures, supply and microvascular beds¶
Both vascular_network and liver_vascular accept reference Robin conditions
under their vascularization Physiology subsection
(Vascularization/Physiology and Liver/Vascularization/Physiology,
respectively). Each entry maps a physical node to a conductance and a source:
set Reference Robin conditions = 30=5e-13:2e-9
The law is \(q_{\mathrm{out}}=\beta p-g\), with outward flow in \(\mathrm{m}^3/\mathrm{s}\), pressure in Pa, and conductance \(\beta\geq0\) in \(\mathrm{m}^3/(\mathrm{Pa}\,\mathrm{s})\). In this example the exterior pressure is \(g/\beta=4000\) Pa and the resistance is \(1/\beta=2\times10^{12}\) \(\mathrm{Pa}\,\mathrm{s}/\mathrm{m}^3\). Use commas between entries. Positive conductance fixes the pressure reference; each connected component needs either a positive Robin conductance or a prescribed pressure. Zero conductance is a prescribed inflow \(g\) and provides no pressure anchor. Robin data may be assigned to ports and unpaired tissue terminals, and cannot be combined with prescribed pressure or inflow at the same node. Paired tissue terminals retain their microvascular coupling and cannot also carry external Robin data. The reported boundary inflows include \(g-\beta p\); dissipated power continues to measure the physical vessels and microvascular beds.
With venous.microvascular_resistance == 0, supply is distributed by relative
territory demand, with total reference_flow through each system. The generator
first solves the two forests under the supplied port pressures and terminal
flows. It derives a separate bed resistance
\(R_i=(p_{A,i}-p_{V,i})/Q_i\) for every territory and requires it to be finite and
strictly positive. It then solves the connected graph using the supplied
external boundary conditions, with the derived beds and no prescribed flows
at paired tissue terminals. The actual bed flow must realize the demand split
within the configured numerical reference tolerance. No inlet receives both
an independent pressure and flow.
A positive venous.microvascular_resistance supplies that resistance directly
for each bed. With reference_flow == 0, pressures and resistances determine
the resulting flows. If a positive reference flow is also requested, the result
must realize its demand split; incompatible resistance/pressure/flow data fail
validation. A bed has no implicit capillary volume or geometric collision
segment.
result.reference reports actual pressures, signed vessel flows, territory
supply, endpoint wall shear, dissipated power and independent mass-balance
checks. target_terminal_flows and target_attachment_flows hold desired
flows, not substitutes for the solved values. The reference problem is steady
Poiseuille flow on complete physical vessels, including imported exterior
parts. It is separate from transient hemodynamics and accepted solver state.
solve_reference() is collective and uses native distributed vectors,
dealii::LinearOperator and CG with Jacobi preconditioning for one or more
processes. It assigns compact algebraic indices independently of physical ids
and numerical cells. Sparse endpoint incidence builds owned graph rows and
neighbor values; neither graph topology nor the solution is gathered.
Pressure offsets follow the nearest explicit anchor in resistive distance,
including across partitions, to retain small vascular drops beside large beds.
Supply each TerminalPair on exactly one process; that process receives its
territory flow. The solver resolves the two terminal identities and their
nodes even when the bed’s process owns no physical entities. Boundary maps may
be distributed or repeated with identical values. Conflicting duplicates and
pressure/flow overprescription fail collectively. Terminal owners must follow
their nodes, and endpoint ghosts must match their authoritative nodes.
Results contain owned vessel, terminal and boundary records, pressures at
locally relevant nodes, and global scalar diagnostics. reference_pressure()
evaluates owned physical vessels.
demand_flow_conditions() constructs collective conditions for an arterial or
unpaired arterial/venous design solve. It uses owned terminals, normalizes
arterial sinks and venous sources independently to the same positive total
flow, and returns conditions on their node owners. Pass its result directly
to solve_reference(). An empty selection includes all terminals; an explicit
selection must contain the same unique physical terminal ids on every rank.
Only explicitly supplied physical ports receive pressure data. Pressure maps
may be distributed or repeated identically. Demand scaling avoids overflow in
the global sum and retains representable flows even when an intermediate
relative demand would underflow. Nonrepresentable terminal flows are rejected.
The public generated-tree reference test exercises this complete path with
nonzero tissue demand and checks actual solved flows against the requested split.
A nonzero dissipation_weight blends generated-vessel volume and dissipated
power, each normalized by the current active forest value. Physical flow is
split globally across the active system’s roots by demand. Frozen imported
vessels contribute to the full reference solve but not to this growth objective.
Complete imported forests¶
Pass the existing NetworkData as initial. Set selected_terminals to
arterial tissue-terminal identifiers and venous.selected_terminals to venous
ones. Empty selections discover available interior leaves of the corresponding
system. If no venous attachment is available, supply a new venous seed radius,
position and viscosity. Multiple arterial attachments and multiple venous
outlets are supported. The territory count must accommodate every selected
attachment in each system.
Every input node and vessel retains its geometry, radii and attributes,
including input entries already marked as generated. Selected tissue-terminal
roles are consumed; attachment_roots maps each consumed input terminal to its
new root vessel. Unselected input terminal records remain unchanged. Explicitly
selected exterior terminals fail instead of extending new vessels outside the
background. Supply independent loads on unpaired input terminals through
reference_conditions; paired terminals cannot also receive external loads.
Volume fraction includes only inside portions of all 1D vessels, with additive volumes at junctions. It excludes beds. Coverage is a separate background quadrature estimate of finite service balls. In arterio-venous mode calibration uses dual coverage: quadrature points must lie within service distance of both systems. This geometric coverage is not a proof of physiological perfusion.
measure_tissue() is collective on the background communicator and returns
all metrics globally on each process. Only owners of terminal records supply
service centres. Centres travel to every background partition within service
reach, independently of the vessel ghost halo; a service ball can cross several
partitions. Each process integrates its own background cells, then the physical
quadrature weights are summed before normalization. Empty processes participate.
The supplied discretize() result already contains the global inside vessel
volume, which is used once. Setting the service radius to zero disables coverage
while retaining the global volume fraction.
Arterial terminal-count calibration is collective. Set a positive
calibration.volume_fraction or calibration.minimum_coverage to enable it.
Every process must supply the same targets, tolerance, terminal-count bounds,
attempt budget and initial terminal_count. Each attempt grows from the
unchanged input with native ownership and measures global volume and coverage.
The finite count search makes the same decisions on every rank without
collecting the graph. It does not assume that volume increases with terminal
count. The defaults are 5% relative volume tolerance and at most 12 attempts.
A volume match still fails when the requested coverage is missing.
When the fixed imported inside volume alone exceeds the upper target tolerance,
calibration reports incompatible before growing. Otherwise, exhausting the
allowed counts or attempt budget reports exhausted, which does not prove
incompatibility. A failed search can retain its best valid network for
diagnostics, but its overall validation remains invalid. Trial sample counts
are global and include discarded attempts. The supplied parameters are not
modified.
Imported arterial startup first refreshes native point/midpoint ownership and endpoint incidence. It validates terminal records at their node owners, adds missing non-port leaf records, and selects interior arterial endpoints. Explicit selections must name unique input terminal ids; an explicitly selected exterior leaf is rejected. Unselected and exterior input branches remain fixed, including input entries marked as generated.
Containing-cell queries for all selected leaves use sparse batches. New roots
inherit the incoming vessel’s distal radius, viscosity and extra physical
properties. Node owners propose roots in concurrent geometry rounds; the same
batch checks and joint commit used during growth certify them before adoption.
New identifiers use disjoint integer ranges. No artificial roots or full graph
replicas are introduced. attachment_roots records each consumed terminal’s
new root id on the consumed terminal’s native node owner; root-vessel ownership
may differ after midpoint redistribution.
The final arterio-venous pairing uses sparse exchanges independently of endpoint ownership. Each common territory has one writer and exactly one arterial and one venous terminal, with distinct positions and matching demand. Pair records and actual bed-flow results stay at the territory writer; the venous terminal need not be present there as a ghost. For derived beds, an unpaired design solve supplies endpoint pressures from their native owners. Only new tissue terminals enter the demand normalization. Boundary conditions on unpaired input terminals remain in both reference solves. The final paired solve must produce positive perfusion and the requested demand split; desired flows are never substituted for its measured flow. This finalization is collective, including processes with no vascular entities or background cells.
A new venous outlet receives its stable node id before arterial CCO insertions. For imported networks, selected venous leaves retain their incoming distal radius, viscosity and additional physical properties. Outlet owners query available common territories in batches; each territory producer returns one nearest candidate. Its writer grants competing requests deterministically. A losing outlet retries after committed claims are removed. A geometric failure tries the remaining Cartesian terminal offsets, then another territory. New branches pass the common geometry selection and joint commit before any territory claim is retained. Input vessels and completed arterial geometry remain fixed obstacles throughout venous startup and growth.
Save and reload physical data¶
Save a validated result in a new directory on a filesystem visible to all ranks:
#include <immersx/vascularization/io.h>
const auto archive_directory = output_directory / "physical-network";
save_network(archive_directory, result.network, result.territories,
result.terminal_pairs, communicator);
const auto restored = load_network(archive_directory, background, mapping,
parameters.clearance, communicator);
save_network() is collective. Each writer stores its owned vessels and
terminals, relevant endpoint nodes, and its supplied territory/bed records.
Files contain at most 1024 owned records of each kind, plus endpoint nodes.
Even an archive with one writer can therefore distribute its parts among many
readers without loading the entire physical graph on one process.
A versioned manifest is published only after every part closes successfully.
The target directory must not already exist. An interrupted or failed write
can leave parts without a manifest; load_network() rejects that directory.
Existing archives are never overwritten.
NetworkArchive contains physical network, territories and
terminal_pairs. Integer identifiers, coordinates, endpoint order, taper,
viscosity, additional named properties, arterial/venous classification,
provenance, port roles and normals, tissue demand and bed resistance survive the
round trip. Numerical subdivisions and solved reference states are recomputed.
No growth history or solver state is restored. To complete a restored network,
pass its physical metadata as well as its graph:
const auto completed = generate(background, mapping, restored.network,
parameters, communicator,
restored.territories, restored.terminal_pairs);
Every input vessel remains fixed. Select both endpoints of any existing bed
that will be extended, using the arterial and venous selected-terminal lists.
The consumed bed and territory are replaced by new tissue territories.
Untouched territories and beds remain in the valid result. Existing beds
require venous growth to be enabled and enter both the design solve and the
final paired solve. Their flows can change when new branches share upstream
vessels; their resistances remain fixed. The specified reference_flow and
returned desired-flow maps apply only to new terminals. Volume and coverage
still include the whole physical network. A half-extended bed is rejected
rather than losing its terminal relation. The same rules hold in terminal-count
calibration and on processes with no local metadata.
Readers need not have the writer count. They read old parts cyclically, check
unique authority and matching endpoint data, then use redistribute() against
the prepared background. Interior points and vessel midpoints acquire their
native owners. Exterior entities retain their assigned reader. Territories and
beds move to the restored arterial endpoint owner. Only file counts and coarse
background boxes are replicated; no process must gather the complete graph.
Empty writers and readers join all phases. Invalid records, missing or truncated
parts, unsupported versions and filesystem errors fail collectively.
Discretize across background partitions¶
discretize() is collective and accepts owned physical vessels with their
endpoint nodes, plus optional ghosts. It sends each owned segment to the
background partitions whose coarse boxes it intersects. Exact tetrahedral or
bilinear hexahedral crossings determine resolution; bounding-box overlaps alone
do not count as crossings. Interface portions contribute volume once.
Numerical intervals stay at their physical vessel owner. Each interval records its execution owner from midpoint containment, choosing the lowest containing rank on interfaces; exterior intervals retain the physical owner. Inside volume and short/exterior-cell counts are global on every rank, and the numerical cell budget applies to the whole network. Processes with no owned background cells or vessels still participate.
Pass a native fully distributed destination to build the line mesh collectively:
#include <deal.II/distributed/fully_distributed_tria.h>
parallel::fullydistributed::Triangulation<1, 3> vessels(communicator);
const auto references = create_triangulation(local_network, subdivision, vessels);
Here subdivision is the result of discretize(). Numerical cells migrate to
their execution owners. A sparse directory of physical endpoint and numerical
cut identities supplies vertex-adjacent ghost cells, including all branches at
a bifurcation. No complete physical graph or mesh is gathered. Compact native
cell and vertex indices are independent of physical identifiers, including
identifiers larger than the exact integer range of a double.
references maps every locally owned or ghost coarse cell to its physical
vessel and parameter interval. Use physical_reference(cell->id(), references)
to recover that interval; skip artificial cells on empty ranks. The physical
network and its attributes remain unchanged. To change distributed numerical
resolution, call discretize() again and build a new mesh: native
parallel::fullydistributed::Triangulation does not implement refinement flags.
The serial destination continues to support native binary refinement.
Measure geometry relative to the background¶
Use mesh_length_scale() to derive physical cell-diameter bounds collectively
from the prepared background. The four policies are meter, local mesh size,
min mesh size and max mesh size. A metre policy uses a one-metre reference;
the other constant policies use the global minimum or maximum active-cell
diameter. Only owned cells contribute to those bounds, including when some
processes own no cells.
#include <immersx/vascularization/mesh_scale.h>
const auto scale = mesh_length_scale(background, "local mesh size", communicator);
const auto measures = measure_scaled_vessels(background, mapping, local_network,
scale, communicator);
measures contains owned physical vessels. Each normalized_length integrates
ds/H over the exact cell crossings used by discretize(). An interface uses
the smallest diameter among all containing cells, including across MPI
partitions. minimum_scale, maximum_scale and the normalized radius bounds
account for every positive-length crossing and both endpoints of each tapered
piece. They do not replace the physical lengths or radii in the network.
Local mesh size is undefined outside the background. In local mode,
exterior_length records those portions, complete is false, and the
normalized measures describe only the inside portions. Constant policies also
normalize exterior portions. No local diameter is assigned to exterior input.
The measurement preserves the network, topology, radii and physical properties.
Pass the same prepared scale as the last argument of measure_tissue() to
measure mesh-relative service coverage. A quadrature point is served when its
nearest tissue terminal is within service_radius*H(point). In local mode,
H(point) is the active tissue cell diameter, not the diameter at the terminal.
Remote terminal routing uses the maximum possible reach. Arterial, venous and
dual coverage retain physical JxW weights, and overlapping service regions
count once. The overload without a scale uses physical metre distances.
Redistribute physical ownership¶
For an existing distributed NetworkData, redistribute() migrates owned
entities against a prepared parallel background and refreshes geometric ghosts.
It is collective, including processes with no physical entities or no owned
background cells. This operation does not generate branches. generate()
supports MPI startup from an empty arterial input and from imported arterial
terminals, together with real venous outlets or selected venous leaves in
arterio-venous mode.
#include <immersx/vascularization/distribution.h>
const auto relevant_network =
redistribute(background, mapping, local_input, clearance + search_reach,
communicator);
Each physical entity has one authoritative input owner. A rank needs its own
entities and the endpoint nodes referenced by its locally relevant vessels;
it does not need the complete network. Node and terminal owners follow their
physical points. The physical vessel owner follows its midpoint, matching the
centre-based policy of ImmersedRepartitioner. Native ParticleHandler search
checks actual background cells. Points on interfaces choose the smallest
containing rank; a bounding-box overlap alone does not assign ownership.
Exterior points and vessel midpoints retain their previous owner. A terminal
always follows its node.
A crossing physical vessel retains one identifier, one endpoint pair and its
original taper and material properties. It is not cut into artificial physical
branches at MPI interfaces. Numerical cell ownership is a separate decision.
Every retained vessel carries its endpoint nodes and any terminal records at
those endpoints. Wire messages retain integer identifiers without converting
them to particle properties of type double.
Ghost exchange routes a vessel to background covering boxes intersecting its segment bounding box expanded by the maximum endpoint radius plus the supplied halo. The halo can include wall clearance and an additional candidate-search reach. It can span several cells or partitions and is independent of mesh adjacency. These are conservative recipient queries; callers must still check actual segment geometry before accepting an insertion. Vessel ghosts also reach both endpoint owners, including exterior endpoints beyond every background covering box. This supplies complete incidence for input-leaf validation and attachment selection.
Authoritative node values replace stale ghost values before midpoint location.
Ghost edits to geometry, radii or demand cannot override owned data. Sparse
migration and ghost exchanges use Utilities::MPI::some_to_some; only coarse
background covering boxes are replicated. Duplicate owned identifiers are
checked through a distributed hash partition of the identifier registry.
Invalid ownership or missing authoritative nodes fail collectively at a phase
boundary. Input data remain unchanged. Repeating the operation after changes
replaces obsolete ghosts, including when vessel radii shrink.
Validate a batch of proposals¶
Before accepting a serial CCO insertion, the engine passes its changed vessel geometry through the same batch validation phase exercised by the MPI tests. Point location verifies all endpoints globally. Each participating process checks its owned physical boundary patches and relevant vessel obstacles; process interfaces are not physical boundaries. An open planar section must match at least one physical boundary face, and every process must certify separation from its remaining boundary patches.
A proposed replacement needs consent from the current vessel owner, including agreement with the recorded endpoint identities and positions. New vessel rays are also checked against other proposals. Eligible proposals use cost, producer rank and producer-local identifier as a deterministic priority. One round selects proposals whose priority is lower than every eligible conflicting neighbor; the result is independent, but need not be maximal. Invalid proposals cannot block valid neighbors. All geometric participants and writing owners receive the same final decisions through sparse exchanges.
This phase evaluates a common snapshot and changes no accepted network or history. The public generation loop uses its decisions to confirm or discard the provisional distributed batch after physiological and geometric checks.
Normalize generated forests¶
Final serial CCO radii and desired flows use a forest evaluation kernel that also operates on distributed ownership. It reconstructs incidence at node owners and exchanges subtree demand and moments from leaves to the declared attachment roots. A reverse sweep distributes each fixed attachment radius and total subtree demand. Each vessel owner then computes its Murray radius; terminal owners receive the desired flows. No process needs the full topology.
The mutable vessel set is explicit: input vessels remain frozen even when their recorded provenance is generated. Each mutable root retains its existing uniform radius. Arterial and venous systems use separate exponents and total demands. Positive reference flow is shared by terminal demand within each system; zero leaves desired flows unspecified. Common physiological controls must agree across processes. Invalid connectivity, unrooted cycles, missing terminal roles and unrepresentable positive demands or radii fail collectively.
Pass the prepared MeshLengthScale to evaluate_forest() with its physiological
controls. Its minimum reference supplies a necessary global minimum-radius
bound; exact along-branch bounds are checked during geometry certification.
One-child physical continuations, such as a venous seed stem, preserve the
Murray radius; branching junctions have two children.
The kernel changes no physical data. The public growth loop uses its result to normalize a tentative batch, then adopts the radii only after validation and the required geometric checks. Tests exercise joint ancestor and descendant splits across partitions through this same commit path.
Confirm or discard a batch¶
The internal batch staging path sends each proposed bifurcation to its physical vessel owner. The owner checks the recorded endpoint identities and positions, explicit mutability and exclusive writes. The serial CCO engine uses the same physical split operation. New vessels inherit the selected vessel’s material properties; frozen input vessels remain unchanged.
Tentative geometry migrates against the background with its private mutability and fixed-root controls. The forest kernel computes the joint demand and radii of all accepted splits, including ancestor and descendant splits in the same tree. Ghosts are refreshed for the new radii. Every mutable vessel is then checked on this final snapshot, using the same distributed geometry checks as proposal selection. A failed check discards the entire batch on every process. The caller can resample within its finite budget; it never adopts part of a failed batch. Staging copies locally relevant data once per batch, rather than once per candidate.
The public generation loop connects this commit path to distributed sampling and nearest-vessel selection. Certifying all mutable vessels after each batch is a correctness check; the growth statistics record its cost separately.
Optimize from distributed ancestor paths¶
The serial engine and distributed candidate evaluation use the same CCO cost and deterministic local search. A candidate needs its selected physical vessel and the path to its actual attachment root. Subtree demand, minimum leaf demand, and moments of the complete rooted tree provide the remaining information. Dissipation moments include each branch’s viscosity. Candidate evaluation changes only the trial ancestral path and never copies the complete network.
MPI path requests go directly to branch owners. Each reply identifies the next ancestor and its owner; subsequent requests are grouped by destination. Shared ancestors are retrieved once per requesting process. Processes without requests still participate and can serve their owned branches. Missing owners, cycles, and inconsistent endpoint geometry fail collectively. Only the union of the requested paths is retained, including integer physical identifiers and port normals. No tree is replaced by a local artificial root.
Geometry checks supplied during optimization are preliminary. A candidate from a local background partition can cross other partitions and must pass final owner consent, joint normalization and global geometry certification. An MPI test performs this sequence and compares the candidate cost with independently recomputed physical volume and dissipation. The public generation loop invokes this candidate evaluation after distributed sampling and nearest-vessel search.
Grow a prepared forest in MPI batches¶
The growth loop operates on already prepared real attachment trees with one or several processes. Arterial proposals use cached owned-cell masses and uniform physical-cell sampling through the native inverse mapping. Rank mass correction preserves the global proposal density without assigning equal terminal counts to differently sized partitions. Native containing-cell queries confirm ownership; shared interfaces use the smallest containing diameter and lowest containing rank. Spatial sampling weight and integrated tissue demand remain separate.
Nearest-vessel search uses owned physical segments. Only one covering box per nonempty vessel owner is replicated. A first sparse query obtains a distance bound; a second queries every remaining box that can improve the nearest candidate list. Replies contain physical identifiers, owners and exact point-to- segment distances. Ancestral paths are then fetched for these candidates.
Every producer chooses its best locally optimized candidate. Sparse geometric votes select independent proposals; shared ancestor demand updates do not lock the whole tree. Joint normalization and final geometry checks confirm or discard the selected batch. Provisional integer identifiers are reserved by process and namespace, independent of numerical cells. The loop never creates roots to activate idle processes. A Murray bound rejects targets that cannot fit the minimum radius under the fixed roots. Failed batches consume a finite consecutive sampling budget, and incomplete runs are reported as failures.
Tests start from one real seed, accept insertions from multiple producers in the same batch, and reproduce the result at a fixed partition and random seed. Empty background processes join the same phases without sampling. This loop serves both public arterial and venous generation and calibration. Distributed venous startup, exclusive common-territory assignment, positive bed derivation and actual paired reference validation use the public MPI path. Seed/import startup, terminal-count calibration, discretization, native mesh construction, tissue measurement and the resistive reference solve support MPI. The tutorial reports completed strong and weak scaling measurements, with three seeds per case and up to 10,000 terminals on eight processes.
Public arterial generation uses the batched engine with an explicit terminal count for each calibration attempt; it does not copy or re-register the parameter object. Leaf identities are returned on their terminal owners even when reference flow is unspecified, so geometric territories do not depend on whether a flow problem was requested.
Before local optimization, each producer retrieves vessel geometry from physical owners over the possible trial region and the fixed ancestor paths. This includes exterior input vessels whose owners may have no nearby background cells. Trial collision checks use proposed Murray radii, including shrinking siblings and growing ancestors. A colliding low-cost trial cannot displace an already feasible position. Preliminary batch votes do not compare these trials against obsolete snapshot radii; final confirmation always checks the complete jointly normalized geometry. A regression places a frozen obstacle outside an old parent’s wall and verifies rollback when an insertion makes that parent’s radius reach the obstacle.
The same execution loop grows veins from their actual outlet or attachment roots. Each common territory has one authoritative producer at the arterial terminal owner. That producer selects an offset endpoint and optimizes across nearby venous physical vessels. The offset may belong to another background process; it is not confined to the producer’s partition. Its demand is copied exactly from the arterial territory, without evaluating the tissue-demand expression again at the offset point.
A distributed identifier registry checks that common territories have unique owners. Existing venous leaves claim their territories through this registry; duplicate claims, missing territories and mismatched demand fail collectively. Successful new insertions remove only their producer’s pending territory and return the venous terminal identity to that owner, even if the terminal migrates to another process. No process needs the complete set of territory positions. This enforces one venous endpoint per territory in addition to geometric collision checks. The previous separate serial growth implementation has been removed. Public paired reference solves check positive beds and actual demand-matched flows through the same distributed execution path.
Standalone application¶
The vascular_network executable uses the same library implementation and needs
no MetricFlowX installation. In Debug builds its name is
vascular_network_debug. In a DebugRelease build both names are available.
Print the complete parameter interface into a new file:
build-debug/vascular_network_debug --write-parameters vascular.prm
The following synthetic SI configuration generates one physical vessel with a fixed radius. It prescribes inlet pressure and total tissue outflow at distinct nodes. It is an analytic CLI test case, not an organ calibration.
subsection Vascularization
subsection Input and ports
set Mode = generate
set Seed position = .2,.5,.5
set Seed direction = 1,0,0
set Seed radius = .005
end
subsection Growth
set Seed length = .4
end
subsection Sampling
set Terminal count = 1
set Terminal count multiplier factor = none
end
subsection Physiology
set Dynamic viscosity = .004
set Total new-tissue reference flow = 2e-9
set Reference pressures = 0:12000
end
subsection Objectives
set Maximum service distance = 2
end
subsection Output
set Directory = generated
end
end
Run it with one or several processes:
mpirun -np 2 build-debug/vascular_network_debug vascular.prm
Input archive and output paths are relative to the parameter file. The output
directory must be new. Missing parameter files are errors and are never created
implicitly. parameters.prm records the effective configuration with resolved
paths. Invalid generation or reference results return a nonzero exit code and
write summary.csv and errors-RANK.txt; they are not saved as valid archives.
Mode = complete requires Input archive = generated/network. It preserves
physical input vessels and extends the selected terminal ids. Supply Seed length for new attachment branches and Total new-tissue reference flow for
the new tissue. Viscosity and attachment radii can be inherited from the input.
Retained territories and paired beds follow the library completion rules above.
Mode = load requires the same archive path, performs no growth and does not
require seed parameters or viscosity defaults. Choose a new background or
Discretization/Target length factor to rediscretize the physical network.
For arterial networks without paired beds, Physiology/Loaded tissue reference flow applies a total demand-based outflow across all loaded terminals. Zero
uses only explicit reference boundary data. Input and ports/Evaluate reference
can disable this additional solve. Generated arterio-venous networks always
validate their paired reference solve.
Domain/Shape selects cube or sphere. Configure the cube origin and side
length or the sphere center and radius. Uniform refinements prepares the
background before growth. For a cube, Refinement focus = edge selects an edge
using Edge axis and two Edge sides (lower=0, upper=1, in transverse coordinate
order). For a sphere, Refinement focus = center selects its center. Each
Refinement fractions entry adds one normal deal.II refinement pass where cell
centers lie within that fraction of the side length or sphere radius. The
initial fractions are .5,.25,.125. The prepared background remains fixed
during generation.
Each rank writes owned physical nodes, vessels, terminals, territories and pairs
as CSV. Vessel tables include full physical volume, resistance and, when solved,
flow and shear. Interval tables describe numerical subdivisions separately;
background tables record owned cells, vertices and quadrature volumes. Separate
quadrature-RANK.csv tables contain native background quadrature positions and
physical weights for spatial statistics. Optional
native parallel VTU files include normal mesh ghost cells. Use owned CSV records
for counts and statistics. summary.csv contains global metrics, calibration
status and phase times. ranks.csv contains ownership counts and per-rank phase
times, resident memory and high-water resident memory in kB. Memory fields are
zero on platforms where deal.II does not provide them. The generation phase
includes its internal discretization, metrics and paired reference solve;
additional_reference_seconds measures only the optional subsequent solve.
Library callers can use ApplicationParameters and the collective run() in
immersx/vascularization/application.h for these prepared domains and outputs.
Supply absolute archive/output paths. Use generate() directly when the calling
application already owns its final background.
GenerationResult::growth_statistics records local growth work in ranges grouped
by the power of two of the leaf count before each batch. Each row identifies the
active arterial or venous system and requested terminal count. Calibration keeps
rows from all trials, including unsuccessful trials. first_terminals and
last_terminals are the actual observed counts; a concurrent insertion can cross
a range boundary. Sampling, proposals and accepted insertions count local
producers. Owned-vessel counts include only mutable vessels of the active system,
not geometric ghosts or frozen input vessels. Summing the first owned counts
across ranks therefore shows the real seed/attachment startup, rather than
assuming that every process owns a root.
The standalone writes these rows to growth-RANK.csv, even when physical table
output is disabled. Initialization covers preparation inside the forest-growth
call; seed/attachment construction before that call remains part of the overall
generation time. Local sampling, search, optimization, selection and commit times
include communication waits. Redistribution, physiology and certification are
subphases of commit and must not be added to commit again. Physiology measures
native distributed forest evaluation and radius updates. Certification includes
its geometry-index construction and joint owner votes. The data cover completed
batch attempts; an error that exits midway through a batch does not add a partial
row. Collection adds no graph gathering or collective reductions. It stores only
a logarithmic number of rows per system and calibration trial.
For communication measurements, configure with -DIMMERSX_BUILD_BENCHMARKS=ON
and build vascular_network_profile (or vascular_network_profile_debug). On
Unix MPI builds, these optional executables run the same standalone entry point
with standard PMPI interposition. The normal application and library do not link
the interposition code. Set IMMERSX_MPI_PROFILE_DIRECTORY to an existing,
absolute output directory before launching. Each process writes mpi-RANK.csv
before MPI finalization; existing files are never overwritten. An unset or empty
value disables profiling. A measurement is incomplete if any rank’s file is
missing or the application reports an output error.
The counters cover the explicitly listed functions in that CSV, from completed MPI initialization to the start of finalization, including application setup and output. They record successful calls, time inside those calls, point-to-point posted send bytes and collective input bytes. Bytes count logical datatype contents; they are not measured wire traffic. Collective input counts each process’s contribution, including the root’s broadcast/scatter input, without estimating the MPI implementation’s communication algorithm. Receive calls and waits contribute call counts and time, not additional send bytes. Nonblocking call time measures initiation; completion time appears in the corresponding wait/test. Unlisted MPI functions and MPI initialization/finalization are not included. Counter updates and datatype queries add profiling overhead to the application’s wall time. Compare measurements made with the same instrumentation.
The focused CTest tests ImmersX.VascularProfile.*.MPI2 compare counters against
known native MPI payloads, including in-place collectives and disabled forwarding.
The profiler adds no mandatory dependency and does not collect physical graphs.
Compare background refinements¶
scripts/run_vascular_mesh_study.py runs the ordinary C++ application at three
uniform and three edge-refined resolutions, each with random seeds 23, 24 and
25. It keeps the physical inputs and growth coefficients fixed. Only the domain
refinement, random seed and output location change. The supplied synthetic
unit-cube experiment declares its inlet position, 0.003 m root radius,
0.004 Pa s viscosity, 12000 Pa inlet pressure and 2e-9 m^3/s total flow in every
input file. These values are an experiment, not an organ calibration.
After configuring a Release build, run:
cmake --build build-release --target vascular_network -j
python3 -m venv /tmp/vascular-study-env
/tmp/vascular-study-env/bin/python -m pip install matplotlib
/tmp/vascular-study-env/bin/python scripts/run_vascular_mesh_study.py \
--executable build-release/vascular_network \
--library build-release/libimmersx_lib.so \
--output /tmp/vascular-mesh-study
The default mesh-based count, density exponent, seed length, minimum normalized
length, clearance, service distance and minimum radius apply without tuning
individual resolutions. common-parameters.prm records the shared controls;
the script checks that each case changes only its declared domain/seed/output
fields. It copies and hashes the executable and shared library and verifies the
library selected by the dynamic loader. The manifest also records source
identity, script hash, CPU affinity and relevant environment variables.
Each native case retains its parameters, log, physical archive, VTK meshes,
owned tables and growth work counters. An independent table oracle checks the
tree, positive integrated demands, unit tissue integral, total flow, Murray
radii and normalized minimum lengths. summary.csv reports per-seed quantiles,
counts, demand and flow per tissue volume, runtime and peak rank memory.
results.json retains the individual normalized lengths/radius bounds and
accepted terminal counts by cell diameter. Demand quadrature cardinalities are
global maxima; operation counts are rank sums, with initialization and
finalization removed in the ordinary-work totals.
The SVG/PNG distribution plots pool the declared seeds for each resolution. Their cell-diameter panels retain the individual seed counts. They show observed variation; the study imposes no statistical-equivalence threshold. Geometric constraints can change the accepted spatial distribution even when proposal weights are comparable per cell.
Use --levels 1 --families uniform --seeds 23 in a new output directory for a
small check before the full study. --ranks 2 or --ranks 4 exercises the same
protocol in MPI; use a separate output directory for each rank count.
--resume checks experiment identity and revalidates completed cases. It leaves
failed or incomplete cases available for inspection and requires an explicit
audit before restarting them. Refinement-study timings include the current
quadrature resolution and one-time final work; compare the operation counts and
reported quadrature sizes when interpreting them.
Reproducible MPI measurements¶
scripts/run_vascular_benchmarks.py launches the instrumented C++ application
sequentially across declared process counts and random seeds. It copies the
executable and ImmersX shared library into a new build-tree run directory so a
later rebuild cannot change an experiment already running. The manifest records
binary and script hashes, source revision and changes, MPI version, linked
libraries, CPU information, affinity and thread settings. Each case retains the
exact command, binding report, input and effective parameters, native tables,
per-rank profiles and independently checked measurements in case.json.
For example, with OpenMPI and a Release build:
cmake -S . -B build-release -DCMAKE_BUILD_TYPE=Release \
-DIMMERSX_BUILD_BENCHMARKS=ON
cmake --build build-release --target vascular_network_profile -j
OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 \
python3 scripts/run_vascular_benchmarks.py \
--executable build-release/benchmarks/vascular_network_profile \
--library build-release/libimmersx_lib.so \
--output build-release/vascular-strong \
--kind strong --terminals 1000,10000 --ranks 1,2,4,8 --seeds 23,24,25
OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 \
python3 scripts/run_vascular_benchmarks.py \
--executable build-release/benchmarks/vascular_network_profile \
--library build-release/libimmersx_lib.so \
--output build-release/vascular-weak \
--kind weak --terminals 1000 --ranks 1,8 --seeds 23,24,25
Use --mpiexec to select the MPI launcher matching deal.II. The default launch
arguments are OpenMPI’s --bind-to core --report-bindings; repeat
--mpi-argument=VALUE to replace them for another MPI implementation. Use enough
allocated cores and avoid competing jobs while measuring. These experiments can
be expensive: the per-case timeout defaults to four hours and is configurable
with --timeout in seconds. Use smaller terminal counts for a driver smoke test.
Strong scaling keeps the unit cube, 512 background cells, terminal count,
fixed seed radius 0.003 m, reference flow 2e-9 m^3/s and other parameters
unchanged across process counts. Weak scaling uses cube process counts
P=1,8,64,...: side length P^(1/3) m, 512P background cells, TP terminals
for the supplied base count T, flow 2e-9 P m^3/s, and seed radius
0.003 P^(1/3) m. This preserves background cell size, terminal density, flow
per terminal and Murray leaf radius across the weak cases. Each individual run
keeps its seed radius fixed and starts from exactly one real inlet at the same
relative position (0.2,0.5,0.5). The seed length remains 0.15 m, service
distance 0.2 m, viscosity 0.004 Pa s, density exponent zero, Murray exponent
three and inlet reference pressure 12000 Pa. These are declared synthetic
experiments, with fixed terminal counts and no volume-fraction calibration or
minimum-coverage target. The growing domain and single-root propagation make
weak scaling a measurement of this protocol, not an assumption of independent
per-rank trees.
Validation checks unique owned identifiers, connected acyclic single-inlet topology, exact leaf counts, Murray normalization, positive lengths, cube wall containment, junction angles, uniform terminal flow, pressure drops and exact volume/resistance preservation after subdivision. It also checks ownership and growth counters and requires an MPI profile from every process. The recorded quality distributions include physical vessel lengths, radii, junction angles and terminal distances from the cube centre. Coverage and volume fraction come from the native background integration. No claim of physiological adequacy is made from those geometric measures.
growth in case.json reports maxima of local phase times in each leaf-count
range and the fraction of rank-batches owning mutable vessels. The separate
growth_rank_totals retain local sums across ranges; use their maximum when
reporting a whole-growth phase time. Do not sum per-range maxima as if they were
a single rank’s elapsed time. mpi.seconds is a sum of rank times inside the
listed MPI functions and can exceed elapsed wall time. ranks retains per-rank
resident and peak resident memory. owned_vessel_max_over_mean measures final
physical ownership imbalance; it does not describe startup activity.
The benchmark publisher also reports cumulative per-rank sample, proposal and accepted-insertion counts and their maximum-to-mean ratios. These measure recorded producer work independently of final physical ownership. The startup proposer fraction counts ranks that produced at least one proposal in the first leaf-count range; it differs from the writer rank-batch fraction and does not measure CPU utilization. The published manifest hashes the native growth counter files used for these summaries.
A new run never overwrites an existing directory. Pass --resume with identical
arguments to recheck completed cases and continue an experiment. A case left
running, timed out or failed requires an explicit audit; the runner does not
infer that its processes have stopped or restart it automatically. Keep failed
runs and their diagnostics when reporting measured outcomes.
Prepare a generated network for native MetricFlowX¶
The generator has no MetricFlowX dependency. Applications that enable the
optional package can include <immersx/vascularization/metric_flow_x.h> and
call prepare_metric_flow_x() after generation succeeds. This requires a
MetricFlowX installation providing BloodFlowSystem::NetworkPhysics and
setup(const NetworkPhysics &). Older installations continue to support the
ordinary imported-input application path; selecting generation reports the
missing native API before growing a network. Build MetricFlowX with
-DMETRIC_FLOW_X_BACKEND=Trilinos for this checkout, whose ImmersX linear
algebra configuration selects Trilinos. Both packages must use the same backend
and deal.II installation. See the supporting native API in
MetricFlowX #86.
Construct the final background first, then generate, validate and prepare the
flow problem. For parameter_snapshot_file, pass an absolute build-tree filename
inside the selected output directory. The native overload creates its parents
and avoids a parameter snapshot in the working directory. For example, with a
configured native Square root wall model:
using namespace ImmersX::Vascularization;
using Flow = MetricFlowX::BloodFlowSystem<1, 3>;
Flow flow(comm);
flow.initialize_params(parameter_file, parameter_snapshot_file);
auto result = generate(background, mapping, {}, parameters, comm);
AssertThrow(result.validation.valid(), dealii::ExcMessage("Generation failed"));
auto transfer = prepare_metric_flow_x(
flow, result,
[](Id, const Vessel &) {
Flow::VesselProperties material;
material.E = 5.e5; // Pa, synthetic
material.h_wall = .0003; // m, synthetic
material.p_d = material.p0 = 0.;
return material;
},
[](Id, const Node &, const Terminal *) {
// R1, R2 [Pa s/m^3], C [m^3/Pa], outlet pressure [Pa], label.
return Flow::RCRPhysics{0., 1.e8, 0., 0., "synthetic"};
});
Both callbacks run only at physical owners. The first receives a physical vessel ID and its immutable attributes; geometry, length, endpoint radii, reference areas and viscosity come from the result. The second receives a physical node ID and a terminal pointer, or a null pointer at a real outlet. Arterial inlets retain native boundary zero and need no terminal callback. Material and wall choices must match the native configured model. These synthetic values illustrate the API; they are not an organ calibration.
For an explicitly closed anatomical end, pass its physical Terminal ID in
the optional fifth argument, a std::set<Id> of sealed terminals. Node IDs
and native boundary IDs are different namespaces. Entries may be distributed
or repeated across ranks. The helper rejects unknown IDs and terminals already
connected through a saved bed. It skips the boundary callback for sealed ends
and transfers an exact zero-flow condition, with no downstream reservoir or
capacitor state. This requires the sealed-terminal API from
MetricFlowX #88 and a
rebuild against its matching headers and library; an empty set also works with
older backends.
A sealed end is a declared model boundary. Record that choice and exclude it from coverage measures intended to describe supplied and drained tissue.
The helper constructs the native fully distributed line mesh with execution owners following the background. It exchanges coefficient tables, identifiers and endpoint boundary flags, without gathering the complete physical geometry or interior incidence. The native solver currently replicates coefficient and terminal tables. This is a remaining memory cost even though the vascular network and numerical mesh stay distributed.
Physical IDs remain unchanged and may exceed native material-ID capacity in
value. transfer.vessel_ids, boundary_ids and terminal_ids map them to
compact execution IDs. Venous cells and tapered radii are reversed into flow
direction in a temporary mesh input; the physical network retains its original
growth orientation. transfer.reference(cell->id()) returns the original
physical interval, including after native binary refinement. In a reversed
vessel, local coordinate zero corresponds to the interval’s end. Supplied
arteriovenous terminal pairs become native terminal resistive connections;
there is no second 0D model or extra line-cell volume.
The helper rejects a reference area at or below the native numerical floor
1e-10 m^2; it does not clip a physical radius. Native setup validates the
material law, terminal records and prepared mesh before distributing flow DoFs.
Use the native rest initial-state source. The native VTK-based steady-network
initializer is not available for this in-memory path. Subsequent native
consistent initialization and the ImmersX IDA adapter use the transferred data.
In metric_flow_x_elastodynamics, set /MetricFlowX coupling/Generate vascular network = true and provide /Vascularization/ parameters. The solid domain is
the generator’s background; standalone domain and output settings do not create
a second background in this application. Supply the coupling parameters
Generated vessel Young modulus, Generated vessel wall thickness,
Generated vessel reference pressure, Generated terminal resistance and
Generated outlet pressure in SI units. Young modulus, thickness and terminal
resistance must be positive. The terminal option supplies native R2 with R1 and
C zero; the library callback permits more general native RCR data. The vessel
reference pressure sets native p_d once; native p0 remains zero.
Generation and validation occur before native flow setup, multiplier-space construction, existing wall coupling and IDA initialization. Constant-radius networks use the native area basis for their reference field, so the initial wall displacement is exactly zero. When a physical vessel tapers, the reference degree is at least two to represent the squared linear radius exactly. The example uses fully internal walls. Adapting the multiplier to couple only internal portions of an imported external network remains deferred.
The configured fixture
gtests/parameters/metric_flow_x_generated_coupled.prm uses an internal
four-terminal tree, explicitly declared synthetic wall and terminal data, and
nonzero time-dependent inflow. The MetricFlowXGenerated application CTest
cases require accepted IDA steps, the requested output times, a nonempty
internal mesh, positive final flow change and positive solid displacement.
The fixture permits 30 nonlinear iterations per IDA corrector, retaining the
same error tolerances. HLL can need more than the imported example’s budget
of eight iterations when IDA reuses its Newton matrix. The pressure term
registers both its pressure-coefficient derivative and, with the supporting
native API, the flow-state derivative of HLL weights at fixed pressure.
The separate MetricFlowX.Vascular integration cases exercise tapered arterial
and reversed venous inputs, stable large physical IDs, native terminal beds and
nonzero consistent inlet flow on one, two and four ranks. Tests requiring the
new native API are skipped explicitly with an older installed package.
Publish measured MPI figures¶
After all planned cases in each experiment have completed, regenerate the selected benchmark assets from their native outputs:
python3 scripts/update_vascular_benchmarks.py \
build-release/vascular-strong build-release/vascular-weak
python3 scripts/update_vascular_benchmarks.py \
build-release/vascular-strong build-release/vascular-weak --check
Supply multiple experiment directories when strong terminal counts were run
separately. A repeated case name is rejected. Each plotted group and seed must
have its one-rank reference case. --destination /path/to/build/figures
creates draft assets outside the documentation. The publisher rejects a
destination that overlaps an input experiment, so it cannot overwrite native
measurement manifests.
The publisher verifies the frozen executable and library hashes, requires every
planned case to be complete, reruns the independent native-table checks, and
compares the recomputed measurements with each saved case.json. Only then
does it render SVGs and the small measurements.csv and manifest.json
summaries. Full meshes, tables, profiles and logs remain in the build tree.
--check rerenders and compares selected asset bytes against the destination;
it rejects a changed figure or changed native measurements. Use the same
measured run directories for this check. A new benchmark execution yields new
timing samples for review and publication, rather than identical timing figures.
Curves show medians across the supplied seeds; whiskers span their complete
observed ranges. The CSV retains each seed and physical-quality quantiles.
Strong and weak protocols are labeled separately. In strong scaling the plotted
T(1)/T(P) is a speedup for this distributed generation protocol. Concurrent
batching changes geometry, accepted work per batch and the number of global
physiological updates across partitions, so it is not an identical-work hardware
scaling experiment. In weak scaling the same ratio compares increasing global
problem sizes under the declared single-inlet protocol.
Resource figures report the maximum per-rank peak RSS, physical ownership imbalance, summed logical posted-send and collective-input bytes, summed rank time inside the selected MPI functions, the maximum rank physiological-update time and startup writer ownership. Logical bytes are not wire traffic; MPI rank time is not elapsed wall time. The sum of per-rank peaks in the CSV is not a simultaneous node-memory measurement and double-counts shared pages. Startup writer ownership counts mutable vessel owners, not every rank that samples or produces proposals. Quality figures show observed volume, coverage, physical lengths and angles, without fitted physiological acceptance thresholds or a claim of statistical equivalence from a small seed sample.
An instrumented forest-growth loop also records growth-mpi-RANK.csv in that
profiling directory. Each row identifies a local run, system, requested count,
power-of-two terminal range and initialization, batch or finalization
phase. Run identifiers follow the collective growth-call order, including
calibration attempts. A zero terminal-range lower bound marks initialization
entered before the accepted count is known. Nested finalization is charged to
finalization alone. Ending initialization leaves no active fallback phase between batches.
posted_sends counts successful point-to-point send calls, including zero-byte
messages and self sends, but excludes MPI_PROC_NULL. Collectives do not add
estimated message counts. Their logical input bytes remain separate from posted
send bytes; add the two byte columns for total recorded logical payload.
reductions counts the recorded reduce/scan calls, including MPI_Iallreduce.
Scope collection adds no MPI calls, reductions or extra serialization. It records
only the MPI functions covered by the profiler, not an estimate of wire traffic.
The phase file is absent when profiling is disabled. A header without rows means
that the executable opened no growth scopes. Seed preparation and other work
outside the forest-growth call remain in mpi-RANK.csv. Compare its totals with
the phase sums to identify that work. Phase data live in this optional profiler;
they are not zero-valued communication counters in GrowthWork.