Inversion module (inversion.h)

Typedefs

typedef PetscErrorCode (*InversionObjGradFn)(Vec X, PetscReal *F, Vec G, void *ctx)

Objective + gradient callback type for L-BFGS.

Signature: (Vec X, PetscReal *F, Vec G, void *ctx). Must be called collectively on all MPI ranks.

Enums

enum ObservedDataMode

Observed-data source mode for the inverse kernel.

Selects which schema loadObservedDataset() reads the [numFreqs x numReceivers] Ex misfit data from:

OBS_EXTERNAL : the bundle’s /observed/Ex compound-complex dataset (the post-noise external observed data; default).
OBS_FM_NATIVE: an fm.csem responses HDF5 file’s native per-source /sources/src{k}/fields/Ex Vecs (source k -> frequency row k). Lets a forward run feed the inverse kernel directly, without the Python reshape step (noise, if wanted, must already be baked into that file).

Values:

enumerator OBS_EXTERNAL

enumerator OBS_FM_NATIVE

Functions

PetscErrorCode readInversionParams(imParams *iparams)

Reads inversion parameters from the PETSc options database.

This function initializes inversion-control parameters and applies optional CLI overrides for optimization, regularization, diagnostics, and debugging behavior.

The FEM basis order (-nord) is optional. The bundled HDF5 input is the primary source of truth, and its value is later propagated into the inversion parameters after loadCsemInputs. A CLI-provided -nord overrides the bundled value. A value of 0 means “not overridden”.

The routine also initializes defaults for L-BFGS settings, Tikhonov regularization, RMS stopping criteria, snapshot generation, and optional smoothing controls.

Transmiter/frequency metadata is not loaded here; it is populated later from the unified input bundle.

Parameters:

iparams – [out] Struct receiving the parsed inversion parameters.
im_Params – [out] Structure containing inversion parameters.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode loadInversionMetaFromBundle(const char *bundleFile, imParams *iparams)

Applies case-property defaults from the bundle to iparams.

Reads the error_level attribute on /observed and the fixed_materials array under /inv_meta and applies them UNLESS the corresponding CLI override was present (see the *FromCLI provenance flags). Safe to call even when the bundle has no such entries - iparams keeps the readInversionParams defaults.

Applies case-property defaults from the bundle to iparams.

This function reads optional inversion-related metadata from the unified PETGEM input file and applies it to the inversion parameter structure. It complements readInversionParams() by filling values that are defined at the case level rather than the CLI.

The following fields may be updated from the bundle:

/observed@error_level → im_Params->errorLevel
/inv_meta/fixed_materials → fixed material IDs list

CLI values always take precedence:

If errorLevelFromCLI is set, /observed is ignored.
If fixedMaterialsFromCLI is set, bundle values are ignored.

Missing entries are silently skipped, leaving defaults unchanged.

Parameters:

bundleFile – [in] Path to the unified PETGEM HDF5 bundle.
iparams – [inout] Inversion parameters updated in place.
bundleFile – [in] Path to PETGEM HDF5 bundle.
im_Params – [inout] Inversion parameter structure to update.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode setupInversionSources(const char *bundleFile, imParams *iparams)

Loads multi-frequency inversion sources from the unified bundle.

Reads the bundle’s /sources group (replacing the legacy text-file format) and populates iparams->numFreqs and invSources[]. bundleFile is the same HDF5 path consumed by loadCsemInputs.

Loads multi-frequency inversion sources from the unified bundle.

This function reads the /sources group from the input bundle and constructs the inversion source list used by the optimization layer. Each entry corresponds to one frequency–dipole configuration, with shared structure between forward and inverse problems.

The data is stored as PETSc Vecs, so loading is performed using PetscViewerHDF5 + VecLoad to preserve correct scalar-type handling.

The function populates:

im_Params->invSources[] (per-frequency dipole definitions)
im_Params->numFreqs (number of loaded entries)

A hard limit (INV_MAX_FREQUENCIES) is enforced to prevent oversized inversion problems.

Parameters:

bundleFile – [in] Path to the unified PETGEM HDF5 bundle.
iparams – [inout] Inversion parameters whose sources are filled.
bundleFile – [in] Path to unified PETGEM HDF5 input file.
im_Params – [out] Inversion parameter structure to populate.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode loadObservedData(const char *bundleFile, PetscInt numFreqs, PetscInt numReceivers, Mat *dObs)

