Forward Solver

The forward solver assembles the global stiffness matrix, applies boundary conditions, and solves the linear system $\mathbf{K}\mathbf{u}=\mathbf{f}$.

DOF Handler Setup

Sentinel.setup_hex27_dofhandler Function

setup_hex27_dofhandler(grid) -> DofHandler

Create a Ferrite DofHandler with a 3-component vector displacement field on a hex27 mesh. Each of the 27 nodes per element gets 3 DOFs (u_x, u_y, u_z).

Returns a closed DofHandler ready for DOF queries and assembly iteration.

source

Stiffness Assembly

The global system $\mathbf{K}\mathbf{u}=\mathbf{f}$ is assembled from element contributions. For mixed formulations (Models 1, 2, 5, 6, 7), the stiffness matrix includes both displacement and pressure DOFs.

Sentinel.allocate_stiffness Function

allocate_stiffness(dh::DofHandler, model::AbstractMaterialModel;
                   condense_pressure::Bool=false) -> SparseMatrixCSC{ComplexF64}

Preallocate the global complex stiffness matrix with the correct sparsity pattern for the given DofHandler and material model.

For mixed (incompressible) models (1, 2, 5, 6, 7), elemental pressure DOFs are appended after the Ferrite displacement DOFs. The global index for element el and local pressure DOF p (1-based) is ndofs(dh) + 4*(el-1) + p.

For pure-displacement models (IsotropicCompressible), uses Ferrite's native allocate_matrix for optimal sparsity structure.

With condense_pressure=true, the elemental pressure DOFs of a mixed model are statically condensed out at assembly time (assemble_stiffness! detects the displacement-only size), so the allocated matrix is the pure-displacement pattern — the condensed Schur complement adds coupling only among each element's own displacement DOFs, which the pattern already contains.

Returned matrix size

Model	Size
Incompressible (1,2,5,6,7)	(ndofs(dh)+4*ne) × same
Incompressible, condensed	ndofs(dh) × ndofs(dh)
Compressible (3)	ndofs(dh) × ndofs(dh)

source

Sentinel.assemble_stiffness! Function

assemble_stiffness!(K, dh, cv_disp, cv_press, model, props::GaussPointMaterial, ω) -> K

Assemble the global complex stiffness matrix using uniform material properties (the same GaussPointMaterial at every Gauss point of every element).

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} from allocate_stiffness(dh, model)

• dh: DofHandler with vector field Lagrange{RefHexahedron, 2}()^3

• cv_disp: CellValues for scalar Lagrange{RefHexahedron, 2}(), used for shape function evaluation inside the element kernels

• cv_press: CellValues for pressure field (or nothing; currently unused)

• model: AbstractMaterialModel singleton for kernel dispatch

• props: GaussPointMaterial applied uniformly at all Gauss points

• ω: angular frequency rad/s

Notes

• K is zeroed on entry; existing values are discarded.

• For incompressible models, pressure DOFs are at global indices ndofs(dh) + 4*(el-1) + p for element el and local DOF p (1..4).

• Static condensation: if K has the displacement-only size ndofs(dh) for a mixed model (allocated with condense_pressure=true), the elemental pressure DOFs are condensed out per element (_condense_pressure!) before scatter — the assembled K is then the exact Schur complement.

• The assembled matrix is complex symmetric (K = Kᵀ), not Hermitian.

Modifies

K in-place; also returns K.

source

assemble_stiffness!(K, dh, cv_disp, cv_press, model, material, gp2mtrs, meshes, ω;
                    ρ=1000.0, fiber_field=nothing) -> K

Assemble the global complex stiffness matrix from a Material struct with spatially varying properties interpolated via the GP2MTR mapping.

For each element and each Gauss point:

1. Maps the Ferrite flat quadrature index to the GP2Mtr (ig, jg, kg) indices

2. Evaluates each material property at the Gauss point location using the material mesh interpolation (evaluate_property_at_gp)

3. Constructs a GaussPointMaterial via model-specific build_gp_material

4. Calls element_stiffness! with the per-GP materials

5. Scatters into the global K

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} from allocate_stiffness

• dh: DofHandler with vector field Lagrange{RefHexahedron, 2}()^3

• cv_disp: CellValues for scalar Lagrange{RefHexahedron, 2}()

