Inverse Solver Pipeline

Overview

The inverse solver in Sentinel.jl reconstructs spatially varying material properties (e.g., shear modulus, loss modulus) from measured displacement data. The algorithm uses zone-based domain decomposition: the global mesh is partitioned into overlapping subdomains ("zones"), each solved as an independent inverse problem, then consolidated via Gauss-point-weighted averaging.

The end-to-end pipeline has three stages:

1. parse_inverse_runfile – Parse a Fortran-compatible .dat runfile into an InverseRunfileConfig

2. setup_inverse_problem – Load mesh, displacement, and boundary files; generate material meshes; build regularizations; return an InverseProblemSetup

3. run_inverse_solve! – Execute the global iteration loop with zone decomposition

The convenience function run_inverse_solve! chains all three steps:

material, disp, state = run_inverse("runfile.dat"; basedir="/data/")

The setup phase (setup_inverse_problem) performs several important preparatory steps:

• Loads the displacement mesh (.nod + .elm files) and builds a Ferrite DofHandler

• Generates structured hex8 material meshes at configured resolutions, with GP-to-material-node mappings (GP2Mtr)

• Loads measured displacement from .dsp files, mapping between Ferrite DOF ordering and mesh-file node ordering

• Reads boundary node sets from .bnd files and builds Dirichlet BCs from measured displacement

• Initializes Gauss-point counts for consolidation weighting

• Constructs all regularization objects from the runfile configuration

Global Iteration Loop

Each global iteration regenerates the zone grid (with a randomized seed element to reduce partition bias), solves all zones, consolidates, applies post-processing regularizations, and checks convergence.

Zone grid regeneration. At each global iteration, the zone grid is regenerated with a different seed element drawn from a Park-Miller RNG. This prevents systematic bias from any single partition and ensures all regions of the mesh are eventually covered by zone interiors (where reconstruction quality is highest).

Iteration structures. The runfile specifies one or more ZoneIterationStructure entries, each defining the number of CG, Gauss-Newton, or L-BFGS iterations per zone. The active structure changes at configurable global iteration thresholds via structure_limits. Typical configurations start with fewer CG iterations and increase as the solution matures.

Parallelism. Zone solves are embarrassingly parallel. Pass pmap_func=Distributed.pmap to run_inverse_solve! for distributed execution, or use Julia's Threads via a threaded map. The default is serial map. BLAS threads are reduced to 1 during zone solves to prevent thread contention and ensure deterministic floating-point results.

Zone boundary conditions. Each zone's boundary nodes receive Dirichlet BCs from the measured displacement (not the previous iteration's calculated displacement). This prevents forward-solve errors from propagating into zone BCs and causing divergence. Interior nodes are free DOFs whose displacement is determined by the forward solve.

Callbacks. An optional callback(giter, material, state) function is called after each global iteration, enabling per-iteration recording (e.g., writing VTK snapshots or logging to file).

Per-Zone CG Optimization

Inside each zone, the optimizer (typically Polak-Ribiere-Plus conjugate gradient) minimizes the displacement misfit between the forward solution and measured data. The inner loop for a single zone proceeds as follows:

Key implementation details:

• evaluate_objective! assembles the global stiffness matrix from current material properties, applies Dirichlet BCs (measured displacement on zone boundary nodes), runs the forward solve, and computes the L2 displacement misfit plus regularization penalties.

• The adjoint solve reuses the cached LU factorization from the forward solve. Since $\mathbf{K}$ is complex symmetric ($\mathbf{K}={\mathbf{K}}^T$), the adjoint system $\overline{\mathbf{K}}\mathit{\lambda }=-({\mathbf{u}}_{\text{calc}}-{\mathbf{u}}_{\text{meas}})$ uses $\overline{\mathbf{K}}={\mathbf{K}}^H$, which is the conjugate of the already-factored matrix.

• Beta computation uses the raw (unscaled) gradient. The CG direction is also computed from raw gradients. Only after beta and the direction are formed is per-property scaling applied.

• The secant line search approximates the directional derivative ${\phi }^{\prime }(\alpha )=\mathrm{\nabla }f(\mathbf{x}+\alpha \mathbf{d}{)}^T\mathbf{d}$ using finite differences with a probe step $\sigma =0.1$. The secant update is:

  $$\mathrm{\Delta }\alpha =-\sigma \cdot \frac{{\phi }^{\prime }(\alpha )}{{\phi }^{\prime }(\alpha +\sigma )-{\phi }^{\prime }(\alpha )}$$

  Each secant iteration requires a full forward + adjoint + gradient evaluation. If secant fails (negative step on first iteration, or zero denominator), the solver falls back to brute force line search, which starts at $\alpha =1$ (capped by CG threshold), then tries $1.5\alpha$ on improvement or $0.5\alpha$ on worsening, for up to 10 iterations. Both line searches track the minimum-error solution encountered across all trial steps.

• BoundsReg is applied before every forward solve during line search to prevent physically invalid material values (e.g., negative stiffness), which would corrupt both the objective landscape and gradient computation.

Optimizer selection

The optimizer is selected per iteration structure: CG iterations > 0 uses conjugate_gradient!, GN iterations > 0 uses gauss_newton!, QN iterations > 0 uses quasi_newton! (L-BFGS). For non-CG optimizers, a fresh forward solve is performed after optimization to obtain the calculated displacement for consolidation.

Consolidation

After all zones are solved, their results are scattered to global accumulators and then consolidated by GP-weighted averaging.

The consolidation uses Gauss-point counts as weights. Each material mesh node accumulates a weighted sum of zone contributions, where the weight is the number of Gauss points that map to that node within that zone. This naturally gives more influence to zones where a node is well-covered by the FE mesh, rather than treating all zone overlaps equally.