Loads observed data from the unified bundle’s /observed/Ex dataset.

Reads the HDF5 compound complex128 dataset of shape [numFreqs, numReceivers] into a dense Mat on PETSC_COMM_SELF. This is the OBS_EXTERNAL backend of loadObservedDataset().

Loads observed data from the unified bundle’s /observed/Ex dataset.

This function reads the /observed/Ex dataset from the unified HDF5 input file and constructs a dense observation matrix of shape [numFreqs × numReceivers].

The dataset is stored in HDF5 as a complex compound type {real, imag} (complex128). Rank 0 performs the file I/O and then broadcasts the data to all MPI ranks.

The result is stored as a PETSc dense matrix (Mat) on PETSC_COMM_SELF for local use on each process.

Parameters:

bundleFile – [in] Path to the unified PETGEM HDF5 bundle.
numFreqs – [in] Number of frequencies (rows).
numReceivers – [in] Number of receivers (columns).
dObs – [out] Dense Mat (numFreqs × numReceivers) of observed Ex.
bundleFile – [in] Path to unified PETGEM HDF5 file.
numFreqs – [in] Number of frequencies (rows).
numReceivers – [in] Number of receivers (columns).
dObs – [out] Dense observation matrix.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode loadObservedFmNative(const char *responsesFile, PetscInt numFreqs, PetscInt numReceivers, Mat *dObs)

Loads observed Ex from an fm.csem responses file’s native schema.

OBS_FM_NATIVE backend of loadObservedDataset(): reads the per-source /sources/src{k}/fields/Ex PETSc Vecs written by the forward kernel (k = 1..numFreqs maps to frequency rows 0..numFreqs-1) into the same [numFreqs × numReceivers] dense Mat layout the external path produces, so the inverse kernel is agnostic to the data origin. Each rank reads the COMM_SELF Vecs independently (no broadcast).

Loads observed Ex from an fm.csem responses file’s native schema.

Counterpart to loadObservedData() that consumes the forward kernel’s native output schema (/sources/src{k}/fields/Ex PETSc Vecs) instead of the pre-reshaped /observed/Ex dataset. Source k (1-based) maps to frequency row k-1 of the returned dense Mat. Every rank opens the COMM_SELF viewer and VecLoads each row independently, so the resulting PETSC_COMM_SELF Mat is identical on all ranks (mirroring loadObservedData’s post-broadcast state).

See the inversion.h prototype for parameter semantics.

Parameters:

responsesFile – [in] Path to the fm.csem responses HDF5 file.
numFreqs – [in] Number of source/frequency rows.
numReceivers – [in] Number of receivers (Ex length / columns).
dObs – [out] Dense Mat (numFreqs × numReceivers) of observed Ex.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

PetscErrorCode loadObservedDataset(const imParams *iparams, PetscInt numReceivers, Mat *dObs)

Observed-data abstraction entry point used by the inverse kernel.

Dispatches to loadObservedData() (OBS_EXTERNAL) or loadObservedFmNative() (OBS_FM_NATIVE) according to iparams->observedMode, resolving the file path from iparams->observedFile (falling back to the unified bundle iparams->fm.inputFile when empty). numFreqs is taken from iparams->numFreqs.

Observed-data abstraction entry point used by the inverse kernel.

Resolves the file path (iparams->observedFile, or the unified bundle iparams->fm.inputFile when empty) and dispatches to the backend selected by iparams->observedMode. Keeps the inverse driver agnostic to data origin.

Parameters:

iparams – [in] Inversion parameters (mode, file, numFreqs, bundle).
numReceivers – [in] Number of receivers (columns).
dObs – [out] Dense Mat (numFreqs × numReceivers) of observed Ex.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

PetscErrorCode buildNeighborSmoothingGraph(const DM dm, const Grid *grid, const imParams *iparams, Vec materialsID, NeighborGraph *graph)

Builds the CSR neighbor smoothing graph from DMPlex topology.

Cells whose material_id matches any entry in iparams->fixedMaterials are treated as fixed (self-referencing only, excluded from smoothing).

For each local cell i: finds all cells sharing at least one vertex via the DMPlex star; excludes fixed cells (material_id in iparams->fixedMaterials); weights w_ij = 1/dist(centroid_i, centroid_j), normalized. Fixed cells with no valid neighbors keep a self-reference.

