Zone-Based Domain Decomposition

Zone-based domain decomposition is the primary strategy Sentinel.jl uses for solving the inverse problem in brain Magnetic Resonance Elastography (MRE). This page explains why zones are methodologically essential, how they are generated and solved, and how zone-level results are consolidated into a global material property map.

Why Zone Decomposition?

In brain MRE, displacement measurements are acquired volumetrically via phase-contrast MRI. The full-brain finite element mesh may contain tens of thousands of elements, but the inverse problem – recovering the spatial distribution of viscoelastic material properties from measured displacements – is ill-posed in ways that zone decomposition directly addresses.

Boundary condition quality. The brain surface presents noisy, incomplete displacement measurements due to partial-volume effects, CSF pulsation, and skull proximity. A single full-brain inverse solve would require prescribing Dirichlet boundary conditions at these unreliable surface nodes. Zone decomposition sidesteps this by partitioning the brain into overlapping subregions (zones) whose boundaries lie in the interior of the brain, where measured displacements are reliable. Each zone's boundary conditions come from neighboring measured displacements, not from the problematic brain surface.

Well-posed subproblems. Each zone is small enough that its boundary-to-interior ratio is favorable: the measured displacement data at interior nodes provides strong constraints on the material parameters, while the prescribed boundary displacements anchor the forward solve. This conditioning is far better than a monolithic solve where the ratio of unknowns to reliable constraints is unfavorable.

Iterative refinement through overlap. Because zones overlap, each material node is typically updated by multiple zones. The GP-weighted averaging in the consolidation step acts as a natural smoothing operator, suppressing zone-boundary artifacts and improving robustness across global iterations.

Not Just Computational Convenience

Unlike classical domain decomposition methods (e.g., Schwarz methods) that partition problems purely for parallelism, zone decomposition in MRE is methodologically essential. Even on a single processor, zones produce better reconstructions than a monolithic solve because they create better-conditioned subproblems with reliable interior boundary conditions.

Zone Generation

Zone centers are placed on a regular 3D grid that covers the mesh bounding box, then pruned to remove empty zones.

Algorithm

1. Seed selection. A seed element determines the grid origin. If seed_element = 0, the mesh bounding box center is used. When randomize_zones = true, each global iteration picks a new seed element via the Park-Miller RNG, shifting the entire zone grid.

2. Regular grid layout. Zone centers radiate outward from the seed with spacing equal to zoneedge in each dimension. The grid extends to cover the full mesh extent (computed from node coordinates, not centroids, to avoid missing boundary elements).

3. Pruning. A zone is kept only if at least one element centroid falls within its core box ($\pm$ zoneedge/2 from center, no overlap). Element assignment for the actual solve uses the extended box that includes the overlap region.

4. Zone growing. Each zone is tested for sufficient interior DOF nodes. The target is numsig = round((N_dof - N_boundary) / N_zones). If a zone has fewer interior nodes than numsig, its overlap is iteratively scaled by ${1.5}^{1/3}\approx 1.145$ until the threshold is met.

Zone Box Geometry

For a zone centered at $\mathbf{c}$ with edge size $\mathbf{e}$ and overlap fraction $\alpha$, an element with centroid $\mathbf{x}$ belongs to the zone if:

$$|x_d-c_d|\le \frac{e_d}{2}+\alpha \cdot e_d,d=1,2,3$$

With the typical 15% overlap ($\alpha =0.15$), each zone extends 65% of the edge size from its center in each direction.

Deterministic Randomization

When randomize_zones = true, zone seeds are generated by a Park-Miller "minimal standard" RNG with Bays-Durham shuffle (ParkMillerRNG). This is a direct port of the Fortran random(idum) function (Numerical Recipes RAN1), initialized with idum = -1 followed by one throwaway call. The deterministic seed sequence ensures exact reproducibility across Julia and Fortran implementations.

Zone Solve Pipeline

Each zone is an independent inverse subproblem. The pipeline extracts a self-contained problem from the global mesh, solves it, and returns updated material properties.

Zone Extraction Details

The build_zone function performs these steps:

1. Grid extraction. A new Ferrite Grid is created containing only the zone's elements, with renumbered nodes and cells.

2. DofHandler setup. A fresh hex27 DofHandler provides zone-local DOF numbering.

3. DOF coordinate matching. Hex27 elements have 27 DOF positions (corners, edge midpoints, face centers, body center), but Ferrite grids store only 8 corner nodes. Zone-to-global DOF mapping is built by matching physical coordinates via a spatial hash table (build_dof_coord_lookup) for O(1) lookup per node.

4. Material mesh extraction. Material meshes are subsetted by following GP2Mtr links from zone elements, not by spatial extent. This ensures correctness when the material mesh is coarser than the displacement mesh.

5. Boundary classification. A DOF node is classified as boundary (Dirichlet) if any of its connected elements in the global mesh is absent from the zone, or if it is a global mesh boundary node (from the .bnd file). See classify_zone_dof_nodes.

