`mom_sommerfeld` — Galerkin MoM with direct Sommerfeld quadrature¶

Physical context¶

Every other layered backend in the family relies on a closed-form representation of the layered Green's function:

image_2layer — geometric Tagg / Sunde series.
image_nlayer — dispatcher to the closed-form variants.
cim — matrix-pencil approximation by complex exponentials.
bem — collocation on top of the CIM kernel.

mom_sommerfeld is methodologically distinct: it evaluates the Sommerfeld integral numerically, point by point, with adaptive Gauss–Kronrod quadrature. The recursive $\Gamma_1(\lambda)$ enters the integrand as is — no expansion, no fit. The price is speed (a single $N \times N$ reaction matrix can take seconds rather than milliseconds), but the result is an absolute reference: any disagreement between the closed-form engines and mom_sommerfeld indicates a fit / series accuracy issue, not a kernel bug.

In ADR-0002 this is the engine that anchors the cross-validation envelope for the $n \ge 3$ regime, and provides the independent methodology that the closed-form layered family lacks.

Governing equation: layered Sommerfeld integral¶

For a top-layer source at depth $z_s$ and a top-layer field point at $(s, z)$, the full layered Green's function is the Sommerfeld integral with multiple-reflection multiplier:

\[ G(s, z, z_s) \;=\; \int_0^{\infty} \frac{ e^{-\lambda |z - z_s|} + e^{-\lambda (z + z_s)} + \Gamma_1(\lambda)\, e^{-\lambda (2 h_1 - |z - z_s|)} + \Gamma_1(\lambda)\, e^{-\lambda (2 h_1 - z - z_s)} } { 1 - \Gamma_1(\lambda)\, e^{-2\lambda h_1} }\, J_0(\lambda s)\, d\lambda. \]

The integrand carries:

Two direct exponentials — the source and its air mirror.
Two reflected exponentials — the source's reflection at $z = h_1$ and the air mirror's reflection at $z = h_1$.
A multiplexion factor $1 / (1 - \Gamma_1(\lambda) e^{-2\lambda h_1})$ that captures the infinite back-and-forth bouncing between the air boundary ($R_{\text{air}} = +1$) and the layer interface ($\Gamma_1$).

For $\Gamma_1 \to 0$ (homogeneous) the integral reduces to $1/r + 1/r_{\text{img}}$. For $\Gamma_1 \equiv K_1$ (2-layer) the geometric expansion of the multiplier reproduces the Tagg / Sunde series. For $n \ge 3$ no closed form is available — the engine just integrates.

Numerical strategy¶

Quadrature truncation¶

The integrand decays exponentially in $\lambda$ at the rate of the fastest decaying exponential $\min(|z - z_s|, z + z_s, 2 h_1 - z - z_s, 2 h_1 - |z - z_s|)$. We bound the upper limit at $\lambda_{\max} = \texttt{lambda\_max\_factor} / \bar h$ with $\bar h = \min(h_1, s + z + z_s + \epsilon)$ — i.e. the characteristic length scale of the geometry, and a default factor of 200. This puts the residual integrand contribution above $\lambda_{\max}$ at $e^{-200} \sim 10^{-87}$, well below any practical tolerance.

Inside that bound, scipy.integrate.quad (Gauss–Kronrod adaptive quadrature with a 21-point rule) handles the integration to absolute / relative tolerances $10^{-9} / 10^{-7}$. The Bessel function $J_0(\lambda s)$ oscillates with period $2\pi/s$; the adaptive rule subdivides until the per-subinterval rule and a higher-order rule agree to within tolerance, so the oscillation is captured automatically.

Why not the Zou et al. complex contour?¶

Zou, Du & Zhou (2015) proposed a deformed integration contour in the complex $\lambda$-plane that breaks the Bessel oscillation and makes the integrand non-oscillatory along the contour. The current implementation uses the simpler real-axis adaptive quadrature. Pragmatic justification:

For $s \lesssim 100\,\text{m}$ (the typical range) the real-axis adaptive rule converges in 200–400 kernel evaluations per pair, which is acceptable for the engine's role as cross-check reference.
The Zou contour requires careful handling of the branch points in $\Gamma_1$ for $n \ge 3$, which would make the implementation significantly larger.
The cross-validation tests show that the real-axis rule already agrees with the closed-form engines to better than 1 % across the typical contrast range.