Parameters:

dm – [in] DMPlex mesh.
grid – [in] Finite-element grid descriptor.
iparams – [in] Inversion parameters (fixed-material list).
materialsID – [in] Per-cell material-id Vec.
graph – [out] Neighbor graph to populate.
dm – [in] DMPlex mesh.
grid – [in] Finite-element grid descriptor.
iparams – [in] Inversion parameters (fixed-material list).
materialsID – [in] Per-cell material-id Vec.
graph – [out] Neighbor graph populated with CSR data and flags.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode computeGradientContribution(const DM dm, const Grid *grid, const Vec conductivity, const Vec xLocal, const Vec nxLocal, PetscScalar constFactor, DM dmInversion, Vec DfDm, const Quadrature3D *quadrature_3d, PetscReal **Me, PetscReal **Ke)

Accumulates the per-element adjoint gradient into DfDm (1 DOF/cell).

Computes DfDm[ie] += real( (-2·constFactor·Me_e·x_e)^T · nx_e ) using a plain transpose (no conjugate), matching MATLAB iG.’*inx. quadrature_3d, Me, Ke are caller-owned workspace buffers (set up once on InversionContext via setupInversionWorkspace).

Accumulates the per-element adjoint gradient into DfDm (1 DOF/cell).

For each local cell ie: idKdm = -2 · constFactor · Me_e, iG = idKdm · x_e, DfDm[ie] += real( iG · nx_e ) (plain transpose, MATLAB iG.’*inx). The 3D quadrature and Me/Ke buffers are provided by the caller (allocated once on the InversionContext via setupInversionWorkspace); the hoist is a pure lifetime change - the values produced by each cell evaluation are byte-identical to the per-callback version.

Parameters:

dm – [in] H(curl) DM.
grid – [in] Finite-element grid descriptor.
conductivity – [in] Current conductivity Vec.
xLocal – [in] Forward solution (local).
nxLocal – [in] Adjoint solution (local).
constFactor – [in] Frequency factor iωμ.
dmInversion – [in] DM for the inversion field (1 DOF/cell).
DfDm – [inout] Gradient accumulator (local, 1 DOF/cell).
quadrature_3d – [in] 3D quadrature workspace.
Me – [inout] Scratch elemental mass buffer.
Ke – [inout] Scratch elemental stiffness buffer.
dm – [in] H(curl) DM.
grid – [in] Finite-element grid descriptor.
conductivity – [in] Current conductivity Vec.
xLocal – [in] Forward solution (local).
nxLocal – [in] Adjoint solution (local).
constFactor – [in] Frequency factor iωμ.
dmInversion – [in] DM for the inversion field (1 DOF/cell).
DfDm – [inout] Gradient accumulator (local, 1 DOF/cell).
quadrature_3d – [in] 3D quadrature workspace.
Me – [inout] Scratch elemental mass buffer.
Ke – [inout] Scratch elemental stiffness buffer.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode applyGaussSeidelSmoothing(const NeighborGraph *graph, PetscReal diagWeight, Vec v)

Applies forward + reverse Gauss-Seidel smoothing to a vector.

On serial / single-rank runs this gathers to rank 0, sweeps, and scatters back; the parallel block-Jacobi path is used when the graph carries the overlap-1 state (see NeighborGraph::hasParallelGraph).

Applies forward + reverse Gauss-Seidel smoothing to a vector.

Each updated value is computed from a frozen snapshot of the array (not in place), so the result is independent of the visiting order. This is the key property that makes the smoother bit-reproducible across any MPI partition: an in-place Gauss-Seidel sweep depends on the (partition-derived) cell order and on stale ghost values, which at high rank counts under-regularizes the many partition boundaries and drives runaway overfitting there (observed ρ up to ~1e6 at 336 ranks). Jacobi removes that order dependence. The legacy name is retained for call-site stability. Operates on a partition-local Vec; non-owned DOFs are not visited.

Parameters:

