Add sparse near-field inductance matrix

CHANGELOG.md

@@ -1,5 +1,15 @@
 # Changelog
 
+## 13.0.0 2026-08-19
+
+* Rust
+    * Add near-field interaction map to hierarchical diagnostics outputs
+    * !Add option for hierarchical evaluator to skip evaluating kernels for near-field, far-field, or both
+    * !Update Rust one-shot hierarchical solver convenience functions for feature parity with Python ones
+    * Add sparse inductance matrix for linear filaments using CSC interaction map
+* Python
+    * Plumb in bindings to new `skip` option and sparse inductance matrix
+
 ## 12.1.0 2026-08-17
 
 * Rust

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "cfsem"
-version = "12.1.0"
+version = "13.0.0"
 edition = "2024"
 authors = ["Commonwealth Fusion Systems <jlogan@cfs.energy>"]
 license = "MIT"

benches/linear_filament.rs

@@ -6,6 +6,7 @@ use cfsem::physics::hierarchical::{
 };
 use cfsem::physics::linear_filament::{
     flux_density_linear_filament, flux_density_linear_filament_par,
+    inductance_linear_filaments_matrix_par, inductance_linear_filaments_sparse_csc_par,
     vector_potential_linear_filament, vector_potential_linear_filament_par,
 };
 use criterion::*;
@@ -186,6 +187,7 @@ fn bench_flux_density_linear_filament(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 false,
+                                false,
                                 (&mut bx, &mut by, &mut bz),
                             )
                             .unwrap(),
@@ -216,6 +218,7 @@ fn bench_flux_density_linear_filament(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 true,
+                                false,
                                 (&mut bx, &mut by, &mut bz),
                             )
                             .unwrap(),
@@ -322,6 +325,7 @@ fn bench_vector_potential_linear_filament(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 false,
+                                false,
                                 (&mut ax, &mut ay, &mut az),
                             )
                             .unwrap(),
@@ -353,6 +357,7 @@ fn bench_vector_potential_linear_filament(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 true,
+                                false,
                                 (&mut ax, &mut ay, &mut az),
                             )
                             .unwrap(),
@@ -366,6 +371,75 @@ fn bench_vector_potential_linear_filament(c: &mut Criterion) {
     group.finish();
 }
 
