Architecture

Overview

Sentinel is a Julia port of MRE-Zone v9.37, a Fortran 90/95 nonlinear inverse solver for Magnetic Resonance Elastography (MRE). The architecture uses Julia's multiple dispatch system to replace Fortran if/elseif chains with type-based dispatch on material model singletons.

Data Flow

MRI Data (DICOM/MAT)
    │
    ▼
Preprocessing ─── MRI Readers (Siemens, WashU, UIUC)
    │                  │
    │              DTI Processing (fiber directions)
    ▼
Mesh Generation ─── HexMeshConfig → HexMeshResult
    │                  │
    │              Material Mesh (hex8, GP2MTR mapping)
    ▼
Forward Problem ─── Assembly → K (complex sparse)
    │                  │
    │              Dirichlet BCs → modified K, f
    │                  │
    │              Linear Solve → u_calc
    ▼
Inverse Problem ─── Objective: Φ = ½‖u_calc − u_meas‖²
    │                  │
    │              Adjoint Solve → λ
    │                  │
    │              Gradient: ∂Φ/∂θ (adjoint method)
    │                  │
    │              Regularization (10 methods)
    │                  │
    │              Optimizer (CG / GN / L-BFGS)
    ▼
Zone Decomposition ─── Zone Grid (random seed placement)
    │                      │
    │                  Zone Extraction (overlapping subdomains)
    │                      │
    │                  Zone Forward/CG Solve (per-zone)
    │                      │
    │                  Consolidation (GP-weighted averaging)
    ▼
Postprocessing ─── VTK Export
                   Convergence Analysis
                   Property Interpolation

FEM Infrastructure

Sentinel uses a hybrid approach with Ferrite.jl:

Ferrite.jl provides:

• Element type definitions and shape function evaluation (CellValues)

• Gauss quadrature rules

• DOF management (DofHandler)

• Sparsity pattern generation and sparse matrix allocation

• Assembly infrastructure (start_assemble, assemble!)

Sentinel provides (custom on top of Ferrite):

• Complex-valued element integration kernels for each material model

• Frequency-dependent material property evaluation at Gauss points

• Multi-mesh material parameterization (GP2MTR mapping)

• Adjoint gradient routines

• Inverse problem optimization loop

• All regularization methods

• Preprocessing (mesh generation from MRI, DTI processing)

• Postprocessing (VTK, convergence, property interpolation)

Material Model Dispatch

abstract type AbstractMaterialModel end
struct IsotropicIncompressible <: AbstractMaterialModel end   # model=1
struct Orthotropic <: AbstractMaterialModel end               # model=2
struct IsotropicCompressible <: AbstractMaterialModel end     # model=3
struct Poroelastic <: AbstractMaterialModel end               # model=4
struct TransverseIsotropicV1 <: AbstractMaterialModel end     # model=5
struct TransverseIsotropicV2 <: AbstractMaterialModel end     # model=6
struct GeneralizedAnisotropic <: AbstractMaterialModel end    # model=7

Each model dispatches to specialized methods for:

• element_stiffness! — element-level stiffness integration

• compute_gradient! — adjoint gradient computation

• has_pressure_dofs / pressure_dof_type — DOF structure queries

Module Organization