graph – [in] Neighbor smoothing graph.
diagWeight – [in] Self-weight applied to the diagonal during sweeps.
v – [inout] Vector smoothed in place.
graph – [in] Neighbor smoothing graph.
diagWeight – [in] Self-weight applied to the diagonal during sweeps.
v – [inout] Vector smoothed in place.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode applyLogToSigma(DM dmInversion, DM dmConductivity, const Vec X, const Vec X0, Vec sigmaModel, const Grid *grid, const NeighborGraph *graph, PetscReal diagWeight, Vec xPostSmoothOut)

Recovers conductivity from a log-conductivity perturbation.

Computes sigma = 1 / exp(X_smooth + X0) (isotropic, written to components 0-2). X and X0 are global Vecs on dmInversion (1 DOF/cell); sigmaModel is a local Vec on dmConductivity. Gauss-Seidel smoothing is applied to a local copy of X before exp, matching MATLAB’s tempX smoothing after each L-BFGS step.

Recovers conductivity from a log-conductivity perturbation.

X and X0 are global Vecs on dmInversion (1 DOF/cell); sigmaModel is a local Vec on dmConductivity (3 DOF/cell), with components 0-2 (res_x, res_y, res_z) set to σ. The smoother (graph/diagWeight) is applied to a local copy of X before exp, matching MATLAB’s tempX smoothing after each L-BFGS step; X itself is NOT modified.

Parameters:

dmInversion – [in] DM for the inversion field (1 DOF/cell).
dmConductivity – [in] DM for the conductivity field.
X – [in] Log-perturbation iterate (global).
X0 – [in] Initial log(rho) (global).
sigmaModel – [out] Recovered conductivity (local).
grid – [in] Finite-element grid descriptor.
graph – [in] Neighbor smoothing graph.
diagWeight – [in] Smoother self-weight.
xPostSmoothOut – [out] Optional (may be NULL): local Vec on dmInversion receiving the smoothed X (pre-X0 add) for VTU diagnostic snapshots.
dmInversion – [in] DM for the inversion field (1 DOF/cell).
dmConductivity – [in] DM for the conductivity field.
X – [in] Log-perturbation iterate (global).
X0 – [in] Initial log(ρ) (global).
sigmaModel – [out] Recovered conductivity (local).
grid – [in] Finite-element grid descriptor.
graph – [in] Neighbor smoothing graph.
diagWeight – [in] Smoother self-weight.
xPostSmoothOut – [out] Optional: local Vec on dmInversion receiving the smoothed X (pre-X0 add) for VTU snapshots.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode createInversionDM(DM dmConductivity, const Grid *grid, DM *dmInv)

Creates a DM with 1 DOF per cell for the inversion variables.

Clones the DMPlex topology from dmConductivity but sets a new section with a single scalar DOF per cell (for X, X0, gradient).

Clones the DMPlex topology from dmConductivity and installs a new PetscSection with a single scalar DOF on each cell. A 1-component PetscFV is registered as field 0 so cell-centered data can be written to VTU (DMPlexVTKWriteAll_VTU traverses registered fields; without one, DMGetField fails).

Parameters:

dmConductivity – [in] DM whose topology is cloned.
grid – [in] Finite-element grid descriptor.
dmInv – [out] New 1-DOF/cell inversion DM.
dmConductivity – [in] DM whose topology is cloned.
grid – [in] Finite-element grid descriptor.
dmInv – [out] New 1-DOF/cell inversion DM.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode createInvKSP(const imParams *iparams, const DM dm, const Mat A, const Mat Gbddc, KSP *ksp)

Creates a KSP for inversion (same setup as solveCsemSystem, no solve).

Gbddc is the forward-formulation topological discrete-gradient operator passed to PCBDDC (distinct from the inversion gradient ∂F/∂X). The caller must KSPDestroy it after the forward + adjoint solves.

Creates a KSP for inversion (same setup as solveCsemSystem, no solve).

The caller invokes solveInvSystem twice (forward + adjoint) then KSPDestroy. The MUMPS factorization is triggered on the first KSPSolve.

Parameters:

iparams – [in] Inversion parameters.
dm – [in] H(curl) DM.
A – [in] System matrix.
Gbddc – [in] Topological discrete-gradient operator for PCBDDC.
ksp – [out] Created KSP bound to A.
iparams – [in] Inversion parameters. Reserved; currently unused (solver options come from KSPSetFromOptions).
dm – [in] H(curl) DM.
A – [in] System matrix.
Gbddc – [in] Topological discrete-gradient operator for PCBDDC.
ksp – [out] Created KSP bound to A.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode solveInvSystem(const KSP ksp, const Vec rhs, Vec sol)