• cv_press: CellValues for pressure field (or nothing)

• model: AbstractMaterialModel singleton for kernel dispatch

• material: Material struct with per-property SingleProperty arrays

• gp2mtrs: vector of GP2Mtr (one per material mesh)

• meshes: vector of MaterialMesh (one per material mesh)

• ω: angular frequency rad/s

Keyword Arguments

• ρ: density kg/m³ (default 1000.0)

• fiber_field: optional per-element fiber direction vectors (for TI models)

Notes

• K is zeroed on entry; existing values are discarded.

• Static condensation: if K has the displacement-only size ndofs(dh) for a mixed model (allocated with condense_pressure=true), the elemental pressure DOFs are condensed out per element (_condense_pressure!) before scatter — the assembled K is then the exact Schur complement.

• Currently requires prop.rmesh == prop.imesh for all properties (i.e., real and imaginary parts on the same material mesh).

source

Boundary Conditions

Dirichlet boundary conditions are applied via symmetric elimination: the constrained row and column are zeroed, and the diagonal is set to 1.

Sentinel.apply_dirichlet! Function

apply_dirichlet!(K, f, constrained_dofs, prescribed_values)

Apply Dirichlet boundary conditions by modifying K and f in-place.

Uses the symmetric elimination method:

1. For each constrained DOF k with prescribed value v:

• Update f[i] -= K[i,k] * v for all rows i (read column k before zeroing)

• Zero column k of K (efficient in CSC format)

• Zero row k of K (scan all columns)

• Set K[k,k] = 1, f[k] = v

This preserves complex symmetry of K (K = Kᵀ) after BC application.

Arguments

• K: preallocated SparseMatrixCSC{ComplexF64} to modify in-place

• f: right-hand side vector, length size(K, 1)

• constrained_dofs: 1-based global DOF indices to constrain

• prescribed_values: corresponding prescribed complex values

source

apply_dirichlet!(K, f, bcs::BoundaryConditions, dispset::Int)

Apply Dirichlet boundary conditions from a BoundaryConditions struct for the given displacement set index.

Iterates bcs.bcind[:, dispset] to find all constrained DOFs with BC_DIRICHLET type, then delegates to the DOF-vector overload.

source

Linear Solvers

Sentinel.AbstractLinearSolver Type

AbstractLinearSolver

Abstract type for sparse linear solver backends. Subtypes must implement linear_solve!(x, K, f, solver).

source

Sentinel.DirectSolver Type

DirectSolver <: AbstractLinearSolver

Julia's built-in sparse direct solver (SuiteSparse/UMFPACK). Works for any sparse complex system without additional dependencies.

source

Sentinel.CachedDirectSolver Type

CachedDirectSolver <: AbstractLinearSolver

UMFPACK solver that caches the LU factorization between calls. When K has not changed since the last solve, the cached factorization is reused for subsequent right-hand sides, eliminating repeated factorization.

Call invalidate_cache!(solver) whenever K changes (e.g. after re-assembly).

This is particularly useful for paired forward + adjoint solves on the same K: the forward solve stores lu(K), and the adjoint reuses it.

source

Sentinel.CPUBatchedSolver Type

CPUBatchedSolver <: AbstractLinearSolver

CPU implementation of the batched-solver interface (create_batched_solver!, batched_update_and_refactor!, batched_solve!, similar_batched_solver, batched_solver_ready). Holds one cached UMFPACK LU factorization per zone and operates on the whole batch with Threads.@threads :greedy — the same work-stealing schedule the per-zone driver loops use, since per-zone factor/solve cost is uneven.

This is the Apple-Silicon (and any non-CUDA) analogue of CUDSSBatchedSolver: passing one to zone_decomposition_solve_batched_gpu! as the batched solver switches the whole driver onto the batched path — batched forward/adjoint, the single-dispatch batched (Metal/CPU) gradient, and the lockstep batched secant line search — with the linear solves on CPU instead of the GPU.

The line-search payoff over the per-zone secant is symbolic-analysis reuse: the per-zone path re-runs a full lu (symbolic + numeric) on every trial, whereas batched_update_and_refactor! keeps the symbolic analysis and recomputes only the numeric factors via lu!(F, K). The zone set — and hence each zone's sparsity — is fixed within a CG iteration, which is exactly when refactor is used.

Line-search "lag" mode (lag)

