pycsamt.forward.maxwell.mt3d#

Research-only, small-grid 3-D MT finite-difference adapter.

Status: research-only, not for production use. Per docs/source/development/adr/AI-INVERSION-M6-3D-ADR.md (the M6 architecture-decision record), an in-house 3-D Maxwell solver failed the feasibility gate for production training: direct sparse solves do not scale to realistic 3-D mesh sizes (empirically confirmed there; extrapolates to tens of minutes per single solve at production resolution). This module exists for small-grid research and unit validation only. Production 3-D work should use pycsamt.forward.maxwell.modem3d (a trusted external backend) once available.

Physics#

MT3DAdapter discretizes the frequency-domain equation

\[\nabla \times \nabla \times E + i \omega \mu_0 \sigma E = 0\]

on a Cartesian Yee (staggered edge/face) grid, with E on cell edges and H on cell faces. Cell widths may be non-uniform per axis (a padded/graded tensor mesh, like pycsamt.forward.maxwell.mt2d already uses) — see “Non-uniform mesh support” below. The discrete curl operators (_curl_e2h(), _curl_h2e()) are built and verified independently of the physics: applying _curl_e2h() to a uniform field gives exactly zero, to a linear field gives the exact analytic curl, and face_divergence @ curl_e2h is the exact zero matrix (the topological identity div(curl(E)) = 0), on both uniform and non-uniform grids — see pycsamt/forward/tests/test_maxwell_mt3d.py. The overall sign convention (H = -curl(E) / (i omega mu0), assembled as (curl_h2e @ curl_e2h + i*omega*mu0*diag(sigma_edge)) @ E = 0) was anchored empirically against the analytic half-space limit, not assumed from a textbook derivation, because this codebase’s depth-increases-downward coordinate convention does not necessarily share a textbook right-handed frame’s sign. Two independent horizontal polarizations are solved per frequency (boundary Ex driven and boundary Ey driven); receivers combine both to recover the full impedance tensor.

Non-uniform mesh support#

supports_nonuniform_mesh=True. Earlier versions of this module supported uniform cell spacing only, and the documented consequence was that layered_earth_benchmark() failed by 30-45% even after mesh refinement. Investigating that failure (see git history / the AI-inversion project memory for the session this was diagnosed in) found it was not a physics defect in the boundary-condition approximation itself: the per-column decay in _column_decay() is exact for a laterally uniform half-space regardless of domain size (that is why the half-space benchmark always passed), and for a genuinely layered earth its error stays negligible wherever the field has already decayed close to zero by the domain edge — the same “boundary far enough away” argument pycsamt.forward.maxwell.mt2d documents for its own boundary treatment. The real problem was that a uniform-only mesh cannot reach several skin depths of lateral/vertical extent and resolve a few-hundred-metre layer interface within the maximum_cells budget at the same time, because 3-D cell count scales as the cube of resolution — refining resolution while keeping the same small domain does not help, which is exactly the “refining the mesh does not reduce them” symptom that was previously (incorrectly) attributed to the boundary approximation. Generalizing the curl operators (_curl_e2h(), _curl_h2e()) and edge-conductivity averaging (_sigma_on_edges()) to per-axis, non-uniform cell widths lets a padded mesh (fine cells near the receiver/structure, geometrically growing cells outward) reach the same physical extent at a fraction of the cell cost, matching what pycsamt.forward.maxwell.mt2d already does. _curl_h2e() in particular must divide by the dual-grid spacing (the average of the two neighbouring primal cell widths) rather than either cell’s own width — the two coincide on a uniform mesh, which is why the uniform-only implementation could get away without this distinction.

Deliberate scope reductions (all enforced, not just documented)#

  • Small grids only. maximum_cells (default 6,000, overridable) rejects larger problems outright rather than silently taking an impractical amount of time, operationalizing the ADR’s feasibility finding directly in code. A padded, non-uniform mesh (see above) substantially relaxes what this budget can achieve, but the ADR’s underlying conclusion — that a direct solve does not scale to realistic production 3-D mesh sizes — is unaffected: this module remains research/small-grid only.

  • Direct sparse solve only (scipy.sparse.linalg.spsolve()); no iterative solver or preconditioner, per the ADR’s conclusion that a real one is a research problem in itself.

  • Surface receivers only, within the mesh’s horizontal extent, checked in MT3DAdapter.assess().

  • Boundary conditions are a per-column exponential approximation (_column_decay(), the same level of rigor as pycsamt.forward.em2d’s own _ey_1d_profile): each boundary edge uses the exponential decay implied by its nearest cell column’s single local layer, not a full multi-layer recursion. As described above, this is exact for a half-space and negligible-error elsewhere once the domain is wide/deep enough, not a source of the previously measured layered-earth bias. Ez is fixed at zero on every boundary (plane-wave incidence has no driven vertical field far from the domain interior).

  • Isotropic conductivity, vacuum permeability, exp(+iwt) only, no inactive-cell/topography support — same restrictions as MT2DAdapter.

Measured accuracy#

On a padded, non-uniform 16x16x16 (4,096-cell) research-scale grid — fine (150 m / 120 m) core cells near the receiver and shallow interfaces, geometrically padded out to ~30 km laterally and ~8 km vertically, comfortably beyond the deepest layer’s skin depth at the lowest benchmark frequency — this adapter passes the default BenchmarkThresholds for both half_space_benchmark() (~1.9% normalized RMS, ~1.3 degree phase error) and layered_earth_benchmark() (~3.5% normalized RMS, ~4.5% amplitude error, <1 degree phase error); see pycsamt/forward/tests/test_maxwell_mt3d.py. _VERIFIED_BENCHMARKS now includes both.

Classes

MT3DAdapter(*[, version, policy, max_cells])

Research-only 3-D MT adapter (see module docstring for scope).
class pycsamt.forward.maxwell.mt3d.MT3DAdapter(*, version='1.0-research', policy=None, max_cells=6000)[source]

Bases: BaseMaxwellAdapter

Research-only 3-D MT adapter (see module docstring for scope).

Parameters:

  • version (str, default="1.0-research") – Adapter version reported in every ForwardResult.

  • policy (AdapterPolicy or None, optional) – Solver-independent result acceptance policy.

  • max_cells (int, default=6000) – Safety ceiling on total mesh cells (maximum_cells). A direct sparse solve becomes impractically slow well before typical production 3-D mesh sizes; see docs/source/development/adr/AI-INVERSION-M6-3D-ADR.md. Raise this only if you have confirmed the resulting solve time is acceptable for your use.

Examples

>>> import numpy as np
>>> from pycsamt.forward.maxwell import (
...     MaxwellMesh,
...     MaxwellProblem,
...     ReceiverSet,
... )
>>> mesh = MaxwellMesh(
...     np.linspace(0, 4000, 9),
...     np.linspace(0, 3000, 11),
...     np.linspace(0, 4000, 9),
... )
>>> problem = MaxwellProblem(
...     mesh,
...     np.full(mesh.shape, 1.0 / 100.0),
...     [1.0],
...     ReceiverSet([[2000.0, 2000.0, 0.0]], ["S00"]),
...     ("zxy", "zyx"),
... )
>>> result = MT3DAdapter().solve(problem)
>>> result.shape
(1, 1, 2)
assess(problem)[source]

Assess a problem, adding this solver’s research-only checks.

Parameters:

problem (MaxwellProblem) – Candidate simulation problem.

Returns:

The generic capability report from assess(), extended with surface-receiver, horizontal-bounds, and vacuum-permeability checks.

Return type:

CompatibilityReport

Examples

See MT3DAdapter for a complete solve example.