Solves A·sol = rhs with an already-created (factored) KSP.

Parameters:

ksp – [in] KSP previously created by createInvKSP.
rhs – [in] Right-hand side vector.
sol – [out] Solution vector.
ksp – [in] KSP previously created by createInvKSP.
rhs – [in] Right-hand side vector.
sol – [out] Solution vector.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode inversionObjGrad(Vec X, PetscReal *F, Vec G, void *ctx)

Objective/gradient callback implementation for the inverse kernel.

Per call:

Recover σ from X (applyLogToSigma).
Reassemble the LHS using the updated σ.
Frequency loop: forward solve → field interpolation → residual → adjoint RHS → adjoint solve → gradient accumulation.
Chain rule + smoothing + Tikhonov regularization.

Parameters:

X – [in] Current log-perturbation iterate.
F – [out] Objective value at X.
G – [out] Gradient at X.
ctx – [in] InversionContext pointer.
X – [in] Current log-perturbation iterate.
F – [out] Objective value at X.
Gvec – [out] Gradient at X.
ctx – [in] InversionContext pointer (cast from void*).

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode lbfgsOptimize(InversionObjGradFn objgrad, void *ctx, Vec X, PetscInt M, PetscInt maxIter, PetscReal gtol, const PetscReal *rmsPtr, PetscReal rmsTol, PetscInt *numIters, const char **reasonStr)

L-BFGS optimizer (Nocedal 1980, two-loop recursion).

Replaces PETSc TAO, which is unavailable with complex scalars. If rmsPtr is non-NULL and rmsTol > 0, the optimizer exits as soon as *rmsPtr <= rmsTol after an accepted step - matches MATLAB’s rms <= 1.05 exit in Ex_inv.m. The callback is expected to update *rmsPtr before returning.

L-BFGS optimizer (Nocedal 1980, two-loop recursion).

Implements the two-loop recursion for H·g approximation with a backtracking Armijo line search. Designed to work with complex PetscScalar (TAO is unavailable when PETSC_USE_COMPLEX is set). All optimization variables are stored with zero imaginary part, so VecDot/VecAXPY on them reduce to standard real operations.

Two early-stop modes are supported when the caller maintains an RMS value via rmsPtr:

Absolute: stop when *rmsPtr ≤ rmsTol (matches MATLAB’s rms <= 1.05 exit in Ex_inv.m).
Adaptive plateau: stop when the relative drop is below -inv_rms_rtol for -inv_rms_stall_window consecutive iterations (default 1e-3 / 3 iters).

References: Nocedal & Wright, “Numerical Optimization”, Ch. 7; MATLAB Fortran reference at petgem_inv_new/lbfgs_matlab/.

Parameters:

objgrad – [in] Objective/gradient callback.
ctx – [in] Opaque context passed to objgrad.
X – [inout] Initial iterate; final iterate on return.
M – [in] L-BFGS memory size.
maxIter – [in] Maximum number of iterations.
gtol – [in] Gradient-norm convergence tolerance.
rmsPtr – [in] Optional pointer to a caller-updated RMS value.
rmsTol – [in] RMS early-stop threshold (≤0 disables).
numIters – [out] Number of iterations performed.
reasonStr – [out] Human-readable convergence/stop reason.
objgrad – [in] Objective/gradient callback.
ctx – [in] Opaque context passed to objgrad (an InversionContext * in current usage).
X – [inout] Initial iterate; final iterate on return.
M – [in] L-BFGS memory size.
maxIter – [in] Maximum number of iterations.
gtol – [in] Gradient-norm convergence tolerance.
rmsPtr – [in] Optional pointer to a caller-updated RMS value.
rmsTol – [in] RMS early-stop threshold (≤0 disables).
numIters – [out] Number of iterations performed.
reasonStr – [out] Human-readable convergence/stop reason.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode writeInversionSnapshotVTU(const InversionContext *ctx, PetscInt acceptedIter)

Writes a VTU snapshot of the current conductivity model (ρ = 1/σ).

Called after accepted L-BFGS steps when snapshotInterval > 0; the output file is {output_dir}/inv_model_iter{N:05d}.vtu.