Profiling shows the numeric refactorization dominates the Mac solve cost, and the line search refactors every trial. But trial materials are close to the CG point, so the CG-point factorization M = LU(K_cgpoint) is a near-exact preconditioner for each trial system K_trial·x = b. When lag is set (the driver sets it around the batched line search only if lag_ls_enabled), batched_update_and_refactor! does not factor — it just captures the trial K — and batched_solve! solves by Richardson iterative refinement x += M⁻¹(b − K_trial·x): a few cheap (SpMV + triangular-solve) steps instead of a full factorization. Any zone whose IR does not reach ir_tol within ir_maxit falls back to a real lu, so a lagged solve is always as accurate as the direct one (the line-search decisions, and hence the reconstruction, are unchanged). n_lagged/n_fallback count the split for diagnostics.

Lagged IR is OFF by default — a measured regression

On the real MGH problem the secant line search takes large steps (the α=1 probe especially), so trial materials are far from the CG point: ~82% of zone solves failed to converge and fell back to a full factorization (after paying for the wasted IR sweeps), making the run ~38% slower than direct refactor. The flag is retained (opt in with CPUBatchedSolver(lag_ls=true)) because the SpMV + IR scaffold is reusable and a smarter preconditioner-update policy might yet pay off, but it must not be the default.

source

Sentinel.MUMPSSolver Type

MUMPSSolver <: AbstractLinearSolver

MUMPS sparse direct solver. Requires the MUMPS.jl package extension.

Fields

• sym::Int: MUMPS symmetry flag. 2 for complex symmetric (default, used for viscoelastic models 1,2,3,5,6,7), 0 for general unsymmetric (poroelastic model 4).

source

Sentinel.linear_solve! Function

linear_solve!(x, K, f, solver::AbstractLinearSolver) -> x

Solve the sparse linear system K·x = f in-place, storing the result in x.

Arguments

• x::Vector{ComplexF64}: preallocated solution vector (overwritten)

• K::SparseMatrixCSC{ComplexF64}: assembled stiffness matrix (with BCs applied)

• f::Vector{ComplexF64}: right-hand side vector (with BCs applied)

• solver: solver backend

source

Forward Solve

High-level forward solve that handles DOF scattering:

Sentinel.forward_solve! Function

forward_solve!(disp::Displacement, K, f, model::AbstractMaterialModel,
               dispset::Int; solver=DirectSolver()) -> Displacement

Solve the forward problem K·u = f and store the solution in disp.

K and f must already have boundary conditions applied (via apply_dirichlet!). The solution vector is scattered into disp for the given dispset using scatter_solution!.

Arguments

• disp: preallocated Displacement struct (from init_displacement!)

• K: assembled stiffness matrix with BCs applied

• f: right-hand side vector with BCs applied

• model: material model (used for dispatch/documentation only)

• dispset: displacement set index (1-based)

• solver: linear solver backend (default: DirectSolver())

source

Adjoint Solve

The adjoint problem $\bar{\mathbf{K}}\mathit{\lambda }=-({\mathbf{u}}_{\mathrm{calc}}-{\mathbf{u}}_{\mathrm{meas}})$ is solved for gradient computation. Since $\mathbf{K}$ is complex symmetric, $\bar{\mathbf{K}}={\mathbf{K}}^H$.

Sentinel.adjoint_solve! Function

adjoint_solve!(adjnt::Displacement, K, calc::Displacement, meas::Displacement,
               bcs::BoundaryConditions, dispset::Int;
               solver::AbstractLinearSolver = DirectSolver()) -> Displacement

Solve the adjoint problem conj(K)·λ = -(u_calc - u_meas) and store the result in adjnt for the given displacement set.

K should be the stiffness matrix with forward Dirichlet BCs already applied (rows/columns zeroed, diagonal = 1 at constrained DOFs). The adjoint uses homogeneous BCs on the same DOFs by zeroing the RHS there.

Arguments

• adjnt: preallocated Displacement struct (from init_displacement!)

• K: assembled stiffness matrix with forward BCs applied

• calc: calculated displacement (from forward solve)

• meas: measured displacement (experimental data)

• bcs: BoundaryConditions identifying constrained DOFs

• dispset: displacement set index (1-based)

• solver: linear solver backend (default: DirectSolver())

source