mom — Galerkin Method-of-Moments

Physical context

The image-charge backends (image, image_2layer, image_nlayer, cim) all share an additional simplifying assumption: the current distribution along an electrode is uniform per unit length. This assumption is convenient — it lets the cluster constraints be written directly on the per-electrode total current — but it is a genuine approximation. The true current distribution along a thin wire embedded in soil is non-uniform: it is concentrated at the electrode ends and at junctions, with a residual that scales as \(\sim 5\,\%\) of the integrated input impedance for the canonical Sunde rod.

The Galerkin Method-of-Moments (MoM) drops this assumption. It solves a full segment-level linear system for the actual per-segment currents, and only enforces the cluster equipotentiality and sum-of-currents constraints as augmenting rows. The kernel is the same layered Green's function as the matching image backend (homogeneous for \(n = 1\), Tagg / Sunde for \(n = 2\)); only the resolution scheme differs.

mom was introduced as the second engine in ADR-0001 to provide an independent resolution scheme against which the image backends could be cross-validated. The agreement was deliberately designed to be tight (within 2 %) precisely because the kernel is identical — any larger discrepancy would point to a bug in either the cluster constraint logic or the kernel evaluation.

Governing equation: boundary integral

A wire-segment system in soil satisfies the boundary integral equation

\[ \varphi(\mathbf{r}_i) \;=\; \frac{1}{4\pi} \sum_{j=1}^{N} I_j \int_{\Sigma_j} G(\mathbf{r}_i, \mathbf{r}'_j) \,dS_j, \]

with \(G\) the layered Green's function (the same one the image backends use), \(\Sigma_j\) the surface of segment \(j\), and \(\varphi(\mathbf{r}_i)\) the (unknown) cluster potential of the cluster containing segment \(i\). In matrix form, \(Z\,I = \varphi\) with the reaction matrix

\[ Z_{ij} \;=\; \frac{\rho_1}{4\pi} \cdot \frac{1}{L_i} \int_{\text{seg}_i} \! \int_{\Sigma_j} G(\mathbf{r}, \mathbf{r}')\, dS_j\, d\ell. \]

The Galerkin formulation arises by averaging \(\varphi\) over segment \(i\) — i.e. taking the same basis function as the test function. This is one of two natural choices; the alternative (point matching at the segment midpoint) is implemented in bem.

Numerical strategy

Reaction matrix assembly

The reaction matrix is built directly from the same self-kernel factory that the image backend uses:

• For HomogeneousSoil it calls _self_corrected_kernel(seg_points, seg_lengths, wire_radii, identity, rho) — feeding the identity matrix as the "currents" argument extracts the full \(N \times N\) kernel matrix in one vectorised pass.
• For TwoLayerSoil it calls _two_layer_self_kernel_factory(soil, max_terms, tol) — same trick, but with the Tagg / Sunde geometric series.

The diagonal of \(Z\) carries the analytical line self-potential \(2 \ln(L_i / a_i) / L_i \cdot \rho_1 / (4\pi)\) plus the layered image contributions. The off-diagonals use the point-source approximation (segment midpoints) for the \(1/r + 1/r_{\text{air}}\) direct + air-mirror part and for every \(K_1^n\) image term.

Galerkin solve

The cluster constraints are enforced through an augmented system

\[ \begin{bmatrix} Z & -C \\ C^{\top} & 0 \end{bmatrix} \begin{bmatrix} I_{\text{seg}} \\ \varphi_c \end{bmatrix} = \begin{bmatrix} 0 \\ I_{c,\text{in}} \end{bmatrix}, \]

with \(C\) the segment-to-cluster membership matrix (\(C_{ic} = 1\) if segment \(i\) belongs to cluster \(c\), else 0), and \(I_{c,\text{in}}\) the input current of cluster \(c\). The first \(N\) rows enforce \(\varphi_i = \varphi_c \;\forall i \in c\) on the reaction-matrix level; the next \(K\) rows enforce \(\sum_{i \in c} I_i = I_{c,\text{in}}\) for every active cluster.

The system is symmetric (since \(Z\) is symmetric and \(C\) enters once with each sign); numpy.linalg.solve handles it directly.

Inactive segments (those whose cluster has zero input current) are excluded from the unknowns to keep the matrix small.

Real / imaginary split

Because \(Z\) is real and the only complex contribution to the right-hand side is the source phase, the system is solved separately for the real and imaginary part of every cluster current. This avoids complex arithmetic in the LU solve.

Validity envelope

Property	Range / value
Soil model	HomogeneousSoil or TwoLayerSoil
Frequency	quasi-static, \(f < 1\,\text{kHz}\)
Electrode placement	every segment in the upper layer (for TwoLayerSoil)
Wire radius / segment ratio	thin-wire, \(a \ll L_i\)
Mesh size \(N\)	\(\le 1000\) at acceptable runtime

For \(n \ge 3\) soils, mom does not currently apply — the multilayer kernel is delegated to cim, mom_sommerfeld, or bem. ADR-0002 keeps the option of extending mom with the CIM kernel as a future action item, but the existing mom_sommerfeld engine already covers that role.

Convergence and cost

• Resolution accuracy. The Galerkin scheme typically gets within 1 % of the Dwight reference at the canonical 1.5 m / 5 mm rod, eliminating the 4–5 % uniform-current bias of the image backend. The improvement carries over to layered soils: mom is consistently the closest match to the literature reference values among the segment-discretised engines.
• Computational cost. \(O(N^2)\) matrix build, \(O((N + K)^3)\) LU solve. For typical cases with \(N \le 1000\) the runtime is dominated by the matrix build (a few hundred milliseconds in the 2-layer case with \(\sim 200\) image terms).
• Reduction. At \(K_1 = 0\) the engine collapses bit-exactly to the homogeneous Galerkin solution.

Cross-validation notes

Counterpart	Expected agreement	What is checked
image (homogeneous)	\(\le 2\,\%\)	uniform-current vs. Galerkin
image_2layer (\(n = 2\))	\(\le 2\,\%\)	same kernel, different test function
mom_sommerfeld (\(n = 2\))	\(\le 5\,\%\)	quadrature reproduces the same physics
bem (any \(n\))	\(\le 5\,\%\)	Galerkin vs. collocation on the same kernel
Sunde / Dwight closed forms	\(\le 5\,\%\)	tighter bound than the image backend

mom is the independent resolution-scheme reference in the cross-validation matrix. When image and mom agree on a given problem, the kernel is implemented correctly; when they disagree, the disagreement is reproducible and points to a uniform-current bias issue.

References

• Harrington, R. F. (1968). Field Computation by Moment Methods, Macmillan. The MoM textbook.
• Meliopoulos, A. P. S., Xia, F., Joy, E. B., Cokkinides, G. J. (1993). An advanced computer model for grounding system analysis. IEEE PWRD 8(1). MoM for grounding systems with cluster constraints.
• Sunde, E. D. (1968). Earth Conduction Effects in Transmission Systems, Dover, sect. 3.5. The line self-potential (\(2 \ln(L/a) / L\)) used on the diagonal.
• Dawalibi, F. P. & Barbeito, N. (1991). Measurements and computations of the performance of grounding systems buried in multilayer soils. IEEE PWRD 6(4). Reference values for cross-validation.

Example

import groundfield as gf

soil = gf.HomogeneousSoil(resistivity=100.0)
world = gf.create_world(soil=soil)
gf.create_electrode(world, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=1.5)
gf.create_source(world, attached_to="g1", magnitude=1.0)

engine = gf.create_engine(backend="mom",
                          segment_length=0.05,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])

API reference

mom

Method-of-Moments backend (mom).

Independent second engine used to cross-validate the image backends (see ADR-0001). The image backends approximate the current distribution along an electrode by a uniform per unit length profile; this backend solves for the actual non-uniform per-segment current distribution through a Galerkin scheme.

Method

1. Discretise every electrode into wire segments (same routine as the image backends).
2. Assemble the full \(N \times N\) segment-level reaction matrix \(Z\) using the same Green's-function kernel as the corresponding image backend (so for a HomogeneousSoil world mom shares its physics with image; for a TwoLayerSoil world it shares it with image_2layer). The diagonal of \(Z\) carries the analytical line-self-potential.
3. Solve the linear system $$ \begin{bmatrix} Z & -C \ C^{\top} & 0 \end{bmatrix} \begin{bmatrix} I_\text{seg} \ \varphi_c \end{bmatrix} = \begin{bmatrix} 0 \ I_\text{in} \end{bmatrix}, $$ where \(C\) is the segment-to-cluster membership matrix.

The constraint enforces equal segment potential within each galvanic cluster; the additional row sums the segment currents to the input current per cluster. The system is symmetric and dense; :func:numpy.linalg.solve handles it directly.

1. Compute the average cluster potential and produce a :class:FieldResult exactly like the image backends, so that :func:groundfield.compare_engines can compare the two side by side.

What this means in practice

• For homogeneous soil, mom agrees with image to within a few per cent and lies closer to the Sunde reference (the average- potential method of image carries a residual ~5 % bias from the uniform-current assumption, which the Galerkin solve removes).
• For 2-layer soil, mom and image_2layer agree within the same envelope. In the limit \(\rho_2 = \rho_1\) both engines collapse onto the homogeneous solution.

Limitations

• Performance is \(O(N^3)\) for the LU solve plus \(O(N^2)\) for the matrix build. For typical geometries N stays well below 1 000, so the runtime is acceptable. Larger meshes will need a preconditioned iterative solver later.
• This backend deliberately re-uses the existing Green's-function kernels rather than introducing a Sommerfeld quadrature. A truly independent layered Green's function (option B in ADR-0001) would give a stronger cross-check; this lighter version is the pragmatic first step.

solve_mom

solve_mom(
    world: "World",
    engine: "Engine",
    *,
    two_layer_max_terms: int | None = None,
    two_layer_tol: float | None = None
) -> FieldResult

Galerkin Method-of-Moments backend.

Supports :class:HomogeneousSoil and :class:TwoLayerSoil. The Green's-function kernel is taken from the matching image backend; only the resolution scheme (per-segment Galerkin instead of per-electrode average potential) differs.

Parameters:

Name	Type	Description	Default
world	'World'	World to evaluate.	required
engine	'Engine'	Engine configuration; engine.segment_length controls the discretisation.	required
two_layer_max_terms	int | None	Truncation parameters for the Tagg/Sunde series, only used when world.soil is a :class:TwoLayerSoil. None (default) takes engine.image_max_terms / engine.image_series_tol, so the engine knob reaches this backend consistently with image_2layer.	None
two_layer_tol	int | None	Truncation parameters for the Tagg/Sunde series, only used when world.soil is a :class:TwoLayerSoil. None (default) takes engine.image_max_terms / engine.image_series_tol, so the engine knob reaches this backend consistently with image_2layer.	None

Returns:

Type	Description
FieldResult	Result object compatible with :class:FieldResult from the image backends. metadata['solver'] is set to 'galerkin'.

Source code in src/groundfield/solver/mom.py

def solve_mom(
    world: "World",
    engine: "Engine",
    *,
    two_layer_max_terms: int | None = None,
    two_layer_tol: float | None = None,
) -> FieldResult:
    """Galerkin Method-of-Moments backend.

    Supports :class:`HomogeneousSoil` and :class:`TwoLayerSoil`. The
    Green's-function kernel is taken from the matching image backend;
    only the resolution scheme (per-segment Galerkin instead of
    per-electrode average potential) differs.

    Parameters
    ----------
    world
        World to evaluate.
    engine
        Engine configuration; ``engine.segment_length`` controls the
        discretisation.
    two_layer_max_terms, two_layer_tol
        Truncation parameters for the Tagg/Sunde series, only used
        when ``world.soil`` is a :class:`TwoLayerSoil`. ``None``
        (default) takes ``engine.image_max_terms`` /
        ``engine.image_series_tol``, so the engine knob reaches this
        backend consistently with ``image_2layer``.

    Returns
    -------
    FieldResult
        Result object compatible with :class:`FieldResult` from the
        image backends. ``metadata['solver']`` is set to ``'galerkin'``.
    """
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")
    if not isinstance(world.soil, (HomogeneousSoil, TwoLayerSoil)):
        raise TypeError(
            "Backend 'mom' currently supports HomogeneousSoil and "
            f"TwoLayerSoil only. Got: {type(world.soil).__name__}."
        )
    _reject_concrete_shells(world, "mom")
    _warn_ignored_sources(world, "mom")
    if two_layer_max_terms is None:
        two_layer_max_terms = engine.image_max_terms
    if two_layer_tol is None:
        two_layer_tol = engine.image_series_tol

    ds = engine.segment_length
    _log.info(
        "mom: soil=%s, segment_length=%.3f",
        type(world.soil).__name__, ds,
    )

    # 1) Discretisation — identical to the image backends; for
    #    two-layer soil segments are split at the interface
    #    (ADR-0007 step 1, audit 2026-07-08 WP-D2).
    interfaces = (
        (world.soil.h_1,) if isinstance(world.soil, TwoLayerSoil) else None
    )
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=interfaces)
        elec_to_segidx[e.name] = list(range(len(all_segments),
                                            len(all_segments) + len(segs)))
        all_segments.extend(segs)

    # 2) Per-electrode input currents from the configured sources.
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    for src in world.sources:
        if src.kind != "current":
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex

    # 3) Cluster building (ideal conductors only), finite-impedance
    #    branch list, plus the distributed-conductor topology
    #    (ADR-0003) — all fed into the augmented Galerkin system.
    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])

    # 4) Cross-layer support (audit 2026-07-08, WP-D3): the shared
    #    self-kernel factory dispatches to the rigorous ADR-0007
    #    Sommerfeld path when any segment sits at or below h_1, so
    #    the historic hard rejection is gone. mom is the engine that
    #    *solves* the per-segment current distribution — exactly the
    #    tool to quantify the image family's uniform-leakage bias
    #    across the interface.
    if isinstance(world.soil, TwoLayerSoil):
        z_max = seg_points[:, 2].max()
        if z_max >= world.soil.h_1:
            _log.info(
                "mom: cross-layer geometry (z_max = %.3f m >= h_1 = "
                "%.3f m) — using the ADR-0007 Sommerfeld kernel.",
                z_max, world.soil.h_1,
            )

    # 5) Assemble the reaction matrix Z (kernel depends on soil model).
    n_terms_used = 0
    if isinstance(world.soil, HomogeneousSoil):
        Z = _build_Z_homogeneous(
            seg_points, seg_lengths, wire_radii, world.soil.resistivity
        )
    else:
        # TwoLayerSoil
        Z = _build_Z_two_layer(
            seg_points, seg_lengths, wire_radii,
            world.soil, two_layer_max_terms, two_layer_tol,
        )
        # Approximate the term count from the same offset generator.
        from groundfield.solver.image_2layer import _two_layer_image_offsets

        _, _, _, n_terms_used = _two_layer_image_offsets(
            world.soil.reflection_coefficient, world.soil.h_1,
            two_layer_max_terms, two_layer_tol,
        )

    # 6) Frequency loop. With no inductive coupling, solve once and
    #    replicate across frequencies; otherwise loop.
    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    if isinstance(world.soil, HomogeneousSoil):
        rho_self = world.soil.resistivity
        phi_kernel = lambda sc_real: _self_corrected_kernel(
            seg_points, seg_lengths, wire_radii, sc_real, rho_self,
        )
    else:
        self_kernel = _two_layer_self_kernel_factory(
            world.soil, two_layer_max_terms, two_layer_tol,
            allow_cross_layer=True,
        )
        phi_kernel = lambda sc_real: self_kernel(
            seg_points, seg_lengths, wire_radii, sc_real,
        )

    def _solve_at(omega: float) -> tuple[np.ndarray, np.ndarray]:
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        sc, _ = _galerkin_solve(
            Z=Z,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            elec_to_segidx=elec_to_segidx,
            n_segments=n_segments,
            finite_branches=finite_branches,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
        )
        ph = np.zeros(n_segments, dtype=complex)
        if sc.any():
            ph = phi_kernel(sc.real) + 1j * phi_kernel(sc.imag)
        return sc, ph

    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            sc, ph = _solve_at(omega)
            sc_per_freq.append(sc)
            phi_per_freq.append(ph)
    else:
        sc, ph = _solve_at(0.0)
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(n_freq)
        ]
        i_list = [
            complex(sc_per_freq[k][idxs].sum()) for k in range(n_freq)
        ]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    # 8) Point-source list and cluster map for post-processing.
    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]
    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    # 9) Soil-specific metadata.
    if isinstance(world.soil, HomogeneousSoil):
        soil_resistivity = float(world.soil.resistivity)
        soil_meta: dict[str, float | int | bool] = {}
    else:
        K = world.soil.reflection_coefficient
        soil_resistivity = float(world.soil.rho_1)
        converged = (
            _series_tail_bound(abs(K), n_terms_used) < two_layer_tol
            if n_terms_used else True
        )
        if not converged:
            warnings.warn(
                f"mom: the Tagg/Sunde image series was truncated at "
                f"max_terms={two_layer_max_terms} with tail bound > "
                f"tol={two_layer_tol:.2e} (|K| = {abs(K):.3f}). Raise "
                "Engine.image_max_terms for high layer contrasts.",
                SeriesTruncationWarning,
                stacklevel=2,
            )
        soil_meta = {
            "K": float(K),
            "rho_1": float(world.soil.rho_1),
            "rho_2": float(world.soil.rho_2),
            "h_1": float(world.soil.h_1),
            "n_terms_used": n_terms_used,
            "converged": bool(converged),
        }

    return FieldResult(
        backend="mom",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=soil_resistivity,
        soil=world.soil,
        clusters=cluster_members,
        metadata={
            "world_name": world.name,
            "n_segments": n_segments,
            "segment_length": ds,
            "solver": "galerkin",
            "stub": False,
            "earth_inductive_model": earth_inductive_model,
            **soil_meta,
            **_carson_skin_metadata(
                has_inductance, sigma_earth_for_carson, world.soil,
                engine.frequencies,
            ),
            **(
                {
                    "conductor_node_currents": conductor_currents,
                    "conductor_node_potentials": conductor_potentials,
                }
                if conductor_currents else {}
            ),
        },
    )

Related material

• ADR-0001 — original methodology decision.