For material nodes not covered by any zone in the current iteration (e.g., near the mesh boundary), the previous iteration's values are preserved. Similarly, uncovered displacement nodes retain their previous values, matching the Fortran behavior where globcalc is never zeroed between iterations.

See scatter_zone_material!, scatter_zone_displacement!, and consolidate_global! for the implementation.

CG Scalar Weights

CG scalar weights (enabled by regind(7) in the runfile) normalize the search direction so that step sizes are comparable across material properties with very different magnitudes. Without this scaling, a single step size $\alpha$ would over-update soft properties (shear modulus ~3300 Pa) and under-update stiff ones (bulk modulus ~1.6e6 Pa).

The scaling works per property group:

1. Accumulate the total absolute CG direction $\sum _k|d_k|$ and the total absolute material value $\sum _k|{\theta }_k|\cdot s$ (where $s$ is the property scalar) for each property group

2. Compute the scale factor: $\text{scale}=w\cdot \overline{\theta }/\overline{d}$ where $w$ is the CG weight, $\overline{\theta }$ is the average absolute material value, and $\overline{d}$ is the average absolute direction

3. Multiply all direction entries for that property group by the scale factor

The weights are interpolated from initial to final values over global iterations:

$$w(\text{iter})=w_{\text{init}}+\frac{\text{iter}-1}{\text{max\_iter}-1}\cdot (w_{\text{final}}-w_{\text{init}})$$

This schedule uses max_global_iter from the runfile (not any early-stopping override), ensuring the weight ramp is never compressed.

The first property's real CG weight also serves as the CG threshold – a cap on the maximum relative material change per step. During both secant and brute line search, trial steps are rescaled if material_epsilon(current, trial) > cg_threshold.

Typical values

For brain MRE with an isotropic incompressible model (Model 1), typical CG weight initial values are 0.01–0.05 for shear modulus and 0.001–0.01 for loss factor, ramping down to smaller final values over ~10–20 global iterations.

See scale_cg_direction! for the implementation.

Post-Processing Regularizations

Post-processing regularizations are applied at the global level after consolidation, not within individual zone solves. This is methodologically important: zone-local material arrays are subsets of the global mesh, and operations like spatial filtering require the full global neighborhood structure.

The only regularization applied at zone level is BoundsReg (hard property clamping during line search).

Spatial Filter (SpatialFilterReg)

Gaussian smoothing on the structured material mesh with kernel widths $({\sigma }_x,{\sigma }_y,{\sigma }_z)$. The filter widths are interpolated from initial to final values over global iterations, allowing aggressive smoothing early (to suppress noise) that relaxes as the solution converges. A sensitivity mask restricts filtering to material nodes covered by at least one Gauss point, preventing uncovered exterior nodes from diluting interior values.

Van Houten (VanHoutenReg)

Outlier suppression that clips power-law exponents to a configured level. Applied after spatial filtering.

Bounds (BoundsReg)

Hard clamping of property values to physically meaningful ranges. For example, shear modulus is clamped to $[{\mu }_{\text{min}},{\mu }_{\text{max}}]$, and derived Poisson ratios are constrained to $[{\nu }_{\text{min}},{\nu }_{\text{max}}]$. Applied last in the post-processing chain to guarantee all output properties are within bounds.

The order of application matters and matches the Fortran implementation:

1. Spatial filter (smooth)

2. Van Houten (suppress outliers)

3. Bounds (clamp to valid range)

Convergence

Global convergence is measured by material_epsilon, the relative change in material properties between consecutive global iterations:

$$\epsilon =\frac{\sum _i|{\theta }_i^{\text{old}}-{\theta }_i^{\text{new}}|}{\sum _i{\theta }_i^{\text{old}}}$$

where the sums run over all property values (real and imaginary parts, all properties, all mesh nodes). The denominator uses raw (not absolute) values, matching the Fortran implementation.

The convergence tolerance in the runfile is specified as a percent (e.g., 0.1 means 0.1% = 0.001 internal tolerance). Convergence is only checked after global iteration 1 (the first iteration always proceeds).

Within each zone, convergence is checked using the same material_epsilon metric against zone_tol. A zone CG loop also terminates early if the line search fails to find a descent step (alpha == 0).

See material_epsilon for the implementation.

Key Functions

Function	Location	Role
run_inverse_solve!	src/io/inverse_setup.jl	Top-level global iteration driver
run_inverse_solve!	src/io/inverse_setup.jl	Convenience wrapper: parse + setup + solve
setup_inverse_problem	src/io/inverse_setup.jl	Load files, build meshes, regularizations
evaluate_objective!	src/optimization/common.jl	Assemble K, forward solve, compute misfit
compute_full_gradient!	src/optimization/common.jl	Adjoint solve + gradient assembly
conjugate_gradient!	src/optimization/conjugate_gradient.jl	PR+ CG optimizer loop
secant_line_search	src/optimization/secant_line_search.jl	Secant method step size selection
brute_line_search	src/optimization/secant_line_search.jl	Brute force fallback line search
scale_cg_direction!	src/regularization/cg_update_scale.jl	Per-property direction scaling
consolidate_global!	src/distributed/zone_consolidation.jl	GP-weighted averaging of zone results
material_epsilon	src/types/properties.jl	Relative material change metric
build_regularizations	src/io/inverse_setup.jl	Construct regularization objects from config
build_zone	src/distributed/zone_extraction.jl	Extract zone-local problem from global state
solve_zone!	src/distributed/zone_solve.jl	Solve a single zone's inverse problem
scatter_zone_material!	src/distributed/zone_consolidation.jl	Accumulate zone material into global arrays
generate_zone_grid	src/distributed/zone_grid.jl	Generate overlapping zone partition