+fn bench_sparse_inductance(c: &mut Criterion) {
+    const NSEGMENT: usize = 512;
+    const HALF_BANDWIDTH: usize = 8;
+
+    let input = circular_loop_linear_filament_bench_input(NSEGMENT, NSEGMENT);
+    let mut row_indices = Vec::new();
+    let mut column_pointers = Vec::with_capacity(NSEGMENT + 1);
+    column_pointers.push(0);
+    for target in 0..NSEGMENT {
+        for source in 0..NSEGMENT {
+            let separation = source.abs_diff(target);
+            let periodic_separation = separation.min(NSEGMENT - separation);
+            if periodic_separation <= HALF_BANDWIDTH {
+                row_indices.push(source);
+            }
+        }
+        column_pointers.push(row_indices.len());
+    }
+
+    let xyz = (&input.xfil[..], &input.yfil[..], &input.zfil[..]);
+    let dlxyz = (&input.dlxfil[..], &input.dlyfil[..], &input.dlzfil[..]);
+    let sparse_fraction = row_indices.len() as f64 / (NSEGMENT * NSEGMENT) as f64;
+    let mut group = c.benchmark_group("Linear Filament Inductance Matrix");
+    group.sample_size(10);
+    group.measurement_time(Duration::from_secs(5));
+
+    group.throughput(Throughput::Elements((NSEGMENT * NSEGMENT) as u64));
+    group.bench_function("Dense GL3, Parallel", |b| {
+        let mut out = vec![0.0; NSEGMENT * NSEGMENT];
+        b.iter(|| {
+            black_box(
+                inductance_linear_filaments_matrix_par(
+                    xyz,
+                    dlxyz,
+                    xyz,
+                    dlxyz,
+                    &input.wire_radius,
+                    &mut out,
+                )
+                .unwrap(),
+            )
+        });
+    });
+
+    group.throughput(Throughput::Elements(row_indices.len() as u64));
+    group.bench_function(
+        format!("Sparse GL3, Parallel ({:.1}% nnz)", 100.0 * sparse_fraction),
+        |b| {
+            let mut out = vec![0.0; row_indices.len()];
+            b.iter(|| {
+                black_box(
+                    inductance_linear_filaments_sparse_csc_par(
+                        xyz,
+                        dlxyz,
+                        xyz,
+                        dlxyz,
+                        &input.wire_radius,
+                        &row_indices,
+                        &column_pointers,
+                        &mut out,
+                    )
+                    .unwrap(),
+                )
+            });
+        },
+    );
+    group.finish();
+}
+
 criterion_group!(
     group_bench_flux_density_linear_filament,
     bench_flux_density_linear_filament
@@ -374,8 +448,10 @@ criterion_group!(
     group_bench_vector_potential_linear_filament,
     bench_vector_potential_linear_filament
 );
+criterion_group!(group_bench_sparse_inductance, bench_sparse_inductance);
 
 criterion_main!(
     group_bench_flux_density_linear_filament,
-    group_bench_vector_potential_linear_filament
+    group_bench_vector_potential_linear_filament,
+    group_bench_sparse_inductance
 );

benches/point_source.rs

@@ -113,6 +113,7 @@ fn bench_flux_density_dipole(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 false,
+                                false,
                                 (&mut outx, &mut outy, &mut outz),
                             )
                             .unwrap(),
@@ -140,6 +141,7 @@ fn bench_flux_density_dipole(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 true,
+                                false,
                                 (&mut outx, &mut outy, &mut outz),
                             )
                             .unwrap(),
@@ -251,6 +253,7 @@ fn bench_vector_potential_dipole(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 false,
+                                false,
                                 (&mut outx, &mut outy, &mut outz),
                             )
                             .unwrap(),
@@ -278,6 +281,7 @@ fn bench_vector_potential_dipole(c: &mut Criterion) {
                                 BuildMethod::LongestAxis,
                                 HIERARCHICAL_THETA,
                                 true,
+                                false,
                                 (&mut outx, &mut outy, &mut outz),
                             )
                             .unwrap(),

cfsem/__init__.py

@@ -36,6 +36,7 @@
     gs_operator_order4,
     HierarchicalDiagnostics,
     inductance_linear_filaments,
+    inductance_linear_filaments_sparse,
     inductance_piecewise_linear_filaments,
     mutual_inductance_circular_to_linear,
     rotate_filaments_about_path,
@@ -94,6 +95,7 @@
     "solve_flux_axisymmetric",
     "filament_helix_path",
     "inductance_linear_filaments",
+    "inductance_linear_filaments_sparse",
     "inductance_matrix_axisymmetric_coaxial_rectangular_coils",
     "inductance_piecewise_linear_filaments",
     "self_inductance_piecewise_linear_filaments",

cfsem/bindings.py

@@ -12,6 +12,7 @@
 
 from numpy import asarray, ascontiguousarray, column_stack, float64, full, int64, uint64, zeros_like
 from numpy.typing import NDArray
+from scipy.sparse import csc_array, csc_matrix
 
 from cfsem.types import Array3xN
 
@@ -48,6 +49,9 @@
 from .cfsem import (
     inductance_linear_filaments_matrix as em_inductance_linear_filaments_matrix,
 )
+from .cfsem import (
+    inductance_linear_filaments_sparse_csc as em_inductance_linear_filaments_sparse_csc,
+)
 from .cfsem import (
     inductance_piecewise_linear_filaments as em_inductance_piecewise_linear_filaments,
 )
@@ -122,6 +126,7 @@
     "HierarchicalDiagnostics",
     "SolveResult",
     "inductance_linear_filaments",
+    "inductance_linear_filaments_sparse",
     "inductance_piecewise_linear_filaments",
     "mutual_inductance_circular_to_linear",
     "rotate_filaments_about_path",
@@ -1102,6 +1107,74 @@ def inductance_linear_filaments(
     raise ValueError("output must be 'vector' or 'matrix'")
 
 
+def inductance_linear_filaments_sparse(
+    xyzfil_tgt: Array3xN,
+    dlxyzfil_tgt: Array3xN,
+    xyzfil_src: Array3xN,
+    dlxyzfil_src: Array3xN,
+    interaction_map: csc_matrix | csc_array,
+    wire_radius_src: float | NDArray[float64] = 0.0,
+    par: bool = True,
+) -> csc_matrix:
+    """Evaluate selected direct source-target filament inductances.
+
+    Every stored coordinate in ``interaction_map`` is evaluated with the finite-radius source
+    vector-potential kernel and three-point Gauss--Legendre integration over the complete target
+    segment. Map data values are ignored. The result has shape ``(nsrc, ntgt)`` and exactly the
+    same CSC row-index and column-pointer arrays, including entries whose inductance is zero.
+
+    Args:
+        xyzfil_tgt: [m] target filament segment start points
+        dlxyzfil_tgt: [m] target filament segment deltas
+        xyzfil_src: [m] source filament segment start points
+        dlxyzfil_src: [m] source filament segment deltas
+        interaction_map: Canonical CSC interaction pattern with shape ``(nsrc, ntgt)``
+        wire_radius_src: [m] source filament radius, scalar or array of length ``nsrc``
+        par: Whether to evaluate stored interactions in parallel
+
+    Returns:
+        [H] CSC inductance matrix with the supplied sparsity pattern
+
+    Raises:
+        TypeError: If ``interaction_map`` is not a SciPy ``csc_matrix`` or ``csc_array``.
+        ValueError: If the map is non-canonical or has the wrong shape.
+        DimensionalityError: If filament geometry or radius lengths are inconsistent.
+    """
+    if not isinstance(interaction_map, csc_matrix | csc_array):
+        raise TypeError("interaction_map must be a scipy.sparse.csc_matrix or csc_array")
+    if not interaction_map.has_canonical_format:
+        raise ValueError("interaction_map must have sorted, unique row indices in each column")
+
+    xyzfil_tgt = _3tup_contig(xyzfil_tgt)
+    dlxyzfil_tgt = _3tup_contig(dlxyzfil_tgt)
+    xyzfil_src = _3tup_contig(xyzfil_src)
+    dlxyzfil_src = _3tup_contig(dlxyzfil_src)
+    nsrc = xyzfil_src[0].size
+    ntgt = xyzfil_tgt[0].size
+    if interaction_map.shape != (nsrc, ntgt):
+        raise ValueError(f"interaction_map must have shape ({nsrc}, {ntgt}); got {interaction_map.shape}")
+
+    if asarray(wire_radius_src).ndim == 0:
+        wire_radius_src = full(nsrc, float(wire_radius_src))
+    wire_radius_src = ascontiguousarray(wire_radius_src, dtype=float64).ravel()
+    row_indices = ascontiguousarray(interaction_map.indices, dtype=uint64)
+    column_pointers = ascontiguousarray(interaction_map.indptr, dtype=uint64)
+    values = em_inductance_linear_filaments_sparse_csc(
+        xyzfil_tgt,
+        dlxyzfil_tgt,
+        xyzfil_src,
+        dlxyzfil_src,
+        wire_radius_src,
+        row_indices,
+        column_pointers,
+        par,
+    )
+    return csc_matrix(
+        (values, interaction_map.indices.copy(), interaction_map.indptr.copy()),
+        shape=interaction_map.shape,
+    )
+
+
 def gs_operator_order2(rs: NDArray[float64], zs: NDArray[float64]) -> SparseTriplet:
     """Build second-order Grad-Shafranov operator in triplet format.
     Assumes regular grid spacing.