Writes a VTU snapshot of the current conductivity model (ρ = 1/σ).

This function exports the current inversion iterate as a ParaView compatible parallel VTK dataset (.pvtu + per-rank .vtu pieces). It is typically called every accepted L-BFGS step when snapshotting is enabled.

Each MPI rank writes only its locally owned cells (no global gather). Cell data is written in exploded tetrahedral form to avoid global point renumbering.

The snapshot carries a single cell-centered scalar field:

rho_ohm_m : 1 / sigma_x (conductivity-derived resistivity)

Output structure:

inv_model_iterXXXXX.pvtu (master file, rank 0)
inv_model_iterXXXXX_pRRRR.vtu (one file per rank)

The .pvtu file aggregates all per-rank pieces into a single dataset for visualization in ParaView.

Parameters:

ctx – [in] Inversion context.
acceptedIter – [in] Index of the accepted iteration (used in filename).
ctx – [in] Inversion execution context (DM, fields, grid).
acceptedIter – [in] Current accepted optimization iteration.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode runCsemInversion(const imParams *iparams, const DM dm, const Grid *grid, Vec conductivity, Vec materialsID, Vec receivers, PetscLogDouble *tAssemblyOut, PetscLogDouble *tSolverOut)

Top-level inversion driver.

receivers is the serial Vec (PETSC_COMM_SELF, length 3·N_recv) loaded by loadCsemInputs from /receivers in the unified PETGEM input HDF5. It is consumed once by buildReceiverInterpolationMatrices; the caller retains ownership for cleanup.

Builds the inversion DM, neighbor smoothing graph, workspace, and KSPs; runs the custom L-BFGS optimizer over the per-frequency objective/gradient callback; writes the final HDF5 results; and returns accumulated assembly and solver wall-clock times for the kernel timer report. receivers is the serial Vec produced by loadCsemInputs; it is consumed by buildReceiverInterpolationMatrices and the caller retains ownership.

Parameters:

iparams – [in] Inversion parameters.
dm – [in] H(curl) DM.
grid – [in] Finite-element grid descriptor.
conductivity – [in] Initial per-cell conductivity Vec.
materialsID – [in] Per-cell material-id Vec.
receivers – [in] Serial Vec of 3·N_recv receiver coordinates.
tAssemblyOut – [out] Accumulated assembly time (seconds) for reporting.
tSolverOut – [out] Accumulated solver time (seconds) for reporting.
iparams – [in] Inversion parameters.
dm – [in] H(curl) DM.
grid – [in] Finite-element grid descriptor.
conductivity – [in] Initial per-cell conductivity Vec.
materialsID – [in] Per-cell material-id Vec.
receivers – [in] Serial Vec of 3·N_recv receiver coordinates.
tAssemblyOut – [out] Accumulated assembly time (seconds).
tSolverOut – [out] Accumulated solver time (seconds).

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

PetscErrorCode destroyNeighborGraph(NeighborGraph *graph)

Frees the NeighborGraph memory.

Parameters:

graph – [inout] Neighbor graph to destroy.
graph – [inout] Neighbor graph to destroy.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PETSc error code otherwise.

struct InvCsemSource

#include <inversion.h>

Inversion source: a transmitter carrying its own frequency.

Extends the forward CsemSource layout with a per-source frequency so a multi-frequency inversion can store one record per (source, frequency).

Public Members

PetscReal freq: Frequency (Hz).

PetscReal position[3]: Transmitter position (x, y, z).

PetscReal current: Electric current.

PetscReal length: Dipole length.

PetscReal dipAngle: Dip angle.

PetscReal azimuthAngle: Azimuth angle.

struct imParams

#include <inversion.h>

Inversion parameters (read from the PETSc options database).

Unnamed Group

PetscBool errorLevelFromCLI: Provenance flags: PETSC_TRUE iff the field was set on the CLI (and so should NOT be overridden by the bundle reader). Error level came from CLI.

PetscBool fixedMaterialsFromCLI: Fixed materials came from CLI.

Public Members

fmParams fm: Shared base parameters, parsed by the SAME readfmParams() the forward kernel uses (input/output paths, basis order, MPI task count, quiet flag). Unifies the fm.csem / im.csem interface: the inverse-only controls below extend this common base. The basis order is fm.nord and the unified-bundle path (consumed by setupInversionSources, loadInversionMetaFromBundle, loadObservedData) is fm.inputFile - formerly the separate bundleFile mirror.