src/
  Sentinel.jl              # Top-level module, exports, includes
  types/
    materials.jl           # AbstractMaterialModel and 7 concrete types
    properties.jl          # Material, SingleProperty, MaterialGradient
    solution.jl            # Displacement, NodeDisplacement, NodePressure
    problem.jl             # Problem, BoundaryConditions, BCType
    forward_problem.jl     # ForwardProblemSetup
    inverse_problem.jl     # InverseProblemSetup
    global_state.jl        # GlobalState, global problem state management
  elements/
    isotropic_incompressible.jl   # Model 1 (85×85)
    orthotropic.jl                # Model 2 (85×85)
    isotropic_compressible.jl     # Model 3 (81×81)
    poroelastic.jl                # Model 4 (Tet4, 16×16)
    transverse_isotropic_v1.jl    # Model 5 (85×85)
    transverse_isotropic_v2.jl    # Model 6 (85×85)
    generalized_anisotropic.jl    # Model 7 (85×85)
    material_evaluation.jl        # GaussPointMaterial, property evaluation
  forward/
    assembly.jl            # Global stiffness assembly, Dirichlet BCs
    solver.jl              # DirectSolver, MUMPSSolver, forward_solve!
    adjoint.jl             # Adjoint solve for gradient computation
  inverse/
    gradient.jl            # Adjoint gradient ∂Φ/∂θ (all 7 models)
    jacobian.jl            # Jacobian computation (diagonal Hessian)
    objective.jl           # Objective function Φ = ½‖u_calc − u_meas‖²
  optimization/
    common.jl              # ForwardProblemContext, ConvergenceHistory
    line_search.jl         # Armijo backtracking
    conjugate_gradient.jl  # Polak-Ribière-Plus CG
    gauss_newton.jl        # Diagonal-Hessian preconditioned GD
    quasi_newton.jl        # L-BFGS two-loop recursion
  regularization/
    abstract_regularization.jl  # AbstractRegularization, WeightSchedule
    tikhonov.jl                 # L2 deviation penalty
    total_variation.jl          # Edge-preserving TV
    soft_prior.jl               # Region-based Laplacian prior
    exterior_constraint.jl      # Domain boundary penalty
    joachimowicz.jl             # Error-scaled Hessian modifier
    marquardt.jl                # Normalized LM damping
    van_houten.jl               # Exponent clipping post-processor
    spatial_filter.jl           # Gaussian smoothing post-processor
    bounds.jl                   # McGarry Poisson bounds
    cg_update_scale.jl          # CG step size capping
  io/
    file_io.jl             # .nod/.elm/.dsp/.mtr/.pre/.bnd/.bcs/.basis/.meshind
    mesh_io.jl             # read_mesh_legacy, Fortran↔Ferrite permutations
    runfile_parser.jl      # RunfileConfig, parse_runfile, setup_forward_problem
    inverse_runfile.jl     # InverseRunfileConfig, parse_inverse_runfile
    inverse_setup.jl       # setup_inverse_problem, run_inverse_solve!
  preprocessing/
    hex_mesh_generation.jl # HexMeshConfig, generate_hex_mesh
    tet_mesh_generation.jl # Tet4 mesh generation (stub)
    mri_data_conversion.jl # MREData, Siemens/WashU/UIUC readers
    dti_processing.jl      # DTIData, process_dti, fiber interpolation
  postprocessing/
    visualization.jl       # export_vtk, PropertyMapData
    convergence.jl         # print_convergence_table, export_convergence_csv
    property_interpolation.jl # interpolate_properties, interpolate_to_grid
  utils/
    coordinate_transforms.jl  # fiber_rotation_matrix, bond_matrix
    material_mesh_interp.jl   # MaterialMesh, GP2Mtr, hex8_basis
    sensitivity.jl            # Parameter sensitivity analysis
    timing.jl                 # TimingData for profiling
  distributed/
    zone_types.jl          # ZoneConfig, Zone, ZoneGrid, ZoneResult, etc.
    zone_grid.jl           # Zone center placement, pruning, growing
    zone_extraction.jl     # Extract zone mesh, BCs, material from global
    zone_boundary.jl       # Zone boundary node classification, BC construction
    zone_solve.jl          # Per-zone CG/GN/L-BFGS optimization
    zone_consolidation.jl  # GP-weighted material averaging, displacement averaging
    zone_driver.jl         # Full zone decomposition iteration driver
  ext/
    SentinelMUMPSExt.jl    # MUMPS solver extension (weakdep)
    SentinelCUDAExt.jl     # CUDA extension (planned)
    SentinelMPIExt.jl      # MPI extension (planned)
    SentinelMetalExt.jl    # Metal extension (planned)

Zone-Based Domain Decomposition

Sentinel uses a zone-based domain decomposition approach (M6.1, complete) for solving large-scale inverse problems. The global domain is partitioned into overlapping subdomains ("zones"), each of which is solved independently.

Zone Placement and Grid

Zone centers are placed using a random seed algorithm that fills the domain with approximately uniformly distributed zones. The ZoneGrid manages the spatial layout, including zone pruning (removing zones with too few elements) and zone growing (expanding zones until they reach the target size). Each zone overlaps with its neighbors to ensure smooth material property transitions.

Master-Worker Pattern

The zone decomposition driver follows a master-worker pattern:

• Serial mode: Zones are solved sequentially in a single process