The module's docstring still cites Zou et al. as the reference for the idea; switching to a complex contour is a documented future optimisation (ADR-0002 action items).

Reaction matrix assembly¶

mom_sommerfeld builds the $N \times N$ reaction matrix in three modes:

Homogeneous ($n = 1$). Falls back on the closed-form homogeneous self-kernel — quadrature is unnecessary and would introduce numerical noise.
Two-layer ($n = 2$). Uses the Tagg / Sunde self-kernel for the diagonal (the closed-form multi-image self-action) and the direct Sommerfeld quadrature for the off-diagonals. The two are consistent — each off-diagonal entry is the same physics evaluated by quadrature, and the diagonal is the closed-form geometric series for the same physics.
Multilayer ($n \ge 3$). Uses the homogeneous line self-potential plus a layered self-correction obtained from the reflection-only Sommerfeld integral evaluated at $s = 0$, $|z - z_s| = 0$ (the integrand has no $1/r$ singularity in this reflection-only form, so it can be integrated point-on-source).

Galerkin solve¶

The actual linear-system solve reuses the _galerkin_solve helper from mom — the cluster augmenting rows are identical; only the kernel changes.

Validity envelope¶

Property	Range / value
Soil model	`HomogeneousSoil`, `TwoLayerSoil`, `MultiLayerSoil`
Frequency	quasi-static, $f < 1\,\text{kHz}$
Electrode placement	every segment in the upper layer
Quadrature tolerances	$\text{epsabs} = 10^{-9}$, $\text{epsrel} = 10^{-7}$
Quadrature limit	400 subdivisions
$\lambda_{\max}$ factor	200 / characteristic length (default)
Mesh size $N$	practical limit $\sim 200$ segments at acceptable runtime

The runtime cost of the engine grows as $O(N^2)$ Sommerfeld quadrature evaluations, each costing a few hundred kernel evaluations. For $N = 200$ on a 2-layer world the engine completes in a few seconds; for $N = 1000$ it would take of the order of a minute. This is acceptable for cross-validation but not for production sweeps — use image_2layer / cim / bem there.

Convergence and cost¶

Quadrature accuracy. With the default tolerances the per-pair error is $\sim 10^{-8}$ relative to the analytical closed forms. The cluster-impedance error in the assembled result stays below 1 % of the closed-form reference.
Per-pair cost. Adaptive Gauss–Kronrod typically converges in 100–400 kernel evaluations per (field, source) pair.
Reduction. At $\Gamma_1 \to 0$ the engine short-circuits to $1/r + 1/r_{\text{img}}$ — bit-exact match with image.

Cross-validation notes¶

Counterpart	Expected agreement	What is checked
`image` ($n = 1$)	bit-exact	quadrature short-circuits to closed form
`image_2layer` ($n = 2$)	$\le 5\,\%$	quadrature reproduces the geometric series
`cim` (any $n$)	$\le 5\,\%$	matrix-pencil fit accuracy
`bem` (any $n$)	$\le 5\,\%$	collocation on the same physics
`mom` ($n \le 2$)	$\le 2\,\%$	same Galerkin scheme, different kernel evaluation

The engine's role is the absolute reference. Whenever the closed-form layered engines disagree, the disagreement is measured against mom_sommerfeld. The cross-engine validation harness report agreement tables relative to this engine.

References¶

Sommerfeld, A. (1909). Über die Ausbreitung der Wellen in der drahtlosen Telegraphie. Annalen der Physik 28. The original paper introducing the integral.
Zou, J., Du, X. & Zhou, C. (2015). Fast calculation of the Green function of a point current source in a horizontal layered soil with a new complex path. IEEE Trans. Magn. 51(3). The reference for a deformed contour that the current implementation does not use, but which is on the optimisation roadmap.
Piessens, R. et al. (1983). QUADPACK, Springer. The Gauss–Kronrod adaptive quadrature implementation that scipy wraps.
Dwight, H. B. (1936). Calculation of resistances to ground. Reference DC resistances used in the cross-validation tests.

Example¶

import groundfield as gf

soil = gf.TwoLayerSoil(rho_1=100.0, rho_2=500.0, h_1=2.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)