PetscInt maxIter: Max L-BFGS iterations.

PetscInt lbfgsMemory: L-BFGS M parameter.

PetscReal lambda: Tikhonov factor.

PetscReal errorLevel: Relative data error.

PetscReal gtol: Gradient convergence tolerance.

PetscReal rmsTol: RMS early-stop (<=0 off).

PetscReal diagGradientWeight: Self-weight in the smoother.

PetscInt numFixedMaterials: Fixed material IDs: cells whose material_id matches one of these values are excluded from gradient smoothing (treated as self-referencing). Defaults come from the bundle’s /inv_meta/fixed_materials dataset (written by the preprocess from sigmas.txt’s fixed column); -inv_fixed_materials on the CLI is an override. Number of fixed IDs.

PetscInt fixedMaterials[INV_MAX_FIXED_MATERIALS]: Fixed material ID list.

PetscInt numFreqs: Source-frequency entries loaded from the unified bundle’s /sources group (freq, position, current, length, dipAngle, azimuthAngle). Populated by setupInversionSources. Each entry’s frequency lives in invSources[i].freq; no separate frequency array is kept. Number of entries.

InvCsemSource invSources[INV_MAX_FREQUENCIES]: One source record per entry.

PetscInt snapshotInterval: VTU snapshot: write the conductivity model every N accepted L-BFGS steps. 0 (default) disables snapshots. Output dir is taken from -output_dir. Snapshot interval (0 = disabled).

ObservedDataMode observedMode: Observed-data abstraction: which schema/file the misfit data is read from (see ObservedDataMode). Set by readInversionParams from -inv_observed_mode. External vs fm-native source.

char observedFile[PETSC_MAX_PATH_LEN]: Path to the observed-data file. For OBS_EXTERNAL an empty string means “use the unified bundle (fm.inputFile)”; for OBS_FM_NATIVE it is the fm.csem responses HDF5 file. Set by -inv_observed_file. Observed-data file (empty = bundle).

struct NeighborGraph

#include <inversion.h>

Neighbor smoothing graph in CSR format.

Fixed elements (material_id matching any entry in fixedMaterials) are flagged so they are excluded from gradient smoothing.

Public Members

PetscInt numLocalCells: Number of locally owned cells.

PetscInt *neighborStart: CSR row pointer, size numLocalCells+1.

PetscInt *neighborList: Local cell IDs of neighbors.

PetscReal *neighborWeights: Normalized 1/dist weights.

PetscBool *isFixed: PETSC_TRUE for fixed-material cells.

PetscBool hasParallelGraph: Fully-parallel block-Jacobi GS state, set by setupParallelSmoothingGraph. When hasParallelGraph is true (production path for MPI > 1), applyGaussSeidelSmoothing runs forward + reverse Gauss-Seidel on each rank’s owned cells using its own local copy plus one layer of ghost values, exchanging ghosts between sweeps (DMLocalToGlobal / DMGlobalToLocal pair). Cost is O(local cells) per smoother call, no rank-0 bottleneck. NULL/false on single-rank runs (the local sweep in applyGaussSeidelSmoothing handles serial directly). True when the parallel GS state is built.

DM dmInversionOver: overlap=1 1-DOF/cell DM (owned).

PetscInt *oNeighborStart: CSR row pointer [numLocalCells+1].

PetscInt *oNeighborList: Local overlap-1 neighbor indices.

PetscReal *oNeighborWeights: Normalized 1/dist weights.

Vec oLocalScratch: Local Vec on dmInversionOver.

Vec oGlobalScratch: Global Vec on dmInversionOver.

struct InversionContext

#include <inversion.h>

Context passed to the L-BFGS objective/gradient callback.

Bundles the DMs, model/gradient Vecs, observed data, precomputed per-frequency inputs, and cached matrices/solvers reused across every objgrad evaluation. The pre-allocated workspace fields are owned by setupInversionWorkspace and freed by destroyInversionWorkspace.

Unnamed Group

PetscLogDouble tAssembly: Phase timers accumulated across all objgrad evaluations (seconds), so im.csem can report an Assembly/Solver breakdown consistent with fm.csem instead of lumping the whole inversion into one bucket. Ms refill + A = K - iωμ·Ms.