• Parallel mode: Zones are distributed via a configurable pmap_func (e.g., Distributed.pmap or a custom parallel map) for multi-process execution

Zone Solve Pipeline

Each global iteration follows this pipeline:

1. Zone Grid Construction (zone_grid.jl): Place zone centers, prune and grow zones to target sizes

2. Zone Extraction (zone_extraction.jl): Extract per-zone meshes, boundary conditions, and material properties from the global problem

3. Zone Boundary Classification (zone_boundary.jl): Identify boundary nodes and construct Dirichlet BCs from measured displacements

4. Zone Forward/CG Solve (zone_solve.jl): Run per-zone optimization (CG, GN, or L-BFGS) to update local material properties

5. Consolidation (zone_consolidation.jl): Merge zone results back into the global material map using Gauss-point-weighted averaging; average overlapping displacement fields

Key Types

• ZoneConfig: Configuration parameters (zone size, overlap, solver settings)

• Zone: A single subdomain with its local mesh, BCs, and material properties

• ZoneGrid: The collection of all zones with their spatial layout

• ZoneResult: Per-zone solve output (updated material, displacement, convergence)

Key Conventions

Complex Symmetric Matrices

All FEM stiffness matrices are complex symmetric: $K=K^T$ (NOT Hermitian). This arises from the viscoelastic constitutive law using complex moduli.

DOF Layout

• Displacement DOFs: managed by Ferrite's DofHandler, 3 components × 27 nodes = 81 per element

• Pressure DOFs (elemental): appended after all displacement DOFs, indexed as ndofs(dh) + 4*(el-1) + p

• Poroelastic pressure (nodal): 4 nodal pressure DOFs interleaved with displacement

Dual Mesh System

• Displacement mesh: Hex27 (or Tet4) Ferrite grid for FEM

• Material mesh: Structured hex8 grid for material properties, with GP2MTR mapping connecting Gauss points on the displacement mesh to material property elements

Package Extensions

Optional backends are implemented as Julia package extensions (weak dependencies):

Extension	Trigger	Provides
SentinelMUMPSExt	using MUMPS	MUMPSSolver for complex symmetric sparse direct
SentinelMetalExt	using Metal	Metal GPU arrays (Float32), gradient kernel backend
SentinelCUDAExt	using CUDA, CUDSS	CUDA GPU arrays (Float64), cuDSS batched solver, GPU GMRES
SentinelMPIExt	using MPI	MPI-based zone distribution (planned)

Each extension defines backend-specific methods for the generic dispatch points:

• _gpu_float_type(backend) — Float32 for Metal, Float64 for CUDA

• _to_gpu_array(backend, x) — MtlArray for Metal, CuArray for CUDA

• linear_solve!(x, K, f, solver) — per-backend solver implementations

Solver Dispatch Hierarchy

AbstractLinearSolver
├── DirectSolver          # Julia UMFPACK (K \ f)
├── CachedDirectSolver    # LU caching for forward+adjoint pairs
├── MUMPSSolver           # MUMPS via extension (sym=2 complex symmetric)
├── KrylovSolver          # GMRES + ILU(0) via Krylov.jl
├── CUDAKrylovSolver      # GPU GMRES with cuSPARSE (via extension)
├── CUDSSDirectSolver     # GPU sparse direct via cuDSS (via extension)
└── CUDSSBatchedSolver    # Batched GPU factor+solve for all zones

See Solver Backends for details on when to use each.

GPU Gradient Architecture

The gradient computation is the primary GPU-accelerated component, using KernelAbstractions.jl for portable kernels that run on CPU, Metal, and CUDA:

AbstractGradientBackend
├── CPUGradientBackend    # Default, standard CPU gradient
└── KAGradientBackend     # GPU kernel via KernelAbstractions
    ├── CPU()             # KA on CPU (for testing)
    ├── MetalBackend()    # Apple Silicon GPU (Float32)
    └── CUDABackend()     # NVIDIA GPU (Float64)

The batched GPU driver (zone_decomposition_solve_batched_gpu!) uses a bulk-synchronous parallel (BSP) execution model where all zones synchronize at phase boundaries, enabling a single GPU kernel launch for the gradient computation across all ~290 zones simultaneously.

See GPU Acceleration Overview for the full architecture.