6. Boundary conditions. Boundary DOF nodes are prescribed with measured displacement values, providing reliable constraints from MRI data.

Zone-Level Optimization

The zone inverse solve dispatches to the configured optimizer via solve_zone!:

Optimizer	Symbol	Method
Conjugate Gradient	:cg	conjugate_gradient!
Gauss-Newton	:gn	gauss_newton!
Quasi-Newton (L-BFGS)	:lbfgs	quasi_newton!

CG is the default and most commonly used optimizer for brain MRE.

Consolidation

After all zones in a global iteration are solved, their results must be merged back into the global material and displacement arrays. The consolidation pipeline handles overlapping contributions through Gauss-point-weighted averaging.

Three-Phase Process

1. Prepare (prepare_accumulation!): Zero all property values, Gauss-point counts, and zone-coverage counts in the global material. Zero displacement values and counts.

2. Scatter (scatter_zone_material!, scatter_zone_displacement!): For each zone result, accumulate zone-local property values into the global arrays. Each node's contribution is weighted by the number of Gauss points that map to it in that zone (the rgp/igp counts). Displacement components are summed with a simple count.

3. Consolidate (consolidate_global!): Divide accumulated values by their total Gauss-point weights (material) or counts (displacement) to produce the weighted average. Nodes not covered by any zone retain their values from the previous global iteration.

GP-Weighted Averaging

For a material node $n$ covered by zones $z_1,\dots ,z_k$, the consolidated value is:

$${\overline{m}}_n=\frac{\sum _{j=1}^k{g}_n^{(z_j)}\cdot m_n^{(z_j)}}{\sum _{j=1}^k{g}_n^{(z_j)}}$$

where $m_n^{(z_j)}$ is the property value from zone $z_j$ and $g_n^{(z_j)}$ is the Gauss-point count for node $n$ in zone $z_j$. Interior nodes (more Gauss points) naturally receive higher weight than boundary nodes, biasing the average toward the better-constrained zone solutions.

Iteration Structures

The zone decomposition solver uses escalating iteration structures that increase computational effort as the global iteration progresses. Early iterations use few optimizer steps (cheap exploration), while later iterations invest more steps (refinement).

Each structure specifies the number of CG iterations, GN iterations, QN iterations, and line search iterations per zone solve. The active structure is selected by comparing the current global iteration against a set of threshold limits.

Typical Configuration

A common three-structure configuration for brain MRE:

Structure	Global Iterations	CG Iters	LS Iters	Purpose
1	1–10	1	1	Coarse exploration
2	11–150	2	2	Intermediate refinement
3	151+	3	2	Fine convergence

This is configured in the inverse runfile as ZoneIterationStructure entries with corresponding structure_limits. The solver selects the appropriate structure at the start of each global iteration via the parsed InverseRunfileConfig.

Global Convergence

At the end of each global iteration, convergence is assessed by the relative L1 material change:

$$\epsilon =\frac{\sum _i|m_i^{(k)}-m_i^{(k-1)}|}{\sum _i{m}_i^{(k-1)}}$$

computed by material_epsilon. The solver converges when $\epsilon <$ global_tol (after at least two iterations). The global_tol parameter in the runfile is specified as a percentage and divided by 100 internally.

Configuration

Zone decomposition is configured through ZoneConfig for programmatic use, or via the inverse runfile for production runs.

ZoneConfig Parameters

Parameter	Type	Default	Description
zoneedge	SVector{3,Float64}	[1, 1, 1]	Zone box half-extents in each direction (meters)
overlap	Float64	0.5	Overlap fraction (0.0–1.0); overlap region = overlap * zoneedge
optimizer	Symbol	:cg	Zone-level optimizer (:cg, :gn, :lbfgs)
max_zone_iter	Int	10	Maximum optimizer iterations per zone
zone_tol	Float64	1e-6	Convergence tolerance for zone solve
max_global_iter	Int	20	Maximum outer iterations
global_tol	Float64	1e-4	Global convergence tolerance (relative L1 change)
seed_element	Int	0	Element index to seed zone grid (0 = mesh center)
randomize_zones	Bool	false	Regenerate zone grid each global iteration
max_ls_iter	Int	1	Maximum line search iterations per CG step
line_tol	Float64	0.01	Line search convergence tolerance

Entry Points

Two driver functions are available:

• zone_decomposition_solve! – Serial zone loop. Solves zones sequentially in a single process.

• zone_decomposition_solve_parallel! – Accepts a pmap_func argument. Pass Distributed.pmap for multi-process parallelism, or use the default map for serial execution.

For production use with real brain MRE data, the recommended entry point is run_inverse_solve!, which handles runfile parsing, escalating iteration structures, reconstruction indicators, and post-processing regularizations automatically.

See Also

• Mathematical Reference – Governing equations and element kernels

• Zone Decomposition – Full API reference for src/distributed/

• Running an Inversion – User guide for configuring and running inversions