PetscLogDouble tSolver: Factorize + fwd/adjoint solves.

Unnamed Group

PetscReal *MeBuf: Per-cell elemental-matrix work buffers used by computeGradientContribution (numDofInCell × numDofInCell each, zeroed per cell inside the loop). Contiguous backing (numDof²).

PetscReal *KeBuf: Contiguous backing (numDof²).

PetscReal **MeRows: Row-of-pointers view of MeBuf.

PetscReal **KeRows: Row-of-pointers view of KeBuf.

Unnamed Group

Vec bVec: Per-iteration reusable Vecs (global, sized on dm). DMCreateGlobalVector caches storage on the DM, but moving them out of the callback still saves the 4× DMCreateGlobalVector + 2× VecCreate calls per iter and keeps the same buffers warm across iterations. RHS workspace (forward solve).

Vec xVec: Forward solution.

Vec nBvec: Adjoint RHS.

Vec nxVec: Adjoint solution.

Vec ExRecvVec: Ex at receivers, parallel.

Vec wcdtDvec: Weighted residual workspace, parallel.

Unnamed Group

Vec *Bvec_per_freq: Per-frequency precomputed inputs. The RHS Vec, weights Vec, and observed-Ex row depend only on the source, the observed data, and iparams->errorLevel - all constant across L-BFGS iterations. numFreqs, sized on dm.

Vec *Wf_per_freq: numFreqs, sized seq Nrec.

Vec *dObsRow_per_freq: numFreqs, sized seq Nrec.

PetscInt numFreqsAlloc: Size of the arrays above.

Unnamed Group

Mat Kmat: Cached LHS matrices. K (stiffness, σ-independent) and G_BDDC (topological vertex incidence, σ-independent) are built ONCE at setup via assembleCsemKandM and reused across all L-BFGS iters. Ms keeps the sparsity from that build but its values are refilled every callback via assembleCsemMsRefill against the current σ. Curl-curl stiffness (σ-independent).

Mat Msmat: Mass-σ matrix (refilled per callback).

Mat Gmat_BDDC: Topological discrete gradient (BDDC hint).

Unnamed Group

Mat *Avec_per_freq: Persistent per-frequency system matrices and solvers. The sparsity pattern of A_f = K - iωμ·Ms is invariant across both frequency and L-BFGS iteration (SAME_NONZERO_PATTERN), so each Avec_per_freq[f] is allocated ONCE and refilled in place every objgrad evaluation, and each ksp_per_freq[f] is created ONCE bound to it. Reusing the KSP keeps the (expensive) symbolic factorization / BDDC topological setup across all iterations; only the cheap numeric refactorization is redone when the matrix values change. Created by setupInversionWorkspace, freed by destroyInversionWorkspace. numFreqs, sized on dm.

KSP *ksp_per_freq: numFreqs.

Public Members

const imParams *iparams: Inversion parameters.

DM dm: H(curl) DM.

DM dmConductivity: DM for conductivity (3 DOF/cell).

DM dmInversion: DM for inversion (1 DOF/cell).

Grid grid: Finite-element grid descriptor.

Vec conductivity: Current sigma model (local).

Vec X0: Initial log(rho) (global, 1 DOF/cell).

Vec DfDm: Gradient workspace (local, 1 DOF/cell).

Mat dObs: Observed Ex, Nfreq × Nrec dense.

const ReceiverInterpolationMatrices *Q: Receiver-interpolation operators.

const NeighborGraph *graph: Neighbor smoothing graph.

Vec notFixedMaskGlobal: 0 at fixed cells, 1 elsewhere (global, dmInversion).

Vec notFixedMaskLocal: Same as above but local (dmInversion).

PetscReal *allRMS: Output: RMS per iteration.

PetscReal lastRMS: RMS at most-recent callback call.

PetscReal lastDataMisfit: Data term at most-recent eval.

PetscReal lastRegTerm: Tikhonov term at most-recent eval.

PetscInt iterCount: Total objgrad calls.

PetscInt acceptedIter: Accepted L-BFGS steps.

Quadrature3D quad3d: 3D quadrature for elemental mass-matrix integration. Depends only on iparams->fm.nord (constant across the run). 3D quadrature rule.

PetscBool quad3dInited: PETSC_TRUE once quad3d is filled.