# Reference engine — slow but methodologically independent.
engine = gf.create_engine(backend="mom_sommerfeld",
                          segment_length=0.1,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])

API reference¶

mom_sommerfeld ¶

Method-of-Moments backend with direct Sommerfeld quadrature (mom_sommerfeld).

Mathematical / physical model

The other layered backends (image_nlayer, cim, bem) rely on a closed-form representation of the layered Green's function (real image series or complex images). This backend instead evaluates the Sommerfeld integral $$ \varphi(s, z; z_s) \;=\; \frac{\rho_1\, I}{4\pi} !\int_0^{\infty}! \bigl[ e^{-\lambda |z - z_s|} + \Gamma_1(\lambda)\, e^{-\lambda (z + z_s)}\bigr] J_0(\lambda s)\, d\lambda $$ numerically, point by point, with adaptive Gauss–Kronrod quadrature (scipy.integrate.quad). The recursive $\Gamma_1(\lambda)$ from :func:groundfield.solver._layered.reflection_gamma is used as is — no expansion, no fit. The result is therefore independent of the expansion choices in image_nlayer / cim and serves as the absolute reference inside the cross-engine validation.

Following Zou et al. 2015, the integration path is deformed into the complex plane to break the Bessel oscillation: the real-axis quadrature is replaced by a contour that combines (a) a flat real segment $[0, \lambda_0]$, (b) a deformed complex piece that decays exponentially. For the radial distances of interest in grounding (s ≲ 100 m) this contour collapses the integrand to a non-oscillatory shape and makes scipy.integrate.quad converge in a few hundred kernel evaluations.

The MoM resolution itself (Galerkin scheme with one constraint per cluster) is the same as :mod:groundfield.solver.mom; only the underlying Z-matrix kernel differs.

Validity

Quasi-static, $f < 1\,\mathrm{kHz}$.
Slow but methodologically independent. Use it as the reference in cross-engine tests on layered worlds with hard contrasts; the closed-form backends are usually within a fraction of a per cent of the quadrature result.

References

Sommerfeld, A. (1909). Über die Ausbreitung der Wellen in der drahtlosen Telegraphie. Annalen der Physik 28.
Zou, J., Du, X., & Zhou, C. (2015). Fast calculation of the Green function of a point current source in a horizontal layered soil with a new complex path. IEEE Trans. Magn. 51(3).

solve_mom_sommerfeld ¶

solve_mom_sommerfeld(
    world: "World",
    engine: "Engine",
    *,
    lambda_max_factor: float = 200.0,
    epsabs: float = 1e-09,
    epsrel: float = 1e-07
) -> FieldResult

Galerkin MoM with direct Sommerfeld quadrature for layered soil.

This is a methodologically independent backend used as the reference inside :func:groundfield.compare_engines for layered worlds with hard contrasts.

Parameters:

Name	Type	Description	Default
`world`	`'World'`	World to evaluate.	required
`engine`	`'Engine'`	Engine configuration.	required
`lambda_max_factor`	`float`	Upper bound of the quadrature in units of $1 / \bar h$.	`200.0`
`epsabs`	`float`	Quadrature tolerances.	`1e-09`
`epsrel`	`float`	Quadrature tolerances.	`1e-09`

Returns:

Type	Description
`FieldResult`

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

