diff --git a/CHANGELOG.md b/CHANGELOG.md index 0c24b92..265496b 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Cargo.lock b/Cargo.lock index 69b5481..1639f71 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -128,7 +128,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "cfsem" -version = "12.1.0" +version = "13.0.0" dependencies = [ "criterion", "faer", diff --git a/Cargo.toml b/Cargo.toml index 2a77c1d..8fff106 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "12.1.0" +version = "13.0.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" diff --git a/benches/linear_filament.rs b/benches/linear_filament.rs index 21931cf..5a7d505 100644 --- a/benches/linear_filament.rs +++ b/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 ); diff --git a/benches/point_source.rs b/benches/point_source.rs index f9136e9..5e5ef15 100644 --- a/benches/point_source.rs +++ b/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(), diff --git a/cfsem/__init__.py b/cfsem/__init__.py index f0ab930..6c32912 100644 --- a/cfsem/__init__.py +++ b/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", diff --git a/cfsem/bindings.py b/cfsem/bindings.py index dd13ae4..424a82c 100644 --- a/cfsem/bindings.py +++ b/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. diff --git a/cfsem/cfsem.pyi b/cfsem/cfsem.pyi index 1c9e658..090d026 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -1,7 +1,8 @@ -from typing import TypeAlias, TypedDict +from typing import Literal, TypeAlias, TypedDict from numpy import complex128, float32, float64, int64, uint64 from numpy.typing import NDArray +from scipy.sparse import csc_matrix FloatArray: TypeAlias = NDArray[float64] ComplexArray: TypeAlias = NDArray[complex128] @@ -48,6 +49,8 @@ class HierarchicalDiagnostics: def source_tree(self) -> SourceTreeDiagnostics | None: ... @property def accepted_levels(self) -> FloatArray | None: ... + @property + def near_field_interaction_map(self) -> csc_matrix | None: ... class SolveResult: """Field arrays and diagnostics returned by a hierarchical solve.""" @@ -292,6 +295,7 @@ def flux_density_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical magnetic flux density of dipoles in Cartesian coordinates. @@ -312,6 +316,8 @@ def flux_density_dipole_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -329,6 +335,7 @@ def vector_potential_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical magnetic vector potential of dipoles in Cartesian coordinates. @@ -349,6 +356,8 @@ def vector_potential_dipole_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -375,6 +384,7 @@ def flux_density_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for many linear filament segments. @@ -396,6 +406,8 @@ def flux_density_linear_filament_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -437,6 +449,7 @@ def vector_potential_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for many linear filament segments. @@ -458,6 +471,8 @@ def vector_potential_linear_filament_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -502,6 +517,16 @@ def inductance_linear_filaments_matrix( wire_radius_src: FloatArray, par: bool = True, ) -> FloatArray: ... +def inductance_linear_filaments_sparse_csc( + xyzfil_tgt: ArrayTriple, + dlxyzfil_tgt: ArrayTriple, + xyzfil_src: ArrayTriple, + dlxyzfil_src: ArrayTriple, + wire_radius_src: FloatArray, + row_indices: UIntArray, + column_pointers: UIntArray, + par: bool = True, +) -> FloatArray: ... def gs_operator_order2(rs: FloatArray, zs: FloatArray) -> tuple[FloatArray, UIntArray, UIntArray]: ... def gs_operator_order4(rs: FloatArray, zs: FloatArray) -> tuple[FloatArray, UIntArray, UIntArray]: ... def flux_density_triangle_mesh( @@ -528,6 +553,7 @@ def flux_density_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for a triangle mesh with nodal stream-function values. @@ -557,6 +583,8 @@ def flux_density_triangle_mesh_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -574,6 +602,7 @@ def vector_potential_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for a triangle mesh with nodal stream-function values. @@ -603,6 +632,8 @@ def vector_potential_triangle_mesh_hierarchical( par: Whether to evaluate target batches in parallel. out: Optional contiguous and aligned output component arrays to fill. extra_diagnostics: Whether to populate source-tree diagnostics that require extra data. + skip: Interaction class to omit. `"near"` returns far-only, `"far"` returns + near-only, `"both"` returns zero field arrays, and `None` returns both contributions. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in diff --git a/docs/python/boundary_element.md b/docs/python/boundary_element.md index f5cbbf2..2ab9f49 100644 --- a/docs/python/boundary_element.md +++ b/docs/python/boundary_element.md @@ -1,5 +1,11 @@ # Boundary Element +Hierarchical boundary-element solvers accept `skip="near"` for a far-only result and `skip="far"` +for a direct near-only result; `skip=None` evaluates both. `skip="both"` returns zero field arrays +without a field traversal. Extra diagnostics include the resulting direct-interaction pattern as a +canonical `(ntri, ntgt)` SciPy CSC matrix; when both interactions are skipped, only the diagnostic +traversal runs. + ## Fields Triangle B-field evaluation is analytic away from each finite source triangle. At a diff --git a/docs/python/dipole.md b/docs/python/dipole.md index 1b39679..a7b543e 100644 --- a/docs/python/dipole.md +++ b/docs/python/dipole.md @@ -1,5 +1,11 @@ # Dipole +Hierarchical dipole solvers accept an optional interaction filter. `skip="near"` returns the +far-only result, `skip="far"` returns the direct near-only result, `skip="both"` returns zero field +arrays without a field traversal, and the default `skip=None` evaluates both. With +`extra_diagnostics=True`, the diagnostics include the direct-interaction pattern as a canonical +`(nsrc, ntgt)` SciPy CSC matrix; for `skip="both"`, only that diagnostic traversal runs. + ## Fields ::: cfsem.flux_density_dipole diff --git a/docs/python/filament.md b/docs/python/filament.md index a9a03db..65d3362 100644 --- a/docs/python/filament.md +++ b/docs/python/filament.md @@ -1,5 +1,13 @@ # Linear Filament +Hierarchical field solvers accept `skip=None`, `skip="near"`, `skip="far"`, or `skip="both"`. The +value names the interaction class to omit: `skip="near"` evaluates only accepted far-field +summaries, while `skip="far"` evaluates only direct near-field interactions. `skip="both"` returns +zero field arrays without a field traversal; with `extra_diagnostics=True`, only the diagnostic tree +walk runs. A full solve is the sum of the near-only and far-only solves, up to floating-point +roundoff. Filtering happens in the common evaluator, so skipped kernel calculations are not +performed. + ## Fields ::: cfsem.flux_density_linear_filament @@ -10,6 +18,63 @@ ::: cfsem.vector_potential_linear_filament_hierarchical +When `extra_diagnostics=True`, the result includes +`diagnostics.near_field_interaction_map`, a canonical SciPy CSC matrix with shape `(nsrc, ntgt)`. +Source indices are rows and target indices are columns. The map records the direct-interaction +classification for the target points supplied to the hierarchical solve; its stored data values are +structural markers only. + +## Inductance + +::: cfsem.inductance_linear_filaments + +::: cfsem.inductance_linear_filaments_sparse + +The sparse method evaluates every stored coordinate with three-point Gauss--Legendre integration +over the complete target segment. Thus, midpoint targets can classify near interactions cheaply, +while the resulting inductance entries still integrate along the full target segments: + +```python +import cfsem + +midpoints = tuple(start + 0.5 * delta for start, delta in zip(xyzfil, dlxyzfil)) +far = cfsem.vector_potential_linear_filament_hierarchical( + midpoints, + xyzfil, + dlxyzfil, + current, + wire_radius, + theta=0.05, + skip="near", # far-only field + extra_diagnostics=True, +) +near_pattern = far.diagnostics.near_field_interaction_map +near_inductance = cfsem.inductance_linear_filaments_sparse( + xyzfil, + dlxyzfil, + xyzfil, + dlxyzfil, + near_pattern, + wire_radius_src=wire_radius, +) +near_coupling = near_inductance.T @ current +``` + +The pointwise near vector potential and segment-integrated near inductance share a classification +pattern, but they are not numerically interchangeable. The map must be rebuilt if geometry, +`theta`, construction method, kernel acceptance logic, or acceptance-relevant source moments +change. For a self-coupled workflow that needs symmetric structural support, make that modeling +choice explicitly before evaluation: + +```python +symmetric_pattern = (near_pattern + near_pattern.T).astype(bool).tocsc() +symmetric_pattern.sum_duplicates() +symmetric_pattern.sort_indices() +near_inductance = cfsem.inductance_linear_filaments_sparse( + xyzfil, dlxyzfil, xyzfil, dlxyzfil, symmetric_pattern, wire_radius +) +``` + ## Force ::: cfsem.body_force_density_linear_filament diff --git a/src/physics/hierarchical/convenience.rs b/src/physics/hierarchical/convenience.rs index 52adf7a..620af71 100644 --- a/src/physics/hierarchical/convenience.rs +++ b/src/physics/hierarchical/convenience.rs @@ -15,9 +15,10 @@ use super::kernels::{ LinearFilamentSources, LinearFilamentVectorPotentialKernel, }; use super::{ - BuildMethod, ClusterTree, EvaluationScratch, HierarchicalError, HierarchicalKernel, Scalar, - SourceCollection, SourceMomentCollection, SourceNodeSummaries, TargetCollection, eval, - eval_par, scratch_len, scratch_len_par, update_summaries, + BuildMethod, ClusterTree, EvaluationScratch, HierarchicalError, HierarchicalKernel, + NearFieldInteractionMap, Scalar, Skip, SourceCollection, SourceMomentCollection, + SourceNodeSummaries, TargetCollection, TraversalDiagnostics, eval, eval_par, scratch_len, + scratch_len_par, traversal_diagnostics, traversal_diagnostics_par, update_summaries, }; /// Diagnostic information returned by stateless hierarchical solves. @@ -32,6 +33,8 @@ pub struct Diagnostics { pub source_count: usize, /// Number of targets in the solve. pub target_count: usize, + /// Accepted levels and direct near-field sparsity, when requested. + traversal_diagnostics: Option>, } impl Diagnostics { @@ -40,6 +43,33 @@ impl Diagnostics { pub fn source_tree(&self) -> &ClusterTree { &self.source_tree } + + /// Borrow the requested traversal diagnostics, if they were collected. + #[inline] + pub fn traversal_diagnostics(&self) -> Option<&TraversalDiagnostics> { + self.traversal_diagnostics.as_ref() + } + + /// Borrow the mean accepted source-tree level per target, if requested. + #[inline] + pub fn accepted_levels(&self) -> Option<&[K::Scalar]> { + self.traversal_diagnostics() + .map(|diagnostics| diagnostics.accepted_levels.as_slice()) + } + + /// Borrow the direct near-field interaction pattern, if requested. + #[inline] + pub fn near_field_interaction_map(&self) -> Option<&NearFieldInteractionMap> { + self.traversal_diagnostics() + .map(|diagnostics| &diagnostics.near_field_interaction_map) + } + + /// Move traversal diagnostics into an in-crate result adapter. + #[cfg(feature = "python")] + #[inline] + pub(crate) fn take_traversal_diagnostics(&mut self) -> Option> { + self.traversal_diagnostics.take() + } } /// Hierarchical magnetic flux density of dipole sources at Cartesian targets. @@ -61,10 +91,11 @@ impl Diagnostics { /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, tree construction fails, @@ -77,6 +108,7 @@ pub fn flux_density_dipole_hierarchical( construction_method: BuildMethod, theta: T, par: bool, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result>, HierarchicalError> { let sources = DipoleSources::new(loc.0, loc.1, loc.2, outer_radius); @@ -90,6 +122,8 @@ pub fn flux_density_dipole_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, out, ) } @@ -113,10 +147,11 @@ pub fn flux_density_dipole_hierarchical( /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, tree construction fails, @@ -129,6 +164,7 @@ pub fn vector_potential_dipole_hierarchical( construction_method: BuildMethod, theta: T, par: bool, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result>, HierarchicalError> { let sources = DipoleSources::new(loc.0, loc.1, loc.2, outer_radius); @@ -142,6 +178,8 @@ pub fn vector_potential_dipole_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, out, ) } @@ -166,10 +204,11 @@ pub fn vector_potential_dipole_hierarchical( /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, tree construction fails, @@ -183,6 +222,7 @@ pub fn flux_density_linear_filament_hierarchical( construction_method: BuildMethod, theta: T, par: bool, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result>, HierarchicalError> { let sources = LinearFilamentSources::new(xyzfil, dlxyzfil, wire_radius); @@ -195,6 +235,8 @@ pub fn flux_density_linear_filament_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, out, ) } @@ -219,10 +261,11 @@ pub fn flux_density_linear_filament_hierarchical( /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, tree construction fails, @@ -236,6 +279,7 @@ pub fn vector_potential_linear_filament_hierarchical( construction_method: BuildMethod, theta: T, par: bool, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result>, HierarchicalError> { let sources = LinearFilamentSources::new(xyzfil, dlxyzfil, wire_radius); @@ -248,6 +292,8 @@ pub fn vector_potential_linear_filament_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, out, ) } @@ -278,10 +324,11 @@ pub fn vector_potential_linear_filament_hierarchical( /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, mesh conversion fails, @@ -293,6 +340,7 @@ pub fn flux_density_triangle_mesh_hierarchical( construction_method: BuildMethod, theta: f64, par: bool, + extra_diagnostics: bool, out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result>, HierarchicalError> { mesh.validate_nodal_values(s) @@ -308,6 +356,8 @@ pub fn flux_density_triangle_mesh_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, out, ) } @@ -338,10 +388,11 @@ pub fn flux_density_triangle_mesh_hierarchical( /// construction_method: Source-tree construction method. /// theta: Barnes-Hut acceptance angle. Smaller values are more accurate and slower. /// par: Whether to evaluate target batches in parallel. +/// extra_diagnostics: Whether to collect accepted levels and the near-field interaction map. /// out: Output component slices to fill. /// /// Returns: -/// Source-tree diagnostics and construction/evaluation timing on success. +/// Source-tree metadata, optional traversal diagnostics, and timing on success. /// /// Errors: /// Returns [`HierarchicalError`] when input lengths are inconsistent, mesh conversion fails, @@ -353,6 +404,7 @@ pub fn vector_potential_triangle_mesh_hierarchical( construction_method: BuildMethod, theta: f64, par: bool, + extra_diagnostics: bool, out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result>, HierarchicalError> { mesh.validate_nodal_values(s) @@ -368,6 +420,190 @@ pub fn vector_potential_triangle_mesh_hierarchical( construction_method, theta, par, + None, + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`flux_density_dipole_hierarchical`]. +pub fn flux_density_dipole_hierarchical_with_skip( + loc: (&[T], &[T], &[T]), + moment: (&[T], &[T], &[T]), + obs: (&[T], &[T], &[T]), + outer_radius: &[T], + construction_method: BuildMethod, + theta: T, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [T], &mut [T], &mut [T]), +) -> Result>, HierarchicalError> { + let sources = DipoleSources::new(loc.0, loc.1, loc.2, outer_radius); + let moments = DipoleMoments::new(moment.0, moment.1, moment.2); + let targets = DipoleTargets::new(obs.0, obs.1, obs.2); + one_shot_vec3( + DipoleFluxDensityKernel::::new(), + sources, + moments, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`vector_potential_dipole_hierarchical`]. +pub fn vector_potential_dipole_hierarchical_with_skip( + loc: (&[T], &[T], &[T]), + moment: (&[T], &[T], &[T]), + obs: (&[T], &[T], &[T]), + outer_radius: &[T], + construction_method: BuildMethod, + theta: T, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [T], &mut [T], &mut [T]), +) -> Result>, HierarchicalError> { + let sources = DipoleSources::new(loc.0, loc.1, loc.2, outer_radius); + let moments = DipoleMoments::new(moment.0, moment.1, moment.2); + let targets = DipoleTargets::new(obs.0, obs.1, obs.2); + one_shot_vec3( + DipoleVectorPotentialKernel::::new(), + sources, + moments, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`flux_density_linear_filament_hierarchical`]. +pub fn flux_density_linear_filament_hierarchical_with_skip( + xyzp: (&[T], &[T], &[T]), + xyzfil: (&[T], &[T], &[T]), + dlxyzfil: (&[T], &[T], &[T]), + ifil: &[T], + wire_radius: &[T], + construction_method: BuildMethod, + theta: T, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [T], &mut [T], &mut [T]), +) -> Result>, HierarchicalError> { + let sources = LinearFilamentSources::new(xyzfil, dlxyzfil, wire_radius); + let targets = DipoleTargets::new(xyzp.0, xyzp.1, xyzp.2); + one_shot_vec3( + LinearFilamentFluxDensityKernel::::new(), + sources, + ifil, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`vector_potential_linear_filament_hierarchical`]. +pub fn vector_potential_linear_filament_hierarchical_with_skip( + xyzp: (&[T], &[T], &[T]), + xyzfil: (&[T], &[T], &[T]), + dlxyzfil: (&[T], &[T], &[T]), + ifil: &[T], + wire_radius: &[T], + construction_method: BuildMethod, + theta: T, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [T], &mut [T], &mut [T]), +) -> Result>, HierarchicalError> { + let sources = LinearFilamentSources::new(xyzfil, dlxyzfil, wire_radius); + let targets = DipoleTargets::new(xyzp.0, xyzp.1, xyzp.2); + one_shot_vec3( + LinearFilamentVectorPotentialKernel::::new(), + sources, + ifil, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`flux_density_triangle_mesh_hierarchical`]. +pub fn flux_density_triangle_mesh_hierarchical_with_skip( + obs: (&[f64], &[f64], &[f64]), + mesh: &TriangleMeshView<'_>, + s: &[f64], + construction_method: BuildMethod, + theta: f64, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result>, HierarchicalError> { + mesh.validate_nodal_values(s) + .map_err(|_| HierarchicalError::LengthMismatch)?; + let sources = BoundaryElementTriangles::new(mesh); + let moments = BoundaryElementNodalValues::new(sources, s); + let targets = DipoleTargets::new(obs.0, obs.1, obs.2); + one_shot_vec3( + BoundaryElementFluxDensityKernel::::new(), + sources, + moments, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, + out, + ) +} + +/// Filtered variant of [`vector_potential_triangle_mesh_hierarchical`]. +pub fn vector_potential_triangle_mesh_hierarchical_with_skip( + obs: (&[f64], &[f64], &[f64]), + mesh: &TriangleMeshView<'_>, + s: &[f64], + construction_method: BuildMethod, + theta: f64, + par: bool, + skip: Skip, + extra_diagnostics: bool, + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result>, HierarchicalError> { + mesh.validate_nodal_values(s) + .map_err(|_| HierarchicalError::LengthMismatch)?; + let sources = BoundaryElementTriangles::new(mesh); + let moments = BoundaryElementNodalValues::new(sources, s); + let targets = DipoleTargets::new(obs.0, obs.1, obs.2); + one_shot_vec3( + BoundaryElementVectorPotentialKernel::::new(), + sources, + moments, + targets, + construction_method, + theta, + par, + Some(skip), + extra_diagnostics, out, ) } @@ -381,6 +617,8 @@ pub(crate) fn one_shot_vec3( construction_method: BuildMethod, theta: T, par: bool, + skip: Option, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result, HierarchicalError> where @@ -407,63 +645,130 @@ where BuildMethod::LongestAxis => ClusterTree::build(sources)?, BuildMethod::MortonLbvh => ClusterTree::build_morton_lbvh(sources)?, }; - let mut source_summaries = SourceNodeSummaries::::new(source_tree.as_view()); - let mut err = update_summaries( - &kernel, - source_tree.as_view(), - sources, - moments, - &mut source_summaries.node_summaries, - ); - if err != HierarchicalError::Ok { - return Err(err); - } - let construction_seconds = construction_start.elapsed().as_secs_f64(); - - let evaluation_start = Instant::now(); - let scratch_len = match par { - true => scratch_len_par(targets.len()), - false => scratch_len(), - }; - let mut scratch_values = vec![[T::ZERO; 3]; scratch_len]; - let mut scratch = EvaluationScratch { - contribution: &mut scratch_values, - }; - let out_components = [out.0, out.1, out.2]; - err = match par { - true => eval_par( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - sources, - targets, - moments, - theta, - out_components, - &mut scratch, - ), - false => eval( + let source_summaries = if skip == Some(Skip::Both) && !extra_diagnostics { + None + } else { + let mut source_summaries = SourceNodeSummaries::::new(source_tree.as_view()); + let err = update_summaries( &kernel, source_tree.as_view(), - &source_summaries.node_summaries, sources, - targets, moments, - theta, - out_components, - &mut scratch, - ), + &mut source_summaries.node_summaries, + ); + if err != HierarchicalError::Ok { + return Err(err); + } + Some(source_summaries) }; - if err != HierarchicalError::Ok { - return Err(err); + let construction_seconds = construction_start.elapsed().as_secs_f64(); + + let evaluation_start = Instant::now(); + if skip != Some(Skip::Both) { + let source_summaries = source_summaries + .as_ref() + .expect("non-skipped evaluation requires source summaries"); + let scratch_len = match par { + true => scratch_len_par(targets.len()), + false => scratch_len(), + }; + let mut scratch_values = vec![[T::ZERO; 3]; scratch_len]; + let mut scratch = EvaluationScratch { + contribution: &mut scratch_values, + }; + let out_components = [out.0, out.1, out.2]; + // Keep the filter variant visible at this boundary so inlining can remove the terminal-node + // filter checks from unfiltered solves while the evaluator API remains consolidated. + let err = match (par, skip) { + (true, Some(skip)) => eval_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + Some(skip), + out_components, + &mut scratch, + ), + (true, None) => eval_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + None, + out_components, + &mut scratch, + ), + (false, Some(skip)) => eval( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + Some(skip), + out_components, + &mut scratch, + ), + (false, None) => eval( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + None, + out_components, + &mut scratch, + ), + }; + if err != HierarchicalError::Ok { + return Err(err); + } + } else { + out.0.fill(T::ZERO); + out.1.fill(T::ZERO); + out.2.fill(T::ZERO); } let evaluation_seconds = evaluation_start.elapsed().as_secs_f64(); + let traversal_diagnostics = if extra_diagnostics { + let source_summaries = source_summaries + .as_ref() + .expect("requested traversal diagnostics require source summaries"); + Some(match par { + true => traversal_diagnostics_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + targets, + theta, + )?, + false => traversal_diagnostics( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + targets, + theta, + )?, + }) + } else { + None + }; + Ok(Diagnostics { source_tree, construction_seconds, evaluation_seconds, source_count: sources.len(), target_count: targets.len(), + traversal_diagnostics, }) } diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index 88b0557..a20e3c3 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -1,7 +1,8 @@ use super::{ - BoundedGeometry, ClusterTreeView, HierarchicalError, HierarchicalKernel, Scalar, + BoundedGeometry, ClusterTreeView, HierarchicalError, HierarchicalKernel, Scalar, Skip, SourceCollection, SourceMomentCollection, TargetCollection, }; +use rayon::prelude::*; use std::sync::atomic::{AtomicU32, Ordering}; /// CPU-owned source summary storage. @@ -90,8 +91,11 @@ where /// /// This is the public hierarchical evaluation path. Each target is summarized /// as a single target leaf, walked against the source tree, and written directly -/// into caller-provided component slices. The output slice count must match the -/// kernel output dimension `D`. +/// into caller-provided component slices. [`Skip::Near`] retains accepted far-summary +/// contributions only, [`Skip::Far`] retains direct leaf contributions only, and `None` +/// evaluates both interaction classes. [`Skip::Both`] zeroes the output without target +/// summarization, source-tree traversal, or contribution scratch. The output slice count must +/// match the kernel output dimension `D`. #[inline] pub fn eval( kernel: &K, @@ -101,6 +105,7 @@ pub fn eval( targets: C, moments: M, theta: T, + skip: Option, out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError @@ -116,22 +121,66 @@ where if err != HierarchicalError::Ok { return err; } - eval_validated( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - out, - scratch, - ) + match skip { + None => eval_validated::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + scratch, + ), + Some(Skip::Near) => eval_validated::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + scratch, + ), + Some(Skip::Far) => eval_validated::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + scratch, + ), + Some(Skip::Both) => eval_validated::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + scratch, + ), + } } #[inline] /// Evaluate validated source-target rows with the hierarchical tree walk. -fn eval_validated( +fn eval_validated< + K, + T, + S, + M, + C, + const D: usize, + const EVALUATE_NEAR: bool, + const EVALUATE_FAR: bool, +>( kernel: &K, source_tree: ClusterTreeView<'_, T>, source_summaries: &[K::SourceSummary], @@ -164,6 +213,12 @@ where return HierarchicalError::LengthMismatch; } } + if !EVALUATE_NEAR && !EVALUATE_FAR { + for component in out { + component.fill(T::ZERO); + } + return HierarchicalError::Ok; + } if source_summaries.len() < source_tree.n_nodes() || scratch.contribution.is_empty() { return HierarchicalError::ScratchTooSmall; } @@ -175,7 +230,7 @@ where for target_id in 0..targets.len() { let target = targets.target(target_id); - let err = eval_scalar( + let err = eval_scalar::( kernel, source_tree, source_summaries, @@ -200,12 +255,121 @@ where HierarchicalError::Ok } +/// Handle terminal nodes selected by the shared source-tree traversal. +trait TraversalVisitor { + /// Handle a source node accepted through the kernel's far criterion. + fn on_far_accept( + &mut self, + source_node_index: usize, + source_level: u32, + source_summary: &K::SourceSummary, + ); + + /// Handle a rejected source leaf through direct source interactions. + fn on_near_leaf(&mut self, source_node_index: usize, source_level: u32, source_ids: &[u32]); +} + +/// Traverse one target against the source tree and report each terminal node. +/// +/// Field evaluation and traversal diagnostics both use this function, keeping +/// kernel-specific acceptance and leaf fallback behavior identical. +#[inline] +fn traverse_source_tree( + kernel: &K, + source_tree: ClusterTreeView<'_, K::Scalar>, + source_summaries: &[K::SourceSummary], + target: &K::TargetGeometry, + theta: K::Scalar, + active: &mut Vec<(u32, u32)>, + visitor: &mut V, +) where + K: HierarchicalKernel, + V: TraversalVisitor, +{ + active.clear(); + active.push((0_u32, 0_u32)); + let target_aabb = target.aabb(); + while let Some((source_node, source_level)) = active.pop() { + let source_node_index = source_node as usize; + let source_summary = &source_summaries[source_node_index]; + let source_aabb = source_tree.node_aabb[source_node_index]; + if kernel.accept_far(target_aabb, source_aabb, source_summary, theta) { + visitor.on_far_accept(source_node_index, source_level, source_summary); + continue; + } + + let leaf_count = source_tree.leaf_count[source_node_index]; + if leaf_count > 0 { + let start = source_tree.leaf_start[source_node_index] as usize; + let count = leaf_count as usize; + let end = start + count; + let source_ids = &source_tree.sorted_indices[start..end]; + visitor.on_near_leaf(source_node_index, source_level, source_ids); + } else { + let next_level = source_level + 1; + active.push((source_tree.node_left_child[source_node_index], next_level)); + active.push((source_tree.node_right_child[source_node_index], next_level)); + } + } +} + +/// Terminal-node visitor that evaluates far summaries and direct leaf sources. +struct EvaluationTraversalVisitor<'a, K, S, M, const EVALUATE_NEAR: bool, const EVALUATE_FAR: bool> +where + K: HierarchicalKernel, +{ + kernel: &'a K, + sources: S, + target: K::TargetGeometry, + moments: M, + out: &'a mut K::Output, + contribution: &'a mut K::Output, + target_summary: &'a K::TargetSummary, +} + +impl TraversalVisitor + for EvaluationTraversalVisitor<'_, K, S, M, EVALUATE_NEAR, EVALUATE_FAR> +where + K: HierarchicalKernel, + S: SourceCollection, + M: SourceMomentCollection, +{ + #[inline] + fn on_far_accept( + &mut self, + _source_node_index: usize, + _source_level: u32, + source_summary: &K::SourceSummary, + ) { + if EVALUATE_FAR { + self.kernel + .eval_far(self.target_summary, source_summary, self.contribution); + self.kernel.accumulate(self.out, self.contribution); + } + } + + #[inline] + fn on_near_leaf(&mut self, _source_node_index: usize, _source_level: u32, source_ids: &[u32]) { + if !EVALUATE_NEAR { + return; + } + for &source_id in source_ids { + let source_id = source_id as usize; + let source = self.sources.source(source_id); + let moment = self.moments.moment(source_id); + self.kernel + .eval_near(&self.target, &source, &moment, self.contribution); + self.kernel.accumulate(self.out, self.contribution); + } + } +} + /// Evaluate one scalar target against the source tree. /// -/// Serial and parallel vector evaluators both call this helper so the source -/// traversal and acceptance behavior cannot diverge between evaluation modes. +/// Serial and parallel vector evaluators both call this helper, and its +/// terminal-node actions use the same traversal as diagnostics. #[inline] -fn eval_scalar( +fn eval_scalar( kernel: &K, source_tree: ClusterTreeView<'_, K::Scalar>, source_summaries: &[K::SourceSummary], @@ -216,7 +380,7 @@ fn eval_scalar( out: &mut K::Output, contribution: &mut K::Output, target_summary: &mut K::TargetSummary, - active: &mut Vec, + active: &mut Vec<(u32, u32)>, target_ids: &[u32], ) -> HierarchicalError where @@ -232,36 +396,24 @@ where return err; } - active.clear(); - active.push(0_u32); - while let Some(source_node) = active.pop() { - let source_node_index = source_node as usize; - let source_summary = &source_summaries[source_node_index]; - let source_aabb = source_tree.node_aabb[source_node_index]; - if kernel.accept_far(target.aabb(), source_aabb, source_summary, theta) { - kernel.eval_far(target_summary, source_summary, contribution); - kernel.accumulate(out, contribution); - continue; - } - - let leaf_count = source_tree.leaf_count[source_node_index]; - if leaf_count > 0 { - let start = source_tree.leaf_start[source_node_index] as usize; - let count = leaf_count as usize; - let end = start + count; - let source_ids = &source_tree.sorted_indices[start..end]; - for i in 0..source_ids.len() { - let source_id = source_ids[i] as usize; - let source = sources.source(source_id); - let moment = moments.moment(source_id); - kernel.eval_near(&target, &source, &moment, contribution); - kernel.accumulate(out, contribution); - } - } else { - active.push(source_tree.node_left_child[source_node_index]); - active.push(source_tree.node_right_child[source_node_index]); - } - } + let mut visitor = EvaluationTraversalVisitor:: { + kernel, + sources, + target, + moments, + out, + contribution, + target_summary, + }; + traverse_source_tree( + kernel, + source_tree, + source_summaries, + &target, + theta, + active, + &mut visitor, + ); HierarchicalError::Ok } @@ -271,8 +423,10 @@ where /// This is intentionally the simplest parallelization of the single-tree /// solver: each worker owns disjoint target and component output slices and /// runs the serial source-tree evaluator on that slice. It shares the source -/// tree and source summaries between workers, and avoids any cross-thread -/// output accumulation. +/// tree and source summaries between workers, and avoids any cross-thread output accumulation. +/// [`Skip::Near`] retains accepted far-summary contributions only, [`Skip::Far`] retains direct +/// leaf contributions only, and `None` evaluates both interaction classes. [`Skip::Both`] zeroes +/// the output without target summarization, source-tree traversal, or parallel scratch use. #[inline] pub fn eval_par( kernel: &K, @@ -282,6 +436,7 @@ pub fn eval_par( targets: C, moments: M, theta: T, + skip: Option, out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError @@ -311,6 +466,12 @@ where return HierarchicalError::LengthMismatch; } } + if skip == Some(Skip::Both) { + for component in out { + component.fill(T::ZERO); + } + return HierarchicalError::Ok; + } if source_summaries.len() < source_tree.n_nodes() { return HierarchicalError::ScratchTooSmall; } @@ -325,26 +486,64 @@ where } let error_code = AtomicU32::new(HierarchicalError::Ok as u32); - eval_par_chunks( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - out, - &mut scratch.contribution[..chunk_count], - chunk_size, - &error_code, - ); + match skip { + None => eval_par_chunks::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + &mut scratch.contribution[..chunk_count], + chunk_size, + &error_code, + ), + Some(Skip::Near) => eval_par_chunks::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + &mut scratch.contribution[..chunk_count], + chunk_size, + &error_code, + ), + Some(Skip::Far) => eval_par_chunks::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + &mut scratch.contribution[..chunk_count], + chunk_size, + &error_code, + ), + Some(Skip::Both) => unreachable!("Skip::Both returns before parallel evaluation"), + } HierarchicalError::from_u32(error_code.load(Ordering::Relaxed)) } #[inline] /// Evaluate validated output chunks in parallel and preserve the first error code. -fn eval_par_chunks( +fn eval_par_chunks< + K, + T, + S, + M, + C, + const D: usize, + const EVALUATE_NEAR: bool, + const EVALUATE_FAR: bool, +>( kernel: &K, source_tree: ClusterTreeView<'_, T>, source_summaries: &[K::SourceSummary], @@ -373,7 +572,7 @@ fn eval_par_chunks( let mut chunk_scratch = EvaluationScratch { contribution: &mut scratch_contributions[..1], }; - let err = eval_validated( + let err = eval_validated::( kernel, source_tree, source_summaries, @@ -405,7 +604,7 @@ fn eval_par_chunks( rayon::join( || { - eval_par_chunks( + eval_par_chunks::( kernel, source_tree, source_summaries, @@ -420,7 +619,7 @@ fn eval_par_chunks( ); }, || { - eval_par_chunks( + eval_par_chunks::( kernel, source_tree, source_summaries, @@ -454,81 +653,274 @@ fn split_output_components( (left, right) } -/// Compute the source-tree level represented at each target by the terminal traversal nodes. +/// Canonical CSC sparsity for direct source-target interactions selected by a tree walk. /// -/// This is a diagnostic companion to [`eval`]. It mirrors -/// the same source-tree walk but does not evaluate field values. Far-accepted -/// nodes contribute their traversal depth, while direct leaf fallbacks -/// contribute the leaf depth. Each contribution is weighted by the number of -/// original source items represented by each terminal node, giving per-target -/// accepted levels. -#[inline] -pub fn accepted_levels( +/// Rows are original source indices and columns are target indices, so the shape is +/// `(source_count, target_count)`. Row indices are sorted within each column and are unique. +#[derive(Clone, Debug, PartialEq, Eq)] +pub struct NearFieldInteractionMap { + /// Original source indices for stored direct interactions. + pub row_indices: Vec, + /// CSC column offsets, with length `target_count + 1`. + pub column_pointers: Vec, + /// Number of source rows in the sparse pattern. + pub source_count: usize, + /// Number of target columns in the sparse pattern. + pub target_count: usize, +} + +/// Diagnostic data collected in one traversal of the source tree per target. +#[derive(Clone, Debug, PartialEq)] +pub struct TraversalDiagnostics { + /// Source-tree level represented at each target. + pub accepted_levels: Vec, + /// Direct near-field source-target interaction pattern. + pub near_field_interaction_map: NearFieldInteractionMap, +} + +/// Terminal-node visitor that collects one target's traversal diagnostics. +struct DiagnosticTraversalVisitor<'a, K> +where + K: HierarchicalKernel, +{ + source_tree: ClusterTreeView<'a, K::Scalar>, + weighted_level: K::Scalar, + represented_sources: K::Scalar, + row_indices: &'a mut Vec, +} + +impl DiagnosticTraversalVisitor<'_, K> +where + K: HierarchicalKernel, +{ + #[inline] + fn record_terminal_node(&mut self, source_node_index: usize, source_level: u32) { + let source_count = crate::math::cast::( + self.source_tree.node_range_count[source_node_index] as f64, + ); + self.weighted_level = self.weighted_level + + crate::math::cast::(f64::from(source_level)) * source_count; + self.represented_sources = self.represented_sources + source_count; + } + + #[inline] + fn accepted_level(&self) -> K::Scalar { + if self.represented_sources > K::Scalar::ZERO { + self.weighted_level / self.represented_sources + } else { + crate::math::cast::(f64::NAN) + } + } +} + +impl TraversalVisitor for DiagnosticTraversalVisitor<'_, K> +where + K: HierarchicalKernel, +{ + #[inline] + fn on_far_accept( + &mut self, + source_node_index: usize, + source_level: u32, + _source_summary: &K::SourceSummary, + ) { + self.record_terminal_node(source_node_index, source_level); + } + + #[inline] + fn on_near_leaf(&mut self, source_node_index: usize, source_level: u32, source_ids: &[u32]) { + self.record_terminal_node(source_node_index, source_level); + self.row_indices.extend_from_slice(source_ids); + } +} + +/// Traversal diagnostics collected for a contiguous target chunk. +struct TraversalDiagnosticsChunk { + accepted_levels: Vec, + row_indices: Vec, + column_lengths: Vec, +} + +/// Collect traversal diagnostics for a validated contiguous target chunk. +fn traversal_diagnostics_chunk( kernel: &K, source_tree: ClusterTreeView<'_, K::Scalar>, source_summaries: &[K::SourceSummary], targets: C, theta: K::Scalar, - out: &mut [K::Scalar], -) -> HierarchicalError +) -> TraversalDiagnosticsChunk where K: HierarchicalKernel, K::TargetGeometry: Copy, C: TargetCollection, +{ + let mut accepted_levels = Vec::with_capacity(targets.len()); + let mut row_indices = Vec::new(); + let mut column_lengths = Vec::with_capacity(targets.len()); + let mut active = Vec::new(); + + for target_id in 0..targets.len() { + let target = targets.target(target_id); + let column_start = row_indices.len(); + let mut visitor = DiagnosticTraversalVisitor:: { + source_tree, + weighted_level: K::Scalar::ZERO, + represented_sources: K::Scalar::ZERO, + row_indices: &mut row_indices, + }; + traverse_source_tree( + kernel, + source_tree, + source_summaries, + &target, + theta, + &mut active, + &mut visitor, + ); + accepted_levels.push(visitor.accepted_level()); + row_indices[column_start..].sort_unstable(); + column_lengths.push(row_indices.len() - column_start); + } + + TraversalDiagnosticsChunk { + accepted_levels, + row_indices, + column_lengths, + } +} + +/// Validate inputs shared by serial and parallel traversal diagnostics. +fn validate_traversal_diagnostics_inputs( + source_tree: ClusterTreeView<'_, K::Scalar>, + source_summaries: &[K::SourceSummary], + targets: C, +) -> Result<(), HierarchicalError> +where + K: HierarchicalKernel, + C: TargetCollection, { let err = validate_source_tree_layout(source_tree); if err != HierarchicalError::Ok { - return err; + return Err(err); } - if targets.len() != out.len() || !targets.valid_lengths() { - return HierarchicalError::LengthMismatch; + if !targets.valid_lengths() { + return Err(HierarchicalError::LengthMismatch); } if source_summaries.len() < source_tree.n_nodes() { - return HierarchicalError::ScratchTooSmall; + return Err(HierarchicalError::ScratchTooSmall); } + Ok(()) +} - let mut active = Vec::new(); - for target_id in 0..targets.len() { - let target = targets.target(target_id); - let mut weighted_level = K::Scalar::ZERO; - let mut represented_sources = K::Scalar::ZERO; - - active.clear(); - active.push((0_u32, 0_u32)); - while let Some((source_node, source_level)) = active.pop() { - let source_node_index = source_node as usize; - let source_count = crate::math::cast::( - source_tree.node_range_count[source_node_index] as f64, - ); - let source_summary = &source_summaries[source_node_index]; - let source_aabb = source_tree.node_aabb[source_node_index]; - if kernel.accept_far(target.aabb(), source_aabb, source_summary, theta) { - weighted_level = weighted_level - + crate::math::cast::(f64::from(source_level)) * source_count; - represented_sources = represented_sources + source_count; - continue; - } +/// Merge target chunks into canonical CSC traversal diagnostics. +fn merge_traversal_diagnostics_chunks( + chunks: Vec>, + source_count: usize, + target_count: usize, +) -> TraversalDiagnostics { + let entry_count = chunks.iter().map(|chunk| chunk.row_indices.len()).sum(); + let mut accepted_levels = Vec::with_capacity(target_count); + let mut row_indices = Vec::with_capacity(entry_count); + let mut column_pointers = Vec::with_capacity(target_count + 1); + column_pointers.push(0); - let leaf_count = source_tree.leaf_count[source_node_index]; - if leaf_count > 0 { - weighted_level = weighted_level - + crate::math::cast::(f64::from(source_level)) * source_count; - represented_sources = represented_sources + source_count; - } else { - let next_level = source_level + 1; - active.push((source_tree.node_left_child[source_node_index], next_level)); - active.push((source_tree.node_right_child[source_node_index], next_level)); - } + for chunk in chunks { + accepted_levels.extend(chunk.accepted_levels); + row_indices.extend(chunk.row_indices); + for column_length in chunk.column_lengths { + column_pointers.push(column_pointers.last().copied().unwrap() + column_length); } + } - out[target_id] = if represented_sources > K::Scalar::ZERO { - weighted_level / represented_sources - } else { - crate::math::cast::(f64::NAN) - }; + debug_assert_eq!(accepted_levels.len(), target_count); + debug_assert_eq!(column_pointers.len(), target_count + 1); + TraversalDiagnostics { + accepted_levels, + near_field_interaction_map: NearFieldInteractionMap { + row_indices, + column_pointers, + source_count, + target_count, + }, } +} - HierarchicalError::Ok +/// Collect accepted levels and the direct near-field CSC pattern. +/// +/// This uses the same terminal-node traversal as field evaluation. The interaction map records +/// every original source owned by a rejected terminal leaf; far-accepted nodes are omitted. The +/// resulting pattern has shape `(source_count, target_count)`. +pub fn traversal_diagnostics( + kernel: &K, + source_tree: ClusterTreeView<'_, K::Scalar>, + source_summaries: &[K::SourceSummary], + targets: C, + theta: K::Scalar, +) -> Result, HierarchicalError> +where + K: HierarchicalKernel, + K::TargetGeometry: Copy, + C: TargetCollection, +{ + validate_traversal_diagnostics_inputs::(source_tree, source_summaries, targets)?; + let source_count = source_tree.node_range_count[0] as usize; + let target_count = targets.len(); + let chunk = traversal_diagnostics_chunk(kernel, source_tree, source_summaries, targets, theta); + Ok(merge_traversal_diagnostics_chunks( + vec![chunk], + source_count, + target_count, + )) +} + +/// Collect accepted levels and the direct near-field CSC pattern in parallel over targets. +/// +/// Target chunks are traversed independently, then merged in target order to preserve canonical +/// CSC columns and byte-for-byte agreement with [`traversal_diagnostics`]. +pub fn traversal_diagnostics_par( + kernel: &K, + source_tree: ClusterTreeView<'_, K::Scalar>, + source_summaries: &[K::SourceSummary], + targets: C, + theta: K::Scalar, +) -> Result, HierarchicalError> +where + K: HierarchicalKernel + Sync, + K::TargetGeometry: Copy, + C: TargetCollection, +{ + validate_traversal_diagnostics_inputs::(source_tree, source_summaries, targets)?; + let source_count = source_tree.node_range_count[0] as usize; + let target_count = targets.len(); + if target_count == 0 { + return Ok(merge_traversal_diagnostics_chunks( + Vec::new(), + source_count, + target_count, + )); + } + + let chunk_size = crate::chunksize(target_count); + let chunk_count = target_count.div_ceil(chunk_size); + let chunks = (0..chunk_count) + .into_par_iter() + .map(|chunk_id| { + let start = chunk_id * chunk_size; + let end = (start + chunk_size).min(target_count); + traversal_diagnostics_chunk( + kernel, + source_tree, + source_summaries, + targets.slice(start, end), + theta, + ) + }) + .collect(); + Ok(merge_traversal_diagnostics_chunks( + chunks, + source_count, + target_count, + )) } /// Dense exact fallback using nested range loops. diff --git a/src/physics/hierarchical/kernel.rs b/src/physics/hierarchical/kernel.rs index 51976d0..2c60aa6 100644 --- a/src/physics/hierarchical/kernel.rs +++ b/src/physics/hierarchical/kernel.rs @@ -1,5 +1,16 @@ use super::{Aabb, Scalar}; +/// Hierarchical interaction class to omit during a filtered evaluation. +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +pub enum Skip { + /// Omit exact source interactions reached through direct leaves. + Near, + /// Omit accepted source-summary interactions. + Far, + /// Omit both interaction classes, producing zero field output without a tree walk. + Both, +} + /// Runtime error code for hierarchical tree operations. #[repr(u32)] #[derive(Clone, Copy, Debug, PartialEq, Eq)] diff --git a/src/physics/hierarchical/mod.rs b/src/physics/hierarchical/mod.rs index a913ad2..0f286fb 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -21,9 +21,14 @@ pub mod kernels; pub mod tree; pub use convenience::{ - flux_density_dipole_hierarchical, flux_density_linear_filament_hierarchical, - flux_density_triangle_mesh_hierarchical, vector_potential_dipole_hierarchical, - vector_potential_linear_filament_hierarchical, vector_potential_triangle_mesh_hierarchical, + Diagnostics, flux_density_dipole_hierarchical, flux_density_dipole_hierarchical_with_skip, + flux_density_linear_filament_hierarchical, flux_density_linear_filament_hierarchical_with_skip, + flux_density_triangle_mesh_hierarchical, flux_density_triangle_mesh_hierarchical_with_skip, + vector_potential_dipole_hierarchical, vector_potential_dipole_hierarchical_with_skip, + vector_potential_linear_filament_hierarchical, + vector_potential_linear_filament_hierarchical_with_skip, + vector_potential_triangle_mesh_hierarchical, + vector_potential_triangle_mesh_hierarchical_with_skip, }; pub(crate) use crate::math::Scalar; @@ -31,9 +36,13 @@ pub(crate) use aabb::Aabb; #[cfg(test)] pub(crate) use evaluator::eval_dense; pub(crate) use evaluator::{ - EvaluationScratch, SourceNodeSummaries, eval, eval_par, scratch_len, scratch_len_par, - update_summaries, + EvaluationScratch, SourceNodeSummaries, scratch_len, scratch_len_par, update_summaries, }; +pub use evaluator::{ + NearFieldInteractionMap, TraversalDiagnostics, eval, eval_par, traversal_diagnostics, + traversal_diagnostics_par, +}; +pub use kernel::Skip; pub(crate) use kernel::{ BoundedGeometry, BoundedGeometryCollection, HierarchicalError, HierarchicalKernel, SourceCollection, SourceMomentCollection, TargetCollection, geometric_accept_far, diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index 15435f7..58047a7 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -10,6 +10,7 @@ use crate::physics::hierarchical::kernels::{ LinearFilamentVectorPotentialSummary, }; use crate::physics::point_source::segment::flux_density_point_segment_scalar; +use std::sync::atomic::{AtomicUsize, Ordering}; #[derive(Clone, Copy)] struct MockPoint { @@ -43,12 +44,20 @@ struct TargetSummary { struct MockKernel { _marker: core::marker::PhantomData, + target_summary_calls: AtomicUsize, + accept_calls: AtomicUsize, + near_calls: AtomicUsize, + far_calls: AtomicUsize, } impl MockKernel { fn new() -> Self { Self { _marker: core::marker::PhantomData, + target_summary_calls: AtomicUsize::new(0), + accept_calls: AtomicUsize::new(0), + near_calls: AtomicUsize::new(0), + far_calls: AtomicUsize::new(0), } } } @@ -119,6 +128,7 @@ impl HierarchicalKernel for MockKernel { targets: &[Self::TargetGeometry], out: &mut Self::TargetSummary, ) -> HierarchicalError { + self.target_summary_calls.fetch_add(1, Ordering::Relaxed); *out = TargetSummary::default(); for i in 0..target_ids.len() { let id = target_ids[i] as usize; @@ -135,6 +145,17 @@ impl HierarchicalKernel for MockKernel { HierarchicalError::Ok } + fn accept_far( + &self, + target_aabb: Aabb, + source_aabb: Aabb, + _source: &Self::SourceSummary, + theta: T, + ) -> bool { + self.accept_calls.fetch_add(1, Ordering::Relaxed); + geometric_accept_far(target_aabb, source_aabb, theta) + } + fn eval_near( &self, target: &Self::TargetGeometry, @@ -142,6 +163,7 @@ impl HierarchicalKernel for MockKernel { moment: &Self::SourceMoment, out: &mut Self::Output, ) { + self.near_calls.fetch_add(1, Ordering::Relaxed); let r2 = dist2(target.point, source.point); out[0] = *moment / (T::ONE + r2); } @@ -152,6 +174,7 @@ impl HierarchicalKernel for MockKernel { source: &Self::SourceSummary, out: &mut Self::Output, ) { + self.far_calls.fetch_add(1, Ordering::Relaxed); let r2 = dist2(target.centroid, source.centroid); out[0] = source.moment / (T::ONE + r2); } @@ -165,6 +188,366 @@ impl HierarchicalKernel for MockKernel { } } +struct NoFieldTraversalKernel; + +impl HierarchicalKernel for NoFieldTraversalKernel { + type Scalar = f64; + type SourceGeometry = MockPoint; + type TargetGeometry = MockPoint; + type SourceMoment = f64; + type SourceSummary = SourceSummary; + type TargetSummary = TargetSummary; + type Output = [f64; 3]; + + fn summarize_leaf_sources( + &self, + source_ids: &[u32], + sources: S, + moments: M, + out: &mut Self::SourceSummary, + ) -> HierarchicalError + where + S: SourceCollection, + M: SourceMomentCollection, + { + let _ = (source_ids, sources, moments, out); + panic!("Skip::Both must not summarize source leaves") + } + + fn combine_source_summaries( + &self, + _children: &[Self::SourceSummary], + _out: &mut Self::SourceSummary, + ) -> HierarchicalError { + panic!("Skip::Both must not combine source summaries") + } + + fn summarize_leaf_targets( + &self, + _target_ids: &[u32], + _targets: &[Self::TargetGeometry], + _out: &mut Self::TargetSummary, + ) -> HierarchicalError { + panic!("Skip::Both must not summarize field targets") + } + + fn eval_near( + &self, + _target: &Self::TargetGeometry, + _source: &Self::SourceGeometry, + _moment: &Self::SourceMoment, + _out: &mut Self::Output, + ) { + panic!("Skip::Both must not evaluate near interactions") + } + + fn eval_far( + &self, + _target: &Self::TargetSummary, + _source: &Self::SourceSummary, + _out: &mut Self::Output, + ) { + panic!("Skip::Both must not evaluate far interactions") + } + + fn accept_far( + &self, + _target_aabb: Aabb, + _source_aabb: Aabb, + _source: &Self::SourceSummary, + _theta: f64, + ) -> bool { + panic!("Skip::Both must not traverse the source tree") + } + + fn zero_output(&self, _out: &mut Self::Output) { + panic!("Skip::Both must zero component arrays without invoking the kernel") + } + + fn accumulate(&self, _out: &mut Self::Output, _contribution: &Self::Output) { + panic!("Skip::Both must not accumulate field interactions") + } +} + +#[test] +fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output() { + let kernel = MockKernel::::new(); + let sources = points_f64(&[[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]); + let targets = points_f64(&[[0.0, 0.0, 0.0]]); + let moments = [2.0, 3.0]; + let source_tree = ClusterTree::build(sources.as_slice()).unwrap(); + let mut summaries = SourceNodeSummaries::>::new(source_tree.as_view()); + assert_eq!( + update_summaries( + &kernel, + source_tree.as_view(), + sources.as_slice(), + &moments, + &mut summaries.node_summaries, + ), + HierarchicalError::Ok + ); + + let mut contribution = [[0.0; 1]]; + let mut scratch = EvaluationScratch { + contribution: &mut contribution, + }; + let mut full = [0.0]; + assert_eq!( + super::eval( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + None, + [&mut full], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 1); + assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 1); + + kernel.near_calls.store(0, Ordering::Relaxed); + kernel.far_calls.store(0, Ordering::Relaxed); + let mut far_only = [0.0]; + assert_eq!( + super::eval( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Some(Skip::Near), + [&mut far_only], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 1); + + kernel.near_calls.store(0, Ordering::Relaxed); + kernel.far_calls.store(0, Ordering::Relaxed); + let mut near_only = [0.0]; + assert_eq!( + super::eval( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Some(Skip::Far), + [&mut near_only], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 1); + assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 0); + assert!((full[0] - near_only[0] - far_only[0]).abs() < 1.0e-15); + + kernel.near_calls.store(0, Ordering::Relaxed); + kernel.far_calls.store(0, Ordering::Relaxed); + let mut far_only_par = [0.0]; + assert_eq!( + super::eval_par( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Some(Skip::Near), + [&mut far_only_par], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 1); + assert_eq!(far_only_par, far_only); +} + +#[test] +fn skip_both_evaluators_bypass_summaries_scratch_and_traversal() { + let sources = points_f64(&[[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]); + let targets = points_f64(&[[0.0, 0.0, 0.0], [5.0, 0.0, 0.0]]); + let moments = [2.0, 3.0]; + let source_tree = ClusterTree::build(sources.as_slice()).unwrap(); + let mut contributions: [[f64; 3]; 0] = []; + let mut scratch = EvaluationScratch { + contribution: &mut contributions, + }; + let mut out0 = [f64::NAN; 2]; + let mut out1 = [f64::NAN; 2]; + let mut out2 = [f64::NAN; 2]; + + assert_eq!( + super::eval( + &NoFieldTraversalKernel, + source_tree.as_view(), + &[], + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Some(Skip::Both), + [&mut out0, &mut out1, &mut out2], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(out0, [0.0; 2]); + assert_eq!(out1, [0.0; 2]); + assert_eq!(out2, [0.0; 2]); + + out0.fill(f64::NAN); + out1.fill(f64::NAN); + out2.fill(f64::NAN); + assert_eq!( + super::eval_par( + &NoFieldTraversalKernel, + source_tree.as_view(), + &[], + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Some(Skip::Both), + [&mut out0, &mut out1, &mut out2], + &mut scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(out0, [0.0; 2]); + assert_eq!(out1, [0.0; 2]); + assert_eq!(out2, [0.0; 2]); +} + +#[test] +fn one_shot_skip_both_bypasses_field_traversal_and_zeroes_outputs() { + let sources = points_f64(&[[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]); + let targets = points_f64(&[[0.0, 0.0, 0.0], [5.0, 0.0, 0.0]]); + let moments = [2.0, 3.0]; + let mut out0 = [f64::NAN; 2]; + let mut out1 = [f64::NAN; 2]; + let mut out2 = [f64::NAN; 2]; + + let diagnostics = super::convenience::one_shot_vec3( + NoFieldTraversalKernel, + sources.as_slice(), + moments.as_slice(), + targets.as_slice(), + BuildMethod::LongestAxis, + 0.5, + true, + Some(Skip::Both), + false, + (&mut out0, &mut out1, &mut out2), + ) + .unwrap(); + + assert_eq!(out0, [0.0; 2]); + assert_eq!(out1, [0.0; 2]); + assert_eq!(out2, [0.0; 2]); + assert_eq!(diagnostics.source_count, 2); + assert_eq!(diagnostics.target_count, 2); + assert!(diagnostics.traversal_diagnostics().is_none()); + assert!(diagnostics.accepted_levels().is_none()); + assert!(diagnostics.near_field_interaction_map().is_none()); +} + +#[test] +fn one_shot_skip_both_collects_requested_traversal_diagnostics() { + let source_x = [0.0, 10.0]; + let source_yz = [0.0; 2]; + let moment_xz = [0.0; 2]; + let moment_y = [2.0, 3.0]; + let outer_radius = [0.0; 2]; + let target_x = [0.0, 10.0, 5.0]; + let target_yz = [0.0; 3]; + + for par in [false, true] { + let mut out0 = [f64::NAN; 3]; + let mut out1 = [f64::NAN; 3]; + let mut out2 = [f64::NAN; 3]; + let diagnostics = super::vector_potential_dipole_hierarchical_with_skip( + (&source_x, &source_yz, &source_yz), + (&moment_xz, &moment_y, &moment_xz), + (&target_x, &target_yz, &target_yz), + &outer_radius, + BuildMethod::LongestAxis, + 0.5, + par, + Skip::Both, + true, + (&mut out0, &mut out1, &mut out2), + ) + .unwrap(); + + assert_eq!(out0, [0.0; 3]); + assert_eq!(out1, [0.0; 3]); + assert_eq!(out2, [0.0; 3]); + assert_eq!(diagnostics.accepted_levels().unwrap().len(), 3); + let interaction_map = diagnostics.near_field_interaction_map().unwrap(); + assert_eq!(interaction_map.row_indices, vec![0, 1]); + assert_eq!(interaction_map.column_pointers, vec![0, 1, 2, 2]); + } +} + +#[test] +fn traversal_diagnostics_returns_canonical_near_field_csc_pattern() { + let kernel = MockKernel::::new(); + let sources = points_f64(&[[0.0, 0.0, 0.0], [10.0, 0.0, 0.0]]); + let targets = points_f64(&[[0.0, 0.0, 0.0], [10.0, 0.0, 0.0], [5.0, 0.0, 0.0]]); + let moments = [2.0, 3.0]; + let source_tree = ClusterTree::build(sources.as_slice()).unwrap(); + let mut summaries = SourceNodeSummaries::>::new(source_tree.as_view()); + assert_eq!( + update_summaries( + &kernel, + source_tree.as_view(), + sources.as_slice(), + &moments, + &mut summaries.node_summaries, + ), + HierarchicalError::Ok + ); + + let diagnostics = super::traversal_diagnostics( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + targets.as_slice(), + 0.5, + ) + .unwrap(); + let diagnostics_par = super::traversal_diagnostics_par( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + targets.as_slice(), + 0.5, + ) + .unwrap(); + assert_eq!(diagnostics, diagnostics_par); + + let map = diagnostics.near_field_interaction_map; + assert_eq!(map.source_count, 2); + assert_eq!(map.target_count, 3); + assert_eq!(map.row_indices, vec![0, 1]); + assert_eq!(map.column_pointers, vec![0, 1, 2, 2]); +} + #[test] fn hierarchical_error_raw_codes_keep_kernel_slots_first() { assert_eq!(HierarchicalError::KernelError0 as u32, 0); @@ -263,6 +646,7 @@ where targets, moments, theta, + None, column_slices, scratch, ); @@ -310,6 +694,7 @@ where targets, moments, theta, + None, column_slices, scratch, ); diff --git a/src/physics/linear_filament.rs b/src/physics/linear_filament.rs index 26b7f4c..e317f71 100644 --- a/src/physics/linear_filament.rs +++ b/src/physics/linear_filament.rs @@ -17,6 +17,38 @@ use crate::{MU0_OVER_4PI, macros::*}; /// (m) minimum representable nonzero wire thickness. const MIN_WIRE_THICKNESS: f64 = 1e-10; +/// Three-point target-segment integral for one unit-current source segment. +#[inline] +fn inductance_linear_filament_pair( + src_start: (f64, f64, f64), + src_delta: (f64, f64, f64), + wire_radius_src: f64, + tgt_start: (f64, f64, f64), + tgt_delta: (f64, f64, f64), +) -> f64 { + let src_end = ( + src_start.0 + src_delta.0, + src_start.1 + src_delta.1, + src_start.2 + src_delta.2, + ); + let gl3_unit = gauss_legendre_unit_interval_table(GaussLegendreRule::Gauss3); + let mut inductance = 0.0; + for &[tq, wq] in gl3_unit { + let obs = ( + tgt_delta.0.mul_add(tq, tgt_start.0), + tgt_delta.1.mul_add(tq, tgt_start.1), + tgt_delta.2.mul_add(tq, tgt_start.2), + ); + let (ax, ay, az) = vector_potential_linear_filament_scalar( + (src_start, src_end, 1.0), + wire_radius_src, + obs, + ); + inductance += wq * (ax * tgt_delta.0 + ay * tgt_delta.1 + az * tgt_delta.2); + } + inductance +} + /// Estimate the inductive coupling between two piecewise-linear current filaments. /// /// This uses the vector-potential line-integral form @@ -69,35 +101,17 @@ pub fn inductance_piecewise_linear_filaments( let (xfil0, yfil0, zfil0) = xyzfil0; let (dlxfil0, dlyfil0, dlzfil0) = dlxyzfil0; let mut inductance = 0.0; // [H] - let gl3_unit = gauss_legendre_unit_interval_table(GaussLegendreRule::Gauss3); - for j in 0..m { - let dltgt = (dlxfil1[j], dlyfil1[j], dlzfil1[j]); // [m] - for &[tq, wq] in gl3_unit { - let obs = ( - dltgt.0.mul_add(tq, xfil1[j]), // [m] - dltgt.1.mul_add(tq, yfil1[j]), // [m] - dltgt.2.mul_add(tq, zfil1[j]), // [m] + let tgt_start = (xfil1[j], yfil1[j], zfil1[j]); + let tgt_delta = (dlxfil1[j], dlyfil1[j], dlzfil1[j]); + for i in 0..n { + inductance += inductance_linear_filament_pair( + (xfil0[i], yfil0[i], zfil0[i]), + (dlxfil0[i], dlyfil0[i], dlzfil0[i]), + wire_radius[i], + tgt_start, + tgt_delta, ); - let mut ax = 0.0; // [V-s/m] - let mut ay = 0.0; // [V-s/m] - let mut az = 0.0; // [V-s/m] - - for i in 0..n { - let fil0 = (xfil0[i], yfil0[i], zfil0[i]); // [m] - let fil1 = ( - fil0.0 + dlxfil0[i], - fil0.1 + dlyfil0[i], - fil0.2 + dlzfil0[i], - ); // [m] - let (axc, ayc, azc) = - vector_potential_linear_filament_scalar((fil0, fil1, 1.0), wire_radius[i], obs); - ax += axc; // [V-s/m] - ay += ayc; // [V-s/m] - az += azc; // [V-s/m] - } - - inductance += wq * (ax * dltgt.0 + ay * dltgt.1 + az * dltgt.2); // [H] } } @@ -319,6 +333,201 @@ pub fn inductance_linear_filaments_matrix_par( Ok(()) } +/// Validate filament geometry and a canonical CSC interaction pattern. +fn validate_sparse_inductance_inputs( + xyzfil_tgt: (&[f64], &[f64], &[f64]), + dlxyzfil_tgt: (&[f64], &[f64], &[f64]), + xyzfil_src: (&[f64], &[f64], &[f64]), + dlxyzfil_src: (&[f64], &[f64], &[f64]), + wire_radius_src: &[f64], + row_indices: &[usize], + column_pointers: &[usize], + out: &[f64], +) -> Result<(usize, usize), &'static str> { + let ntgt = xyzfil_tgt.0.len(); + check_length!( + ntgt, + xyzfil_tgt.0, + xyzfil_tgt.1, + xyzfil_tgt.2, + dlxyzfil_tgt.0, + dlxyzfil_tgt.1, + dlxyzfil_tgt.2 + ); + + let nsrc = xyzfil_src.0.len(); + check_length!( + nsrc, + xyzfil_src.0, + xyzfil_src.1, + xyzfil_src.2, + dlxyzfil_src.0, + dlxyzfil_src.1, + dlxyzfil_src.2, + wire_radius_src + ); + check_length!(row_indices.len(), out); + + if column_pointers.len() != ntgt + 1 { + return Err("CSC column pointer length must equal target count plus one"); + } + if column_pointers.first() != Some(&0) { + return Err("CSC column pointers must start at zero"); + } + if column_pointers.last() != Some(&row_indices.len()) { + return Err("CSC final column pointer must equal the stored-entry count"); + } + if column_pointers + .windows(2) + .any(|pointers| pointers[0] > pointers[1]) + { + return Err("CSC column pointers must be nondecreasing"); + } + + for column in 0..ntgt { + let rows = &row_indices[column_pointers[column]..column_pointers[column + 1]]; + if rows.iter().any(|&row| row >= nsrc) { + return Err("CSC row index exceeds the source count"); + } + if rows.windows(2).any(|pair| pair[0] >= pair[1]) { + return Err("CSC row indices must be sorted and unique within each column"); + } + } + + Ok((nsrc, ntgt)) +} + +/// Evaluate one stored CSC source-target interaction. +#[inline] +fn inductance_linear_filaments_sparse_entry( + xyzfil_tgt: (&[f64], &[f64], &[f64]), + dlxyzfil_tgt: (&[f64], &[f64], &[f64]), + xyzfil_src: (&[f64], &[f64], &[f64]), + dlxyzfil_src: (&[f64], &[f64], &[f64]), + wire_radius_src: &[f64], + target: usize, + source: usize, +) -> f64 { + inductance_linear_filament_pair( + ( + xyzfil_src.0[source], + xyzfil_src.1[source], + xyzfil_src.2[source], + ), + ( + dlxyzfil_src.0[source], + dlxyzfil_src.1[source], + dlxyzfil_src.2[source], + ), + wire_radius_src[source], + ( + xyzfil_tgt.0[target], + xyzfil_tgt.1[target], + xyzfil_tgt.2[target], + ), + ( + dlxyzfil_tgt.0[target], + dlxyzfil_tgt.1[target], + dlxyzfil_tgt.2[target], + ), + ) +} + +/// Evaluate selected source-target filament inductances into CSC value storage. +/// +/// `row_indices` and `column_pointers` describe a canonical CSC matrix with shape +/// `(nsrc, ntgt)`, where source segments are rows and target segments are columns. Each stored +/// coordinate is evaluated with the finite-radius source kernel and three-point Gauss--Legendre +/// integration over the complete target segment. `out` has one value per stored coordinate and +/// retains explicit numerical zeros. +pub fn inductance_linear_filaments_sparse_csc( + xyzfil_tgt: (&[f64], &[f64], &[f64]), + dlxyzfil_tgt: (&[f64], &[f64], &[f64]), + xyzfil_src: (&[f64], &[f64], &[f64]), + dlxyzfil_src: (&[f64], &[f64], &[f64]), + wire_radius_src: &[f64], + row_indices: &[usize], + column_pointers: &[usize], + out: &mut [f64], +) -> Result<(), &'static str> { + let (_, ntgt) = validate_sparse_inductance_inputs( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + row_indices, + column_pointers, + out, + )?; + + for target in 0..ntgt { + for entry in column_pointers[target]..column_pointers[target + 1] { + out[entry] = inductance_linear_filaments_sparse_entry( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + target, + row_indices[entry], + ); + } + } + + Ok(()) +} + +/// Parallel variant of [`inductance_linear_filaments_sparse_csc`]. +pub fn inductance_linear_filaments_sparse_csc_par( + xyzfil_tgt: (&[f64], &[f64], &[f64]), + dlxyzfil_tgt: (&[f64], &[f64], &[f64]), + xyzfil_src: (&[f64], &[f64], &[f64]), + dlxyzfil_src: (&[f64], &[f64], &[f64]), + wire_radius_src: &[f64], + row_indices: &[usize], + column_pointers: &[usize], + out: &mut [f64], +) -> Result<(), &'static str> { + validate_sparse_inductance_inputs( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + row_indices, + column_pointers, + out, + )?; + + // Partition by stored entries instead of columns so a single dense column still uses all + // workers. Each chunk locates its first target once, then follows the CSC boundaries linearly. + let entry_chunk_size = chunksize(out.len()); + out.par_chunks_mut(entry_chunk_size) + .enumerate() + .for_each(|(chunk_id, values)| { + let entry_start = chunk_id * entry_chunk_size; + let mut target = column_pointers.partition_point(|&pointer| pointer <= entry_start) - 1; + for (offset, value) in values.iter_mut().enumerate() { + let entry = entry_start + offset; + while column_pointers[target + 1] <= entry { + target += 1; + } + *value = inductance_linear_filaments_sparse_entry( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + target, + row_indices[entry], + ); + } + }); + + Ok(()) +} + /// Biot-Savart calculation for B-field contribution from many current filament /// segments to many observation points. /// @@ -2848,4 +3057,91 @@ mod test { assert!(approx(out[j], contracted, 1e-12, 1e-15)); } } + + #[test] + fn test_sparse_csc_inductance_matches_selected_dense_entries() { + let xyzsrc = (&[0.0, 1.0, 2.0][..], &[0.0, 0.1, -0.1][..], &[0.0; 3][..]); + let dlxyzsrc = (&[0.0; 3][..], &[0.0; 3][..], &[0.8; 3][..]); + let wire_radius = [0.01, 0.02, 0.03]; + let xyztgt = ( + &[0.2, 1.2, 2.2, 3.2][..], + &[0.3, -0.2, 0.1, 0.0][..], + &[0.1, 0.2, -0.1, 0.3][..], + ); + let dlxyztgt = (&[0.1; 4][..], &[0.05; 4][..], &[0.4, 0.4, 0.0, 0.4][..]); + let row_indices = [0, 2, 1, 0, 1, 2]; + let column_pointers = [0, 2, 2, 3, 6]; + + let mut dense = vec![0.0; 3 * 4]; + inductance_linear_filaments_matrix( + xyztgt, + dlxyztgt, + xyzsrc, + dlxyzsrc, + &wire_radius, + &mut dense, + ) + .unwrap(); + + let mut sparse = vec![f64::NAN; row_indices.len()]; + let mut sparse_par = vec![f64::NAN; row_indices.len()]; + inductance_linear_filaments_sparse_csc( + xyztgt, + dlxyztgt, + xyzsrc, + dlxyzsrc, + &wire_radius, + &row_indices, + &column_pointers, + &mut sparse, + ) + .unwrap(); + inductance_linear_filaments_sparse_csc_par( + xyztgt, + dlxyztgt, + xyzsrc, + dlxyzsrc, + &wire_radius, + &row_indices, + &column_pointers, + &mut sparse_par, + ) + .unwrap(); + + assert_eq!(sparse, sparse_par); + for target in 0..4 { + for entry in column_pointers[target]..column_pointers[target + 1] { + let source = row_indices[entry]; + assert!(approx( + sparse[entry], + dense[source * 4 + target], + 1e-14, + 1e-18 + )); + } + } + assert_eq!(sparse[2], 0.0); + } + + #[test] + fn test_sparse_csc_inductance_rejects_noncanonical_pattern() { + let xyz = (&[0.0, 1.0][..], &[0.0; 2][..], &[0.0; 2][..]); + let dlxyz = (&[0.0; 2][..], &[0.0; 2][..], &[1.0; 2][..]); + let mut out = [0.0; 2]; + let err = inductance_linear_filaments_sparse_csc( + xyz, + dlxyz, + xyz, + dlxyz, + &[0.01; 2], + &[1, 1], + &[0, 2, 2], + &mut out, + ) + .unwrap_err(); + assert_eq!( + err, + "CSC row indices must be sorted and unique within each column" + ); + } } diff --git a/src/python.rs b/src/python.rs index 282b5d2..70258eb 100644 --- a/src/python.rs +++ b/src/python.rs @@ -127,6 +127,7 @@ struct HierarchicalDiagnostics { target_count: usize, source_tree: Option>, accepted_levels: Option>>, + near_field_interaction_map: Option>, } #[pymethods] @@ -168,6 +169,14 @@ impl HierarchicalDiagnostics { .as_ref() .map(|value| value.clone_ref(py)) } + + #[getter] + /// Return the direct near-field interaction pattern as a SciPy CSC matrix. + fn near_field_interaction_map(&self, py: Python<'_>) -> Option> { + self.near_field_interaction_map + .as_ref() + .map(|value| value.clone_ref(py)) + } } #[pyclass(module = "cfsem")] @@ -565,6 +574,7 @@ fn solve_result_from_field( target_count: usize, source_tree: Option>, accepted_levels: Option>>, + near_field_interaction_map: Option>, ) -> PyResult> { let field = PyTuple::new(py, [field.0, field.1, field.2])? .unbind() @@ -578,6 +588,7 @@ fn solve_result_from_field( target_count, source_tree, accepted_levels, + near_field_interaction_map, }, )?; Py::new(py, SolveResult { field, diagnostics }) @@ -625,93 +636,72 @@ fn parse_build_method( } } -/// Compute accepted source-node levels for hierarchical diagnostic output. -fn accepted_levels_diagnostic( - kernel: K, - source_tree: &physics::hierarchical::tree::ClusterTree, - sources: S, - targets: C, - moments: M, - theta: f64, -) -> PyResult> -where - K: physics::hierarchical::kernel::HierarchicalKernel + Sync, - S: physics::hierarchical::kernel::SourceCollection + Copy, - M: physics::hierarchical::kernel::SourceMomentCollection + Copy, - K::TargetGeometry: Copy, - C: physics::hierarchical::kernel::TargetCollection, -{ - let mut source_summaries = - physics::hierarchical::evaluator::SourceNodeSummaries::::new(source_tree.as_view()); - let mut err = physics::hierarchical::evaluator::update_summaries( - &kernel, - source_tree.as_view(), - sources, - moments, - &mut source_summaries.node_summaries, - ); - if err != physics::hierarchical::kernel::HierarchicalError::Ok { - return Err(py_hierarchical_error("source summary update", err)); - } - - let mut out = vec![0.0; physics::hierarchical::kernel::TargetCollection::::len(targets)]; - err = physics::hierarchical::evaluator::accepted_levels( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - targets, - theta, - &mut out, - ); - if err != physics::hierarchical::kernel::HierarchicalError::Ok { - return Err(py_hierarchical_error("source-level diagnostic", err)); +/// Parse the optional hierarchical interaction class to omit. +fn parse_hierarchical_skip( + skip: Option<&str>, +) -> PyResult> { + match skip { + None => Ok(None), + Some("near") => Ok(Some(physics::hierarchical::kernel::Skip::Near)), + Some("far") => Ok(Some(physics::hierarchical::kernel::Skip::Far)), + Some("both") => Ok(Some(physics::hierarchical::kernel::Skip::Both)), + Some(value) => Err(PyInteropError::ValueError { + msg: format!( + "Unsupported hierarchical skip value: {value}. Expected 'near', 'far', 'both', or None." + ), + } + .into()), } - Ok(out) } -type OptionalDiagnosticsPy = (Option>, Option>>); - -struct HierarchicalDiagnosticRequest<'a, K, S, C, M> { - kernel: K, - source_tree: &'a physics::hierarchical::tree::ClusterTree, - sources: S, - targets: C, - moments: M, - theta: f64, -} +/// Move a CSC interaction pattern into a SciPy sparse matrix. +fn near_field_interaction_map_object( + py: Python<'_>, + map: physics::hierarchical::NearFieldInteractionMap, +) -> PyResult> { + let data: Py = PyArray1::from_vec(py, vec![true; map.row_indices.len()]) + .unbind() + .into(); + let row_indices: Py = PyArray1::from_vec(py, map.row_indices).unbind().into(); + let column_pointers: Py = PyArray1::from_vec(py, map.column_pointers).unbind().into(); + let csc_arrays = PyTuple::new(py, [data, row_indices, column_pointers])?; + let constructor = py.import("scipy.sparse")?.getattr("csc_matrix")?; + Ok(constructor + .call1((csc_arrays, (map.source_count, map.target_count)))? + .unbind()) +} + +type OptionalDiagnosticsPy = ( + Option>, + Option>>, + Option>, +); -/// Return diagnostics only when requested by the Python caller. -fn optional_hierarchical_diagnostics( +/// Convert optional diagnostics collected by a Rust one-shot solve for Python. +fn hierarchical_diagnostics_objects( py: Python<'_>, - extra_diagnostics: bool, - request: HierarchicalDiagnosticRequest<'_, K, S, C, M>, + diagnostics: &mut physics::hierarchical::Diagnostics, ) -> PyResult where - K: physics::hierarchical::kernel::HierarchicalKernel + Sync, - S: physics::hierarchical::kernel::SourceCollection + Copy, - M: physics::hierarchical::kernel::SourceMomentCollection + Copy, - K::TargetGeometry: Copy, - C: physics::hierarchical::kernel::TargetCollection, + K: physics::hierarchical::kernel::HierarchicalKernel, { - if !extra_diagnostics { - return Ok((None, None)); - } - - let levels = accepted_levels_diagnostic( - request.kernel, - request.source_tree, - request.sources, - request.targets, - request.moments, - request.theta, - )?; + let Some(traversal) = diagnostics.take_traversal_diagnostics() else { + return Ok((None, None, None)); + }; Ok(( - Some(source_tree_diagnostics_object(py, request.source_tree)?), - Some(PyArray1::from_vec(py, levels).unbind()), + Some(source_tree_diagnostics_object( + py, + &diagnostics.source_tree, + )?), + Some(PyArray1::from_vec(py, traversal.accepted_levels).unbind()), + Some(near_field_interaction_map_object( + py, + traversal.near_field_interaction_map, + )?), )) } -#[pyfunction(signature = (loc, moment, obs, outer_radius, theta=0.01, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (loc, moment, obs, outer_radius, theta=0.01, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate dipole flux density with the hierarchical solver from Python inputs. fn flux_density_dipole_hierarchical( py: Python<'_>, @@ -740,54 +730,45 @@ fn flux_density_dipole_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let loc = read_xyz_tuple(py, &loc, "loc")?; let moment = read_xyz_tuple(py, &moment, "moment")?; let outer_radius = read_float_input_array1(py, &outer_radius, "outer_radius")?; let obs = read_xyz_tuple(py, &obs, "obs")?; let construction_method = parse_build_method(construction_method)?; - let (field, diagnostics) = + let skip = parse_hierarchical_skip(skip)?; + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.len(), "flux_density", |out| { - physics::hierarchical::flux_density_dipole_hierarchical( - loc.as_tuple(), - moment.as_tuple(), - obs.as_tuple(), - outer_radius.as_slice(), - construction_method, - theta, - par, - out, - ) + match skip { + Some(skip) => physics::hierarchical::flux_density_dipole_hierarchical_with_skip( + loc.as_tuple(), + moment.as_tuple(), + obs.as_tuple(), + outer_radius.as_slice(), + construction_method, + theta, + par, + skip, + extra_diagnostics, + out, + ), + None => physics::hierarchical::flux_density_dipole_hierarchical( + loc.as_tuple(), + moment.as_tuple(), + obs.as_tuple(), + outer_radius.as_slice(), + construction_method, + theta, + par, + extra_diagnostics, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical dipole flux density", err)) })?; - let sources = physics::hierarchical::kernels::DipoleSources::new( - loc.as_tuple().0, - loc.as_tuple().1, - loc.as_tuple().2, - outer_radius.as_slice(), - ); - let targets = physics::hierarchical::kernels::DipoleTargets::new( - obs.as_tuple().0, - obs.as_tuple().1, - obs.as_tuple().2, - ); - let moments = physics::hierarchical::kernels::DipoleMoments::new( - moment.as_tuple().0, - moment.as_tuple().1, - moment.as_tuple().2, - ); - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::DipoleFluxDensityKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -797,10 +778,11 @@ fn flux_density_dipole_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } -#[pyfunction(signature = (loc, moment, obs, outer_radius, theta=0.01, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (loc, moment, obs, outer_radius, theta=0.01, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate dipole vector potential with the hierarchical solver from Python inputs. fn vector_potential_dipole_hierarchical( py: Python<'_>, @@ -829,54 +811,47 @@ fn vector_potential_dipole_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let loc = read_xyz_tuple(py, &loc, "loc")?; let moment = read_xyz_tuple(py, &moment, "moment")?; let outer_radius = read_float_input_array1(py, &outer_radius, "outer_radius")?; let obs = read_xyz_tuple(py, &obs, "obs")?; let construction_method = parse_build_method(construction_method)?; - let (field, diagnostics) = + let skip = parse_hierarchical_skip(skip)?; + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.len(), "vector_potential", |out| { - physics::hierarchical::vector_potential_dipole_hierarchical( - loc.as_tuple(), - moment.as_tuple(), - obs.as_tuple(), - outer_radius.as_slice(), - construction_method, - theta, - par, - out, - ) + match skip { + Some(skip) => { + physics::hierarchical::vector_potential_dipole_hierarchical_with_skip( + loc.as_tuple(), + moment.as_tuple(), + obs.as_tuple(), + outer_radius.as_slice(), + construction_method, + theta, + par, + skip, + extra_diagnostics, + out, + ) + } + None => physics::hierarchical::vector_potential_dipole_hierarchical( + loc.as_tuple(), + moment.as_tuple(), + obs.as_tuple(), + outer_radius.as_slice(), + construction_method, + theta, + par, + extra_diagnostics, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical dipole vector potential", err)) })?; - let sources = physics::hierarchical::kernels::DipoleSources::new( - loc.as_tuple().0, - loc.as_tuple().1, - loc.as_tuple().2, - outer_radius.as_slice(), - ); - let targets = physics::hierarchical::kernels::DipoleTargets::new( - obs.as_tuple().0, - obs.as_tuple().1, - obs.as_tuple().2, - ); - let moments = physics::hierarchical::kernels::DipoleMoments::new( - moment.as_tuple().0, - moment.as_tuple().1, - moment.as_tuple().2, - ); - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::DipoleVectorPotentialKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -886,10 +861,11 @@ fn vector_potential_dipole_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } -#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate linear-filament flux density with the hierarchical solver from Python inputs. fn flux_density_linear_filament_hierarchical( py: Python<'_>, @@ -919,6 +895,7 @@ fn flux_density_linear_filament_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let xyzp = read_xyz_tuple(py, &xyzp, "xyzp")?; let xyzfil = read_xyz_tuple(py, &xyzfil, "xyzfil")?; @@ -926,43 +903,42 @@ fn flux_density_linear_filament_hierarchical( let ifil = read_float_input_array1(py, &ifil, "ifil")?; let wire_radius = read_float_input_array1(py, &wire_radius, "wire_radius")?; let construction_method = parse_build_method(construction_method)?; - let (field, diagnostics) = + let skip = parse_hierarchical_skip(skip)?; + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, xyzp.len(), "flux_density", |out| { - physics::hierarchical::flux_density_linear_filament_hierarchical( - xyzp.as_tuple(), - xyzfil.as_tuple(), - dlxyzfil.as_tuple(), - ifil.as_slice(), - wire_radius.as_slice(), - construction_method, - theta, - par, - out, - ) + match skip { + Some(skip) => { + physics::hierarchical::flux_density_linear_filament_hierarchical_with_skip( + xyzp.as_tuple(), + xyzfil.as_tuple(), + dlxyzfil.as_tuple(), + ifil.as_slice(), + wire_radius.as_slice(), + construction_method, + theta, + par, + skip, + extra_diagnostics, + out, + ) + } + None => physics::hierarchical::flux_density_linear_filament_hierarchical( + xyzp.as_tuple(), + xyzfil.as_tuple(), + dlxyzfil.as_tuple(), + ifil.as_slice(), + wire_radius.as_slice(), + construction_method, + theta, + par, + extra_diagnostics, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical linear-filament flux density", err)) })?; - let sources = physics::hierarchical::kernels::LinearFilamentSources::new( - xyzfil.as_tuple(), - dlxyzfil.as_tuple(), - wire_radius.as_slice(), - ); - let targets = physics::hierarchical::kernels::DipoleTargets::new( - xyzp.as_tuple().0, - xyzp.as_tuple().1, - xyzp.as_tuple().2, - ); - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::LinearFilamentFluxDensityKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments: ifil.as_slice(), - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -972,10 +948,11 @@ fn flux_density_linear_filament_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } -#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate linear-filament vector potential with the hierarchical solver from Python inputs. fn vector_potential_linear_filament_hierarchical( py: Python<'_>, @@ -1005,6 +982,7 @@ fn vector_potential_linear_filament_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let xyzp = read_xyz_tuple(py, &xyzp, "xyzp")?; let xyzfil = read_xyz_tuple(py, &xyzfil, "xyzfil")?; @@ -1012,46 +990,44 @@ fn vector_potential_linear_filament_hierarchical( let ifil = read_float_input_array1(py, &ifil, "ifil")?; let wire_radius = read_float_input_array1(py, &wire_radius, "wire_radius")?; let construction_method = parse_build_method(construction_method)?; - let (field, diagnostics) = + let skip = parse_hierarchical_skip(skip)?; + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, xyzp.len(), "vector_potential", |out| { - physics::hierarchical::vector_potential_linear_filament_hierarchical( - xyzp.as_tuple(), - xyzfil.as_tuple(), - dlxyzfil.as_tuple(), - ifil.as_slice(), - wire_radius.as_slice(), - construction_method, - theta, - par, - out, - ) + match skip { + Some(skip) => { + physics::hierarchical::vector_potential_linear_filament_hierarchical_with_skip( + xyzp.as_tuple(), + xyzfil.as_tuple(), + dlxyzfil.as_tuple(), + ifil.as_slice(), + wire_radius.as_slice(), + construction_method, + theta, + par, + skip, + extra_diagnostics, + out, + ) + } + None => physics::hierarchical::vector_potential_linear_filament_hierarchical( + xyzp.as_tuple(), + xyzfil.as_tuple(), + dlxyzfil.as_tuple(), + ifil.as_slice(), + wire_radius.as_slice(), + construction_method, + theta, + par, + extra_diagnostics, + out, + ), + } .map_err(|err| { py_hierarchical_error("hierarchical linear-filament vector potential", err) }) })?; - let sources = physics::hierarchical::kernels::LinearFilamentSources::new( - xyzfil.as_tuple(), - dlxyzfil.as_tuple(), - wire_radius.as_slice(), - ); - let targets = physics::hierarchical::kernels::DipoleTargets::new( - xyzp.as_tuple().0, - xyzp.as_tuple().1, - xyzp.as_tuple().2, - ); - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::LinearFilamentVectorPotentialKernel::::new( - ), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments: ifil.as_slice(), - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1061,10 +1037,11 @@ fn vector_potential_linear_filament_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } -#[pyfunction(signature = (obs, nodes, triangles, s, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (obs, nodes, triangles, s, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate triangle-mesh flux density with the hierarchical solver from Python inputs. fn flux_density_triangle_mesh_hierarchical( py: Python<'_>, @@ -1081,6 +1058,7 @@ fn flux_density_triangle_mesh_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let obs = read_matrix3_input(py, &obs, "obs")?; let nodes = read_matrix3_input(py, &nodes, "nodes")?; @@ -1088,11 +1066,12 @@ fn flux_density_triangle_mesh_hierarchical( let mesh = borrowed_triangle_mesh_view(&nodes, &triangles)?; let s = read_float_input_array1(py, &s, "s")?; let construction_method = parse_build_method(construction_method)?; + let skip = parse_hierarchical_skip(skip)?; let sources = physics::hierarchical::kernels::BoundaryElementTriangles::new(&mesh); let targets = physics::hierarchical::kernels::DipoleTargetRows::new(obs.as_slice()); let moments = physics::hierarchical::kernels::BoundaryElementNodalValues::new(sources, s.as_slice()); - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.nrows(), "flux_density", |out| { physics::hierarchical::convenience::one_shot_vec3( physics::hierarchical::kernels::BoundaryElementFluxDensityKernel::::new(), @@ -1102,22 +1081,14 @@ fn flux_density_triangle_mesh_hierarchical( construction_method, theta, par, + skip, + extra_diagnostics, out, ) .map_err(|err| py_hierarchical_error("hierarchical triangle-mesh flux density", err)) })?; - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::BoundaryElementFluxDensityKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1127,10 +1098,11 @@ fn flux_density_triangle_mesh_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } -#[pyfunction(signature = (obs, nodes, triangles, s, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false))] +#[pyfunction(signature = (obs, nodes, triangles, s, theta=0.05, construction_method="longest_axis", par=true, out=None, extra_diagnostics=false, skip=None))] /// Evaluate triangle-mesh vector potential with the hierarchical solver from Python inputs. fn vector_potential_triangle_mesh_hierarchical( py: Python<'_>, @@ -1147,6 +1119,7 @@ fn vector_potential_triangle_mesh_hierarchical( PyReadwriteArray1, )>, extra_diagnostics: bool, + skip: Option<&str>, ) -> PyResult> { let obs = read_matrix3_input(py, &obs, "obs")?; let nodes = read_matrix3_input(py, &nodes, "nodes")?; @@ -1154,11 +1127,12 @@ fn vector_potential_triangle_mesh_hierarchical( let mesh = borrowed_triangle_mesh_view(&nodes, &triangles)?; let s = read_float_input_array1(py, &s, "s")?; let construction_method = parse_build_method(construction_method)?; + let skip = parse_hierarchical_skip(skip)?; let sources = physics::hierarchical::kernels::BoundaryElementTriangles::new(&mesh); let targets = physics::hierarchical::kernels::DipoleTargetRows::new(obs.as_slice()); let moments = physics::hierarchical::kernels::BoundaryElementNodalValues::new(sources, s.as_slice()); - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.nrows(), "vector_potential", |out| { physics::hierarchical::convenience::one_shot_vec3( physics::hierarchical::kernels::BoundaryElementVectorPotentialKernel::::new(), @@ -1168,25 +1142,16 @@ fn vector_potential_triangle_mesh_hierarchical( construction_method, theta, par, + skip, + extra_diagnostics, out, ) .map_err(|err| { py_hierarchical_error("hierarchical triangle-mesh vector potential", err) }) })?; - let (source_tree, accepted_levels) = optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: - physics::hierarchical::kernels::BoundaryElementVectorPotentialKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - }, - )?; + let (source_tree, accepted_levels, near_field_interaction_map) = + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1196,6 +1161,7 @@ fn vector_potential_triangle_mesh_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -3275,6 +3241,64 @@ fn inductance_linear_filaments_matrix( Ok(PyArray1::from_vec(py, out).unbind()) } +#[pyfunction(signature = (xyzfil_tgt, dlxyzfil_tgt, xyzfil_src, dlxyzfil_src, wire_radius_src, row_indices, column_pointers, par=true))] +fn inductance_linear_filaments_sparse_csc( + py: Python<'_>, + xyzfil_tgt: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), + dlxyzfil_tgt: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), + xyzfil_src: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), + dlxyzfil_src: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), + wire_radius_src: PyReadonlyArray1, + row_indices: PyReadonlyArray1, + column_pointers: PyReadonlyArray1, + par: bool, +) -> PyResult>> { + _3tup_slice_ro!(xyzfil_tgt); + _3tup_slice_ro!(dlxyzfil_tgt); + _3tup_slice_ro!(xyzfil_src); + _3tup_slice_ro!(dlxyzfil_src); + let wire_radius_src = wire_radius_src.as_slice()?; + let row_indices = read_usize_indices("row_indices", row_indices)?; + let column_pointers = read_usize_indices("column_pointers", column_pointers)?; + let mut out = vec![0.0; row_indices.len()]; + + let func = match par { + true => physics::linear_filament::inductance_linear_filaments_sparse_csc_par, + false => physics::linear_filament::inductance_linear_filaments_sparse_csc, + }; + func( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + &row_indices, + &column_pointers, + &mut out, + ) + .map_err(|msg| PyInteropError::DimensionalityError { + msg: msg.to_string(), + })?; + + Ok(PyArray1::from_vec(py, out).unbind()) +} + /// Python bindings for cfsemrs::physics::gradshafranov::gs_operator_order2 #[pyfunction] fn gs_operator_order2( @@ -4596,6 +4620,10 @@ fn _cfsem<'py>(_py: Python, m: Bound<'py, PyModule>) -> PyResult<()> { inductance_linear_filaments_matrix, m.clone() )?)?; + m.add_function(wrap_pyfunction!( + inductance_linear_filaments_sparse_csc, + m.clone() + )?)?; m.add_function(wrap_pyfunction!( body_force_density_linear_filament, m.clone() diff --git a/test/test_electromagnetics.py b/test/test_electromagnetics.py index 9b27563..3db2506 100644 --- a/test/test_electromagnetics.py +++ b/test/test_electromagnetics.py @@ -1,6 +1,7 @@ """Tests of standalone electromagnetics calcs""" import numpy as np +import scipy.sparse as sparse from pytest import approx, mark, raises import cfsem @@ -840,9 +841,9 @@ def test_vector_potential_linear_self_inductance_against_wien(ndiscr, par): par=par, ) l_wien = float(cfsem.self_inductance_circular_ring_wien(major_radius, minor_radius)) # [H] - assert l_from_a == approx(l_wien, rel=8e-2), ( - f"ndiscr={ndiscr}, L_from_A={l_from_a:.6e}, L_wien={l_wien:.6e}" - ) + assert l_from_a == approx( + l_wien, rel=8e-2 + ), f"ndiscr={ndiscr}, L_from_A={l_from_a:.6e}, L_wien={l_wien:.6e}" @mark.parametrize("ndiscr_coarse", [100, 200, 400]) @@ -880,9 +881,7 @@ def test_linear_filament_self_inductance_against_wien(r, a, n): y = major_radius * np.sin(phi) z = np.zeros_like(x) - l_self = float( - cfsem.self_inductance_piecewise_linear_filaments((x, y, z), wire_radius=minor_radius) - ) + l_self = float(cfsem.self_inductance_piecewise_linear_filaments((x, y, z), wire_radius=minor_radius)) l_wien = float(cfsem.self_inductance_circular_ring_wien(major_radius, minor_radius)) assert l_self == approx(l_wien, rel=8e-2) @@ -1055,6 +1054,98 @@ def test_inductance_linear_filaments_matrix_contracts_to_vector(par): ) +@mark.parametrize("par", [True, False]) +@mark.parametrize("csc_type", [sparse.csc_matrix, sparse.csc_array]) +def test_inductance_linear_filaments_sparse_preserves_csc_pattern(par, csc_type): + xyzfil_src = ( + np.array([0.0, 1.0, 2.0]), + np.array([0.0, 0.1, -0.1]), + np.zeros(3), + ) + dlxyzfil_src = (np.zeros(3), np.zeros(3), np.full(3, 0.8)) + xyzfil_tgt = ( + np.array([0.2, 1.2, 2.2, 3.2]), + np.array([0.3, -0.2, 0.1, 0.0]), + np.array([0.1, 0.2, -0.1, 0.3]), + ) + dlxyzfil_tgt = ( + np.full(4, 0.1), + np.full(4, 0.05), + np.array([0.4, 0.4, 0.0, 0.4]), + ) + row_indices = np.array([0, 2, 1, 0, 1, 2], dtype=np.int32) + column_pointers = np.array([0, 2, 2, 3, 6], dtype=np.int32) + interaction_map = csc_type( + (np.full(row_indices.size, np.nan), row_indices, column_pointers), + shape=(3, 4), + ) + + actual = cfsem.inductance_linear_filaments_sparse( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + interaction_map, + wire_radius_src=np.array([0.01, 0.02, 0.03]), + par=par, + ) + dense = cfsem.inductance_linear_filaments( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src=np.array([0.01, 0.02, 0.03]), + par=par, + output="matrix", + ) + + assert sparse.isspmatrix_csc(actual) + assert actual.shape == interaction_map.shape + np.testing.assert_array_equal(actual.indices, interaction_map.indices) + np.testing.assert_array_equal(actual.indptr, interaction_map.indptr) + for target in range(actual.shape[1]): + start, end = actual.indptr[target : target + 2] + np.testing.assert_allclose( + actual.data[start:end], + dense[actual.indices[start:end], target], + rtol=1e-14, + atol=1e-18, + ) + assert actual.nnz == interaction_map.nnz + assert actual.data[2] == 0.0 + + scalar_radius = cfsem.inductance_linear_filaments_sparse( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + interaction_map, + wire_radius_src=0.01, + par=par, + ) + assert scalar_radius.shape == interaction_map.shape + + +def test_inductance_linear_filaments_sparse_validates_map(): + xyz = (np.array([0.0, 1.0]), np.zeros(2), np.zeros(2)) + dlxyz = (np.zeros(2), np.zeros(2), np.ones(2)) + csr_map = sparse.eye(2, format="csr") + with raises(TypeError, match="must be a scipy.sparse.csc_matrix or csc_array"): + cfsem.inductance_linear_filaments_sparse(xyz, dlxyz, xyz, dlxyz, csr_map) + + wrong_shape = sparse.eye(3, format="csc") + with raises(ValueError, match="must have shape"): + cfsem.inductance_linear_filaments_sparse(xyz, dlxyz, xyz, dlxyz, wrong_shape) + + duplicate_rows = sparse.csc_matrix( + (np.ones(2), np.array([0, 0]), np.array([0, 2, 2])), + shape=(2, 2), + ) + assert not duplicate_rows.has_canonical_format + with raises(ValueError, match="sorted, unique row indices"): + cfsem.inductance_linear_filaments_sparse(xyz, dlxyz, xyz, dlxyz, duplicate_rows) + + @mark.parametrize("ndiscr", [128, 200]) @mark.parametrize("par", [True, False]) def test_flux_density_linear_matrix_contracts_to_vector(ndiscr, par): diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 01a8ccb..7fe905d 100644 --- a/test/test_hierarchical.py +++ b/test/test_hierarchical.py @@ -2,6 +2,7 @@ import numpy as np import pytest +import scipy.sparse as sp import cfsem @@ -30,6 +31,15 @@ def _assert_returns_output_views(returned, out): assert np.shares_memory(returned_component, out_component) +def _assert_vec_zero(result): + for component in result.field: + np.testing.assert_array_equal(component, np.zeros_like(component)) + + +def _add_vec3(lhs, rhs): + return tuple(left + right for left, right in zip(lhs, rhs, strict=True)) + + def _assert_diagnostics(result, nsource, ntarget): assert result.diagnostics.construction_time >= 0.0 assert result.diagnostics.evaluation_time >= 0.0 @@ -42,6 +52,10 @@ def _assert_diagnostics(result, nsource, ntarget): assert result.diagnostics.source_tree[8].shape == result.diagnostics.source_tree[0].shape assert result.diagnostics.accepted_levels is not None assert result.diagnostics.accepted_levels.shape == (ntarget,) + interaction_map = result.diagnostics.near_field_interaction_map + assert sp.isspmatrix_csc(interaction_map) + assert interaction_map.shape == (nsource, ntarget) + assert interaction_map.has_canonical_format def test_hierarchical_dipoles_match_direct(): @@ -74,6 +88,10 @@ def test_hierarchical_dipoles_match_direct(): _assert_vec_close(result_a, direct_a) _assert_diagnostics(result_b, nsource=3, ntarget=4) _assert_diagnostics(result_a, nsource=3, ntarget=4) + interaction_map = result_a.diagnostics.near_field_interaction_map + assert interaction_map.nnz == 12 + np.testing.assert_array_equal(interaction_map.indices, np.tile(np.arange(3), 4)) + np.testing.assert_array_equal(interaction_map.indptr, np.arange(0, 13, 3)) out = (np.empty_like(obs[0]), np.empty_like(obs[0]), np.empty_like(obs[0])) returned = cfsem.flux_density_dipole_hierarchical( @@ -81,6 +99,7 @@ def test_hierarchical_dipoles_match_direct(): ) _assert_returns_output_views(returned, out) _assert_vec_close(out, direct_b) + assert returned.diagnostics.near_field_interaction_map is None def test_hierarchical_linear_filaments_match_direct(): @@ -116,6 +135,258 @@ def test_hierarchical_linear_filaments_match_direct(): _assert_diagnostics(result_a, nsource=2, ntarget=3) +@pytest.mark.parametrize("par", [False, True]) +def test_hierarchical_skip_decomposes_near_and_far_fields(par): + loc = ( + np.array([0.0, 0.1, 0.2, 10.0, 10.1, 10.2]), + np.zeros(6), + np.zeros(6), + ) + moment = (np.zeros(6), np.ones(6), np.ones(6)) + obs = (np.array([0.4]), np.array([0.3]), np.array([0.2])) + # Finite source bounds force the nearby leaves down the direct path while + # the compact cluster near x=10 is still accepted as far field. + outer_radius = np.full(6, 0.05) + + full = cfsem.vector_potential_dipole_hierarchical( + loc, moment, obs, outer_radius, theta=0.2, par=par, extra_diagnostics=True + ) + far_only = cfsem.vector_potential_dipole_hierarchical( + loc, + moment, + obs, + outer_radius, + theta=0.2, + par=par, + skip="near", + extra_diagnostics=True, + ) + near_only = cfsem.vector_potential_dipole_hierarchical( + loc, + moment, + obs, + outer_radius, + theta=0.2, + par=par, + skip="far", + extra_diagnostics=True, + ) + diagnostics_only = cfsem.vector_potential_dipole_hierarchical( + loc, + moment, + obs, + outer_radius, + theta=0.2, + par=par, + skip="both", + extra_diagnostics=True, + ) + + _assert_vec_close(full, _add_vec3(far_only.field, near_only.field)) + assert any(np.any(component != 0.0) for component in far_only.field) + assert any(np.any(component != 0.0) for component in near_only.field) + _assert_vec_zero(diagnostics_only) + for filtered in (far_only, near_only, diagnostics_only): + np.testing.assert_array_equal( + filtered.diagnostics.near_field_interaction_map.indices, + full.diagnostics.near_field_interaction_map.indices, + ) + np.testing.assert_array_equal( + filtered.diagnostics.near_field_interaction_map.indptr, + full.diagnostics.near_field_interaction_map.indptr, + ) + + zero_only = cfsem.vector_potential_dipole_hierarchical( + loc, moment, obs, outer_radius, theta=0.2, par=par, skip="both" + ) + _assert_vec_zero(zero_only) + assert zero_only.diagnostics.near_field_interaction_map is None + + +def test_hierarchical_skip_rejects_unknown_value(): + loc = (np.array([0.0]), np.array([0.0]), np.array([0.0])) + moment = (np.array([0.0]), np.array([0.0]), np.array([1.0])) + obs = (np.array([1.0]), np.array([0.0]), np.array([0.0])) + + with pytest.raises(ValueError, match="Unsupported hierarchical skip value"): + cfsem.vector_potential_dipole_hierarchical(loc, moment, obs, np.zeros(1), skip="not-an-interaction") + + +def test_near_field_interaction_map_uses_original_source_rows(): + loc = ( + np.array([10.0, 0.0]), + np.zeros(2), + np.zeros(2), + ) + moment = (np.zeros(2), np.ones(2), np.ones(2)) + obs = ( + np.array([0.0, 10.0, 5.0]), + np.zeros(3), + np.zeros(3), + ) + result = cfsem.vector_potential_dipole_hierarchical( + loc, + moment, + obs, + np.zeros(2), + theta=0.5, + par=False, + extra_diagnostics=True, + ) + + interaction_map = result.diagnostics.near_field_interaction_map + np.testing.assert_array_equal(interaction_map.indices, np.array([1, 0])) + np.testing.assert_array_equal(interaction_map.indptr, np.array([0, 1, 2, 2])) + np.testing.assert_array_equal(interaction_map.data, np.ones(2, dtype=bool)) + + +def test_all_hierarchical_methods_accept_skip_both(): + loc = (np.array([0.0]), np.array([0.0]), np.array([0.0])) + moment = (np.array([0.0]), np.array([0.0]), np.array([1.0])) + obs = (np.array([1.0]), np.array([0.2]), np.array([0.3])) + radius = np.zeros(1) + xyzfil = loc + dlxyzfil = (np.array([0.0]), np.array([0.0]), np.array([0.5])) + current = np.ones(1) + nodes = np.array([[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]) + triangles = np.array([[0, 1, 2]], dtype=np.int64) + stream_function = np.array([0.0, 1.0, 0.25]) + obs_rows = np.array([[0.2, 0.2, 0.5]]) + + results = ( + cfsem.flux_density_dipole_hierarchical(loc, moment, obs, radius, theta=0.0, skip="both"), + cfsem.vector_potential_dipole_hierarchical(loc, moment, obs, radius, theta=0.0, skip="both"), + cfsem.flux_density_linear_filament_hierarchical( + obs, xyzfil, dlxyzfil, current, radius, theta=0.0, skip="both" + ), + cfsem.vector_potential_linear_filament_hierarchical( + obs, xyzfil, dlxyzfil, current, radius, theta=0.0, skip="both" + ), + cfsem.flux_density_triangle_mesh_hierarchical( + obs_rows, nodes, triangles, stream_function, theta=0.0, skip="both" + ), + cfsem.vector_potential_triangle_mesh_hierarchical( + obs_rows, nodes, triangles, stream_function, theta=0.0, skip="both" + ), + ) + for result in results: + _assert_vec_zero(result) + + +def test_skip_both_zeroes_and_returns_supplied_output_arrays(): + loc = (np.array([0.0]), np.array([0.0]), np.array([0.0])) + moment = (np.array([0.0]), np.array([0.0]), np.array([1.0])) + obs = (np.array([1.0, 2.0]), np.zeros(2), np.zeros(2)) + out = tuple(np.full(2, np.nan) for _ in range(3)) + + result = cfsem.vector_potential_dipole_hierarchical( + loc, moment, obs, np.zeros(1), skip="both", out=out + ) + + _assert_returns_output_views(result, out) + _assert_vec_zero(result) + + +@pytest.mark.parametrize("par", [False, True]) +def test_near_field_map_drives_sparse_filament_inductance(par): + nsegment = 16 + phi = np.linspace(0.0, 2.0 * np.pi, nsegment + 1) + points = np.column_stack((np.cos(phi), np.sin(phi), np.zeros_like(phi))) + starts = points[:-1] + deltas = np.diff(points, axis=0) + midpoints = starts + 0.5 * deltas + xyzfil = _tuple_columns(starts) + dlxyzfil = _tuple_columns(deltas) + targets = _tuple_columns(midpoints) + current = np.linspace(0.8, 1.2, nsegment) + wire_radius = np.full(nsegment, 0.02) + + full = cfsem.vector_potential_linear_filament_hierarchical( + targets, + xyzfil, + dlxyzfil, + current, + wire_radius, + theta=0.35, + par=par, + ) + diagnostics_only = cfsem.vector_potential_linear_filament_hierarchical( + targets, + xyzfil, + dlxyzfil, + current, + wire_radius, + theta=0.35, + par=par, + extra_diagnostics=True, + skip="both", + ) + far_only = cfsem.vector_potential_linear_filament_hierarchical( + targets, + xyzfil, + dlxyzfil, + current, + wire_radius, + theta=0.35, + par=par, + skip="near", + ) + near_only = cfsem.vector_potential_linear_filament_hierarchical( + targets, + xyzfil, + dlxyzfil, + current, + wire_radius, + theta=0.35, + par=par, + skip="far", + ) + _assert_vec_close(full, _add_vec3(far_only.field, near_only.field)) + _assert_vec_zero(diagnostics_only) + + interaction_map = diagnostics_only.diagnostics.near_field_interaction_map + assert 0 < interaction_map.nnz < nsegment * nsegment + sparse_inductance = cfsem.inductance_linear_filaments_sparse( + xyzfil, + dlxyzfil, + xyzfil, + dlxyzfil, + interaction_map, + wire_radius_src=wire_radius, + par=par, + ) + dense_inductance = cfsem.inductance_linear_filaments( + xyzfil, + dlxyzfil, + xyzfil, + dlxyzfil, + wire_radius_src=wire_radius, + par=par, + output="matrix", + ) + for target in range(nsegment): + start, end = sparse_inductance.indptr[target : target + 2] + rows = sparse_inductance.indices[start:end] + np.testing.assert_allclose( + sparse_inductance.data[start:end], + dense_inductance[rows, target], + rtol=1e-14, + atol=1e-18, + ) + + expected_contraction = np.zeros(nsegment) + for target in range(nsegment): + start, end = interaction_map.indptr[target : target + 2] + rows = interaction_map.indices[start:end] + expected_contraction[target] = dense_inductance[rows, target] @ current[rows] + np.testing.assert_allclose( + sparse_inductance.T @ current, + expected_contraction, + rtol=1e-14, + atol=1e-18, + ) + + def test_hierarchical_construction_method_is_exposed(): xyzfil = ( np.array([0.0, 0.5, -0.2]),