def solve_mom_sommerfeld(
    world: "World",
    engine: "Engine",
    *,
    lambda_max_factor: float = 200.0,
    epsabs: float = 1e-9,
    epsrel: float = 1e-7,
) -> FieldResult:
    """Galerkin MoM with direct Sommerfeld quadrature for layered soil.

    This is a methodologically independent backend used as the
    reference inside :func:`groundfield.compare_engines` for layered
    worlds with hard contrasts.

    Parameters
    ----------
    world
        World to evaluate.
    engine
        Engine configuration.
    lambda_max_factor
        Upper bound of the quadrature in units of $1 / \\bar h$.
    epsabs, epsrel
        Quadrature tolerances.

    Returns
    -------
    FieldResult
    """
    if not isinstance(world.soil, (HomogeneousSoil, TwoLayerSoil, MultiLayerSoil)):
        raise TypeError(
            "Backend 'mom_sommerfeld' supports HomogeneousSoil, "
            "TwoLayerSoil, and MultiLayerSoil. "
            f"Got: {type(world.soil).__name__}."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")
    _reject_concrete_shells(world, "mom_sommerfeld")
    _warn_ignored_sources(world, "mom_sommerfeld")

    stack = as_layer_stack(world.soil)
    ds = engine.segment_length
    _log.info(
        "mom_sommerfeld: n_layers=%d, segment_length=%.3f, lam_max_fac=%.1f",
        stack.n_layers, ds, lambda_max_factor,
    )

    # 1) Discretisation.
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    interfaces = (
        (float(stack.h[0]),) if stack.n_layers >= 2 else None
    )
    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.
    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

    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    # 2b) Distributed-conductor topology (ADR-0003) + ADR-0004
    #     inductive coupling.
    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])

    if stack.n_layers >= 3:
        z_max = seg_points[:, 2].max()
        h_1 = float(stack.h[0])
        if z_max >= h_1:
            # ADR-0007 Phase B (n≥3): not implemented — the
            # reflection-only diagonal correction below assumes all
            # segments in the top layer. Hard error instead of a
            # warning followed by wrong physics.
            raise ValueError(
                f"mom_sommerfeld: cross-layer geometry on "
                f"n_layers={stack.n_layers} is not supported "
                f"(z_max = {z_max:.3f} m >= h_1 = {h_1:.3f} m). "
                "Use backend='image_2layer' for n=2 cross-layer "
                "worlds; for n>=3 thicken the upper layer."
            )

    # 3) Z-matrix via direct Sommerfeld quadrature (series knobs for
    #    the n=2 diagonal from the engine).
    Z = _build_Z_sommerfeld(
        seg_points, seg_lengths, wire_radii, stack,
        lambda_max_factor=lambda_max_factor,
        epsabs=epsabs, epsrel=epsrel,
        max_terms=engine.image_max_terms, tol=engine.image_series_tol,
    )

    # 4) Frequency loop (Galerkin solve + Z · I_seg for phi).
    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}

    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 = Z @ sc.real + 1j * (Z @ 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

    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
        )

    metadata = {
        "world_name": world.name,
        "n_segments": n_segments,
        "segment_length": ds,
        "n_layers": int(stack.n_layers),
        "rhos": stack.rhos.tolist(),
        "h": stack.h.tolist(),
        "lambda_max_factor": float(lambda_max_factor),
        "solver": "galerkin",
        "stub": False,
        "earth_inductive_model": earth_inductive_model,
    }
    if has_inductance:
        from groundfield.coupling.carson import skin_depth

        sigma_ref = (
            sigma_earth_for_carson
            if sigma_earth_for_carson is not None
            else 1.0 / float(stack.rhos[0])
        )
        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_ref)
            for f in engine.frequencies
        }
    if conductor_currents:
        metadata["conductor_node_currents"] = conductor_currents
        metadata["conductor_node_potentials"] = conductor_potentials

    return FieldResult(
        backend="mom_sommerfeld",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=float(stack.rhos[0]),
        soil=world.soil,
        clusters=cluster_members,
        metadata=metadata,
    )

sommerfeld_kernel_value ¶

sommerfeld_kernel_value(
    stack: LayerStack,
    s: float,
    z: float,
    z_s: float,
    *,
    lambda_max_factor: float = 200.0,
    epsabs: float = 1e-09,
    epsrel: float = 1e-07
) -> float

Evaluate the layered-soil Sommerfeld kernel at one point.

Returns $$ G(s, z, z_s) \;=\; \int_0^{\infty}!\bigl[ e^{-\lambda |z - z_s|} + \Gamma_1(\lambda)\, e^{-\lambda (z + z_s)} \bigr] J_0(\lambda s)\, d\lambda. $$ The full potential is then ρ_1 · I · G / (4 π).

For a homogeneous soil (stack.n_layers == 1) $\Gamma_1 \equiv 0$ and the integral collapses to $1/r + 1/r_{\text{img}}$, with $r = \sqrt{s^2 + (z - z_s)^2}$ and $r_{\text{img}} = \sqrt{s^2 + (z + z_s)^2}$. We short-circuit that case to keep the homogeneous limit bit-exact.