From 9987a9ec9613de1e2cb5edd524036cd21e100e15 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:14:42 -0400 Subject: [PATCH 01/20] [nearfield-separation] Phase 1: filtered evaluator Add required Near/Far skip modes to serial and parallel hierarchical evaluation while preserving existing entry points. Plan: /Users/jlogan/git/cfsem-py/nearfield_separation_plan.md Co-Authored-By: Codex --- src/physics/hierarchical/evaluator.rs | 166 +++++++++++++++++++++++++- src/physics/hierarchical/kernel.rs | 9 ++ src/physics/hierarchical/mod.rs | 5 +- src/physics/hierarchical/tests.rs | 114 ++++++++++++++++++ 4 files changed, 289 insertions(+), 5 deletions(-) diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index 88b0557..32ef7ea 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -1,5 +1,5 @@ use super::{ - BoundedGeometry, ClusterTreeView, HierarchicalError, HierarchicalKernel, Scalar, + BoundedGeometry, ClusterTreeView, HierarchicalError, HierarchicalKernel, Scalar, Skip, SourceCollection, SourceMomentCollection, TargetCollection, }; use std::sync::atomic::{AtomicU32, Ordering}; @@ -104,6 +104,81 @@ pub fn eval( out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError +where + K: HierarchicalKernel, + T: Scalar, + K::TargetGeometry: Copy, + S: SourceCollection, + M: SourceMomentCollection, + C: TargetCollection, +{ + eval_optional_skip( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + None, + out, + scratch, + ) +} + +/// Evaluate vector-valued targets while omitting one interaction class. +/// +/// The source tree is still traversed normally so acceptance decisions do not +/// change. [`Skip::Near`] returns accepted far-summary contributions only; +/// [`Skip::Far`] returns direct leaf contributions only. +#[inline] +pub fn eval_with_skip( + kernel: &K, + source_tree: ClusterTreeView<'_, T>, + source_summaries: &[K::SourceSummary], + sources: S, + targets: C, + moments: M, + theta: T, + skip: Skip, + out: [&mut [T]; D], + scratch: &mut EvaluationScratch<'_, [T; D]>, +) -> HierarchicalError +where + K: HierarchicalKernel, + T: Scalar, + K::TargetGeometry: Copy, + S: SourceCollection, + M: SourceMomentCollection, + C: TargetCollection, +{ + eval_optional_skip( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + Some(skip), + out, + scratch, + ) +} + +#[inline] +fn eval_optional_skip( + kernel: &K, + source_tree: ClusterTreeView<'_, T>, + source_summaries: &[K::SourceSummary], + sources: S, + targets: C, + moments: M, + theta: T, + skip: Option, + out: [&mut [T]; D], + scratch: &mut EvaluationScratch<'_, [T; D]>, +) -> HierarchicalError where K: HierarchicalKernel, T: Scalar, @@ -124,6 +199,7 @@ where targets, moments, theta, + skip, out, scratch, ) @@ -139,6 +215,7 @@ fn eval_validated( targets: C, moments: M, theta: T, + skip: Option, out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError @@ -183,6 +260,7 @@ where target, moments, theta, + skip, &mut target_out, &mut scratch.contribution[0], &mut target_summary, @@ -213,6 +291,7 @@ fn eval_scalar( target: K::TargetGeometry, moments: M, theta: K::Scalar, + skip: Option, out: &mut K::Output, contribution: &mut K::Output, target_summary: &mut K::TargetSummary, @@ -239,13 +318,18 @@ where 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); + if skip != Some(Skip::Far) { + 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 { + if skip == Some(Skip::Near) { + continue; + } let start = source_tree.leaf_start[source_node_index] as usize; let count = leaf_count as usize; let end = start + count; @@ -285,6 +369,77 @@ pub fn eval_par( out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError +where + K: HierarchicalKernel + Sync, + T: Scalar, + K::TargetGeometry: Copy, + S: SourceCollection, + M: SourceMomentCollection, + C: TargetCollection, +{ + eval_par_optional_skip( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + None, + out, + scratch, + ) +} + +/// Evaluate vector-valued targets in parallel while omitting one interaction class. +#[inline] +pub fn eval_par_with_skip( + kernel: &K, + source_tree: ClusterTreeView<'_, T>, + source_summaries: &[K::SourceSummary], + sources: S, + targets: C, + moments: M, + theta: T, + skip: Skip, + out: [&mut [T]; D], + scratch: &mut EvaluationScratch<'_, [T; D]>, +) -> HierarchicalError +where + K: HierarchicalKernel + Sync, + T: Scalar, + K::TargetGeometry: Copy, + S: SourceCollection, + M: SourceMomentCollection, + C: TargetCollection, +{ + eval_par_optional_skip( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + Some(skip), + out, + scratch, + ) +} + +#[inline] +fn eval_par_optional_skip( + kernel: &K, + source_tree: ClusterTreeView<'_, T>, + source_summaries: &[K::SourceSummary], + sources: S, + targets: C, + moments: M, + theta: T, + skip: Option, + out: [&mut [T]; D], + scratch: &mut EvaluationScratch<'_, [T; D]>, +) -> HierarchicalError where K: HierarchicalKernel + Sync, T: Scalar, @@ -333,6 +488,7 @@ where targets, moments, theta, + skip, out, &mut scratch.contribution[..chunk_count], chunk_size, @@ -352,6 +508,7 @@ fn eval_par_chunks( targets: C, moments: M, theta: T, + skip: Option, out: [&mut [T]; D], scratch_contributions: &mut [[T; D]], chunk_size: usize, @@ -381,6 +538,7 @@ fn eval_par_chunks( targets, moments, theta, + skip, out, &mut chunk_scratch, ); @@ -413,6 +571,7 @@ fn eval_par_chunks( left_targets, moments, theta, + skip, left_out, left_scratch, chunk_size, @@ -428,6 +587,7 @@ fn eval_par_chunks( right_targets, moments, theta, + skip, right_out, right_scratch, chunk_size, diff --git a/src/physics/hierarchical/kernel.rs b/src/physics/hierarchical/kernel.rs index 51976d0..bc33d90 100644 --- a/src/physics/hierarchical/kernel.rs +++ b/src/physics/hierarchical/kernel.rs @@ -1,5 +1,14 @@ 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, +} + /// 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..a3ec3c8 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -31,9 +31,10 @@ 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::{eval, eval_par, eval_par_with_skip, eval_with_skip}; +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..5574c64 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,16 @@ struct TargetSummary { struct MockKernel { _marker: core::marker::PhantomData, + near_calls: AtomicUsize, + far_calls: AtomicUsize, } impl MockKernel { fn new() -> Self { Self { _marker: core::marker::PhantomData, + near_calls: AtomicUsize::new(0), + far_calls: AtomicUsize::new(0), } } } @@ -142,6 +147,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 +158,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 +172,113 @@ impl HierarchicalKernel for MockKernel { } } +#[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, + [&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_with_skip( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + 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_with_skip( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + 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_with_skip( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + 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 hierarchical_error_raw_codes_keep_kernel_slots_first() { assert_eq!(HierarchicalError::KernelError0 as u32, 0); From df72501f21a172003b00f4657a9540d0c70fec5d Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:20:07 -0400 Subject: [PATCH 02/20] [nearfield-separation] Phase 2: expose skip in solver APIs --- cfsem/cfsem.pyi | 20 ++- src/physics/hierarchical/convenience.rs | 211 +++++++++++++++++++++++- src/physics/hierarchical/mod.rs | 11 +- src/python.rs | 187 +++++++++++++++------ test/test_hierarchical.py | 43 +++++ 5 files changed, 415 insertions(+), 57 deletions(-) diff --git a/cfsem/cfsem.pyi b/cfsem/cfsem.pyi index 92149d4..3ad0fdb 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -1,4 +1,4 @@ -from typing import TypeAlias, TypedDict +from typing import Literal, TypeAlias, TypedDict from numpy import float32, float64, int64, uint64 from numpy.typing import NDArray @@ -281,6 +281,7 @@ def flux_density_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical magnetic flux density of dipoles in Cartesian coordinates. @@ -301,6 +302,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -318,6 +321,7 @@ def vector_potential_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical magnetic vector potential of dipoles in Cartesian coordinates. @@ -338,6 +342,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -364,6 +370,7 @@ def flux_density_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for many linear filament segments. @@ -385,6 +392,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -426,6 +435,7 @@ def vector_potential_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for many linear filament segments. @@ -447,6 +457,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -517,6 +529,7 @@ def flux_density_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for a triangle mesh with nodal stream-function values. @@ -546,6 +559,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in @@ -563,6 +578,7 @@ def vector_potential_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, + skip: Literal["near", "far"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for a triangle mesh with nodal stream-function values. @@ -592,6 +608,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, and `None` returns both. Returns: Field component arrays and diagnostics. If `out` is provided, returns `out` in diff --git a/src/physics/hierarchical/convenience.rs b/src/physics/hierarchical/convenience.rs index 52adf7a..77cea19 100644 --- a/src/physics/hierarchical/convenience.rs +++ b/src/physics/hierarchical/convenience.rs @@ -16,8 +16,8 @@ use super::kernels::{ }; use super::{ BuildMethod, ClusterTree, EvaluationScratch, HierarchicalError, HierarchicalKernel, Scalar, - SourceCollection, SourceMomentCollection, SourceNodeSummaries, TargetCollection, eval, - eval_par, scratch_len, scratch_len_par, update_summaries, + Skip, SourceCollection, SourceMomentCollection, SourceNodeSummaries, TargetCollection, eval, + eval_par, eval_par_with_skip, eval_with_skip, scratch_len, scratch_len_par, update_summaries, }; /// Diagnostic information returned by stateless hierarchical solves. @@ -90,6 +90,7 @@ pub fn flux_density_dipole_hierarchical( construction_method, theta, par, + None, out, ) } @@ -142,6 +143,7 @@ pub fn vector_potential_dipole_hierarchical( construction_method, theta, par, + None, out, ) } @@ -195,6 +197,7 @@ pub fn flux_density_linear_filament_hierarchical( construction_method, theta, par, + None, out, ) } @@ -248,6 +251,7 @@ pub fn vector_potential_linear_filament_hierarchical( construction_method, theta, par, + None, out, ) } @@ -308,6 +312,7 @@ pub fn flux_density_triangle_mesh_hierarchical( construction_method, theta, par, + None, out, ) } @@ -368,6 +373,177 @@ pub fn vector_potential_triangle_mesh_hierarchical( construction_method, theta, par, + None, + 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, + 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), + 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, + 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), + 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, + 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), + 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, + 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), + 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, + 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), + 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, + 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), out, ) } @@ -381,6 +557,7 @@ pub(crate) fn one_shot_vec3( construction_method: BuildMethod, theta: T, par: bool, + skip: Option, out: (&mut [T], &mut [T], &mut [T]), ) -> Result, HierarchicalError> where @@ -430,8 +607,31 @@ where contribution: &mut scratch_values, }; let out_components = [out.0, out.1, out.2]; - err = match par { - true => eval_par( + err = match (par, skip) { + (true, Some(skip)) => eval_par_with_skip( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + skip, + out_components, + &mut scratch, + ), + (true, None) => eval_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + out_components, + &mut scratch, + ), + (false, Some(skip)) => eval_with_skip( &kernel, source_tree.as_view(), &source_summaries.node_summaries, @@ -439,10 +639,11 @@ where targets, moments, theta, + skip, out_components, &mut scratch, ), - false => eval( + (false, None) => eval( &kernel, source_tree.as_view(), &source_summaries.node_summaries, diff --git a/src/physics/hierarchical/mod.rs b/src/physics/hierarchical/mod.rs index a3ec3c8..c44abde 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, + 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; diff --git a/src/python.rs b/src/python.rs index 179880e..10d99ce 100644 --- a/src/python.rs +++ b/src/python.rs @@ -621,6 +621,23 @@ fn parse_build_method( } } +/// 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(value) => Err(PyInteropError::ValueError { + msg: format!( + "Unsupported hierarchical skip value: {value}. Expected 'near', 'far', or None." + ), + } + .into()), + } +} + /// Compute accepted source-node levels for hierarchical diagnostic output. fn accepted_levels_diagnostic( kernel: K, @@ -707,7 +724,7 @@ where )) } -#[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<'_>, @@ -736,24 +753,39 @@ 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 skip = parse_hierarchical_skip(skip)?; let (field, 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, + 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, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical dipole flux density", err)) })?; let sources = physics::hierarchical::kernels::DipoleSources::new( @@ -796,7 +828,7 @@ fn flux_density_dipole_hierarchical( ) } -#[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<'_>, @@ -825,24 +857,41 @@ 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 skip = parse_hierarchical_skip(skip)?; let (field, 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, + 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, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical dipole vector potential", err)) })?; let sources = physics::hierarchical::kernels::DipoleSources::new( @@ -885,7 +934,7 @@ fn vector_potential_dipole_hierarchical( ) } -#[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<'_>, @@ -915,6 +964,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")?; @@ -922,19 +972,36 @@ 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 skip = parse_hierarchical_skip(skip)?; let (field, 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, + 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, + out, + ), + } .map_err(|err| py_hierarchical_error("hierarchical linear-filament flux density", err)) })?; let sources = physics::hierarchical::kernels::LinearFilamentSources::new( @@ -971,7 +1038,7 @@ fn flux_density_linear_filament_hierarchical( ) } -#[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<'_>, @@ -1001,6 +1068,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")?; @@ -1008,19 +1076,36 @@ 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 skip = parse_hierarchical_skip(skip)?; let (field, 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, + 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, + out, + ), + } .map_err(|err| { py_hierarchical_error("hierarchical linear-filament vector potential", err) }) @@ -1060,7 +1145,7 @@ fn vector_potential_linear_filament_hierarchical( ) } -#[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<'_>, @@ -1077,6 +1162,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")?; @@ -1084,6 +1170,7 @@ 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 = @@ -1098,6 +1185,7 @@ fn flux_density_triangle_mesh_hierarchical( construction_method, theta, par, + skip, out, ) .map_err(|err| py_hierarchical_error("hierarchical triangle-mesh flux density", err)) @@ -1126,7 +1214,7 @@ fn flux_density_triangle_mesh_hierarchical( ) } -#[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<'_>, @@ -1143,6 +1231,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")?; @@ -1150,6 +1239,7 @@ 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 = @@ -1164,6 +1254,7 @@ fn vector_potential_triangle_mesh_hierarchical( construction_method, theta, par, + skip, out, ) .map_err(|err| { diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 01a8ccb..44f77ac 100644 --- a/test/test_hierarchical.py +++ b/test/test_hierarchical.py @@ -30,6 +30,10 @@ def _assert_returns_output_views(returned, out): assert np.shares_memory(returned_component, out_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 @@ -116,6 +120,45 @@ 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 + ) + far_only = cfsem.vector_potential_dipole_hierarchical( + loc, moment, obs, outer_radius, theta=0.2, par=par, skip="near" + ) + near_only = cfsem.vector_potential_dipole_hierarchical( + loc, moment, obs, outer_radius, theta=0.2, par=par, skip="far" + ) + + _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) + + +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_hierarchical_construction_method_is_exposed(): xyzfil = ( np.array([0.0, 0.5, -0.2]), From 05c6e2437f1434eae5ab0b3b2d9ce0dd87a359e0 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:22:52 -0400 Subject: [PATCH 03/20] [nearfield-separation] Phase 3: expose near-field CSC diagnostics --- cfsem/cfsem.pyi | 3 + src/physics/hierarchical/evaluator.rs | 120 ++++++++++++++ src/physics/hierarchical/mod.rs | 5 +- src/physics/hierarchical/tests.rs | 48 ++++++ src/python.rs | 227 ++++++++++++++++---------- test/test_hierarchical.py | 33 ++++ 6 files changed, 345 insertions(+), 91 deletions(-) diff --git a/cfsem/cfsem.pyi b/cfsem/cfsem.pyi index 3ad0fdb..245026a 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -2,6 +2,7 @@ from typing import Literal, TypeAlias, TypedDict from numpy import float32, float64, int64, uint64 from numpy.typing import NDArray +from scipy.sparse import csc_matrix FloatArray: TypeAlias = NDArray[float64] Float32Array: TypeAlias = NDArray[float32] @@ -46,6 +47,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.""" diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index 32ef7ea..0cf2f5c 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -614,6 +614,126 @@ fn split_output_components( (left, right) } +/// Canonical CSC sparsity for direct source-target interactions selected by a tree walk. +/// +/// 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, +} + +/// Collect accepted levels and the direct near-field CSC pattern in one diagnostic walk. +/// +/// 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, +{ + let err = validate_source_tree_layout(source_tree); + if err != HierarchicalError::Ok { + return Err(err); + } + if !targets.valid_lengths() { + return Err(HierarchicalError::LengthMismatch); + } + if source_summaries.len() < source_tree.n_nodes() { + return Err(HierarchicalError::ScratchTooSmall); + } + + let target_count = targets.len(); + let source_count = source_tree.node_range_count[0] as usize; + let mut accepted_levels = Vec::with_capacity(target_count); + let mut row_indices = Vec::new(); + let mut column_pointers = Vec::with_capacity(target_count + 1); + let mut active = Vec::new(); + column_pointers.push(0); + + for target_id in 0..target_count { + let target = targets.target(target_id); + let mut weighted_level = K::Scalar::ZERO; + let mut represented_sources = K::Scalar::ZERO; + let column_start = row_indices.len(); + + 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_at_node = 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_at_node; + represented_sources = represented_sources + source_count_at_node; + continue; + } + + 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_at_node; + represented_sources = represented_sources + source_count_at_node; + + let start = source_tree.leaf_start[source_node_index] as usize; + let end = start + leaf_count as usize; + row_indices.extend_from_slice(&source_tree.sorted_indices[start..end]); + } 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)); + } + } + + row_indices[column_start..].sort_unstable(); + column_pointers.push(row_indices.len()); + accepted_levels.push(if represented_sources > K::Scalar::ZERO { + weighted_level / represented_sources + } else { + crate::math::cast::(f64::NAN) + }); + } + + Ok(TraversalDiagnostics { + accepted_levels, + near_field_interaction_map: NearFieldInteractionMap { + row_indices, + column_pointers, + source_count, + target_count, + }, + }) +} + /// Compute the source-tree level represented at each target by the terminal traversal nodes. /// /// This is a diagnostic companion to [`eval`]. It mirrors diff --git a/src/physics/hierarchical/mod.rs b/src/physics/hierarchical/mod.rs index c44abde..d72907a 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -38,7 +38,10 @@ pub(crate) use evaluator::eval_dense; pub(crate) use evaluator::{ EvaluationScratch, SourceNodeSummaries, scratch_len, scratch_len_par, update_summaries, }; -pub use evaluator::{eval, eval_par, eval_par_with_skip, eval_with_skip}; +pub use evaluator::{ + NearFieldInteractionMap, TraversalDiagnostics, eval, eval_par, eval_par_with_skip, + eval_with_skip, traversal_diagnostics, +}; pub use kernel::Skip; pub(crate) use kernel::{ BoundedGeometry, BoundedGeometryCollection, HierarchicalError, HierarchicalKernel, diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index 5574c64..aff77fe 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -279,6 +279,54 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( assert_eq!(far_only_par, far_only); } +#[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 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]); + + let mut accepted_levels = vec![0.0; targets.len()]; + assert_eq!( + super::evaluator::accepted_levels( + &kernel, + source_tree.as_view(), + &summaries.node_summaries, + targets.as_slice(), + 0.5, + &mut accepted_levels, + ), + HierarchicalError::Ok + ); + assert_eq!(diagnostics.accepted_levels, accepted_levels); +} + #[test] fn hierarchical_error_raw_codes_keep_kernel_slots_first() { assert_eq!(HierarchicalError::KernelError0 as u32, 0); diff --git a/src/python.rs b/src/python.rs index 10d99ce..2b2a19a 100644 --- a/src/python.rs +++ b/src/python.rs @@ -123,6 +123,7 @@ struct HierarchicalDiagnostics { target_count: usize, source_tree: Option>, accepted_levels: Option>>, + near_field_interaction_map: Option>, } #[pymethods] @@ -164,6 +165,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")] @@ -561,6 +570,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() @@ -574,6 +584,7 @@ fn solve_result_from_field( target_count, source_tree, accepted_levels, + near_field_interaction_map, }, )?; Py::new(py, SolveResult { field, diagnostics }) @@ -638,15 +649,15 @@ fn parse_hierarchical_skip( } } -/// Compute accepted source-node levels for hierarchical diagnostic output. -fn accepted_levels_diagnostic( +/// Compute traversal diagnostics after rebuilding source summaries for the returned tree. +fn traversal_diagnostics_for_python( kernel: K, source_tree: &physics::hierarchical::tree::ClusterTree, sources: S, targets: C, moments: M, theta: f64, -) -> PyResult> +) -> PyResult> where K: physics::hierarchical::kernel::HierarchicalKernel + Sync, S: physics::hierarchical::kernel::SourceCollection + Copy, @@ -656,7 +667,7 @@ where { let mut source_summaries = physics::hierarchical::evaluator::SourceNodeSummaries::::new(source_tree.as_view()); - let mut err = physics::hierarchical::evaluator::update_summaries( + let err = physics::hierarchical::evaluator::update_summaries( &kernel, source_tree.as_view(), sources, @@ -667,22 +678,38 @@ where 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( + physics::hierarchical::traversal_diagnostics( &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)); - } - Ok(out) + ) + .map_err(|err| py_hierarchical_error("hierarchical traversal diagnostic", err)) +} + +/// Convert an owned 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>>); +type OptionalDiagnosticsPy = ( + Option>, + Option>>, + Option>, +); struct HierarchicalDiagnosticRequest<'a, K, S, C, M> { kernel: K, @@ -707,10 +734,10 @@ where C: physics::hierarchical::kernel::TargetCollection, { if !extra_diagnostics { - return Ok((None, None)); + return Ok((None, None, None)); } - let levels = accepted_levels_diagnostic( + let diagnostics = traversal_diagnostics_for_python( request.kernel, request.source_tree, request.sources, @@ -720,7 +747,11 @@ where )?; Ok(( Some(source_tree_diagnostics_object(py, request.source_tree)?), - Some(PyArray1::from_vec(py, levels).unbind()), + Some(PyArray1::from_vec(py, diagnostics.accepted_levels).unbind()), + Some(near_field_interaction_map_object( + py, + diagnostics.near_field_interaction_map, + )?), )) } @@ -804,18 +835,19 @@ fn flux_density_dipole_hierarchical( 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) = + optional_hierarchical_diagnostics( + py, + extra_diagnostics, + HierarchicalDiagnosticRequest { + kernel: physics::hierarchical::kernels::DipoleFluxDensityKernel::::new(), + source_tree: &diagnostics.source_tree, + sources, + targets, + moments, + theta, + }, + )?; solve_result_from_field( py, field, @@ -825,6 +857,7 @@ fn flux_density_dipole_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -910,18 +943,19 @@ fn vector_potential_dipole_hierarchical( 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) = + optional_hierarchical_diagnostics( + py, + extra_diagnostics, + HierarchicalDiagnosticRequest { + kernel: physics::hierarchical::kernels::DipoleVectorPotentialKernel::::new(), + source_tree: &diagnostics.source_tree, + sources, + targets, + moments, + theta, + }, + )?; solve_result_from_field( py, field, @@ -931,6 +965,7 @@ fn vector_potential_dipole_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -1014,18 +1049,20 @@ fn flux_density_linear_filament_hierarchical( 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) = + 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, + }, + )?; solve_result_from_field( py, field, @@ -1035,6 +1072,7 @@ fn flux_density_linear_filament_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -1120,19 +1158,21 @@ fn vector_potential_linear_filament_hierarchical( 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) = + 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, + }, + )?; solve_result_from_field( py, field, @@ -1142,6 +1182,7 @@ fn vector_potential_linear_filament_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -1190,18 +1231,20 @@ fn flux_density_triangle_mesh_hierarchical( ) .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) = + optional_hierarchical_diagnostics( + py, + extra_diagnostics, + HierarchicalDiagnosticRequest { + kernel: + physics::hierarchical::kernels::BoundaryElementFluxDensityKernel::::new(), + source_tree: &diagnostics.source_tree, + sources, + targets, + moments, + theta, + }, + )?; solve_result_from_field( py, field, @@ -1211,6 +1254,7 @@ fn flux_density_triangle_mesh_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } @@ -1261,19 +1305,21 @@ fn vector_potential_triangle_mesh_hierarchical( 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) = + optional_hierarchical_diagnostics( + py, + extra_diagnostics, + HierarchicalDiagnosticRequest { + kernel: + physics::hierarchical::kernels::BoundaryElementVectorPotentialKernel::::new( + ), + source_tree: &diagnostics.source_tree, + sources, + targets, + moments, + theta, + }, + )?; solve_result_from_field( py, field, @@ -1283,6 +1329,7 @@ fn vector_potential_triangle_mesh_hierarchical( diagnostics.target_count, source_tree, accepted_levels, + near_field_interaction_map, ) } diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 44f77ac..72e7506 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 @@ -46,6 +47,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(): @@ -159,6 +164,34 @@ def test_hierarchical_skip_rejects_unknown_value(): ) +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_hierarchical_construction_method_is_exposed(): xyzfil = ( np.array([0.0, 0.5, -0.2]), From 5704b47d6684093a6e4c43e1fa626a9674442847 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:25:26 -0400 Subject: [PATCH 04/20] [nearfield-separation] Phase 4: add sparse GL3 inductance core --- src/physics/linear_filament.rs | 338 ++++++++++++++++++++++++++++++--- 1 file changed, 310 insertions(+), 28 deletions(-) diff --git a/src/physics/linear_filament.rs b/src/physics/linear_filament.rs index 26b7f4c..fd70e57 100644 --- a/src/physics/linear_filament.rs +++ b/src/physics/linear_filament.rs @@ -1,7 +1,9 @@ //! Magnetics calculations for piecewise-linear current filaments. use rayon::{ - iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator}, + iter::{ + IndexedParallelIterator, IntoParallelIterator, IntoParallelRefMutIterator, ParallelIterator, + }, slice::{ParallelSlice, ParallelSliceMut}, }; @@ -17,6 +19,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 +103,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 +335,185 @@ 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 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 { + let tgt_start = ( + xyzfil_tgt.0[target], + xyzfil_tgt.1[target], + xyzfil_tgt.2[target], + ); + let tgt_delta = ( + dlxyzfil_tgt.0[target], + dlxyzfil_tgt.1[target], + dlxyzfil_tgt.2[target], + ); + for entry in column_pointers[target]..column_pointers[target + 1] { + let source = row_indices[entry]; + out[entry] = 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], + tgt_start, + tgt_delta, + ); + } + } + + 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, + )?; + + out.par_iter_mut().enumerate().for_each(|(entry, value)| { + let target = column_pointers.partition_point(|&pointer| pointer <= entry) - 1; + let source = row_indices[entry]; + *value = 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], + ), + ); + }); + + Ok(()) +} + /// Biot-Savart calculation for B-field contribution from many current filament /// segments to many observation points. /// @@ -2848,4 +3043,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" + ); + } } From 415e03a57771d6448dea5f524e99f34a874bdbd9 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:27:45 -0400 Subject: [PATCH 05/20] [nearfield-separation] Phase 5: expose sparse inductance in Python --- cfsem/__init__.py | 2 + cfsem/bindings.py | 73 ++++++++++++++++++++++++++++ cfsem/cfsem.pyi | 10 ++++ src/python.rs | 62 ++++++++++++++++++++++++ test/test_electromagnetics.py | 91 ++++++++++++++++++++++++++++++++--- 5 files changed, 232 insertions(+), 6 deletions(-) diff --git a/cfsem/__init__.py b/cfsem/__init__.py index 618037a..ce3648b 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..01318e6 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_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, + 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``. + 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): + raise TypeError("interaction_map must be a scipy.sparse.csc_matrix") + 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 245026a..7677199 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -506,6 +506,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( diff --git a/src/python.rs b/src/python.rs index 2b2a19a..28fb88f 100644 --- a/src/python.rs +++ b/src/python.rs @@ -3409,6 +3409,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 = row_indices.as_slice()?; + let column_pointers = column_pointers.as_slice()?; + 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( @@ -4546,6 +4604,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..264e33a 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,86 @@ def test_inductance_linear_filaments_matrix_contracts_to_vector(par): ) +@mark.parametrize("par", [True, False]) +def test_inductance_linear_filaments_sparse_preserves_csc_pattern(par): + 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 = sparse.csc_matrix( + (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 + + +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"): + 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): From ea53571fe4995c619c3211c9f7cc90a3d7a42b9d Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:28:42 -0400 Subject: [PATCH 06/20] [nearfield-separation] Phase 6: verify frozen-far workflow --- test/test_hierarchical.py | 165 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 159 insertions(+), 6 deletions(-) diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 72e7506..97e9f6f 100644 --- a/test/test_hierarchical.py +++ b/test/test_hierarchical.py @@ -31,6 +31,11 @@ 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)) @@ -83,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( @@ -90,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(): @@ -139,18 +149,41 @@ def test_hierarchical_skip_decomposes_near_and_far_fields(par): outer_radius = np.full(6, 0.05) full = cfsem.vector_potential_dipole_hierarchical( - loc, moment, obs, outer_radius, theta=0.2, par=par + 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" + 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" + loc, + moment, + obs, + outer_radius, + theta=0.2, + par=par, + skip="far", + 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) + for filtered in (far_only, near_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, + ) def test_hierarchical_skip_rejects_unknown_value(): @@ -159,9 +192,7 @@ def test_hierarchical_skip_rejects_unknown_value(): 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" - ) + 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(): @@ -192,6 +223,128 @@ def test_near_field_interaction_map_uses_original_source_rows(): np.testing.assert_array_equal(interaction_map.data, np.ones(2, dtype=bool)) +def test_all_hierarchical_methods_accept_skip(): + 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="near"), + cfsem.vector_potential_dipole_hierarchical(loc, moment, obs, radius, theta=0.0, skip="near"), + cfsem.flux_density_linear_filament_hierarchical( + obs, xyzfil, dlxyzfil, current, radius, theta=0.0, skip="near" + ), + cfsem.vector_potential_linear_filament_hierarchical( + obs, xyzfil, dlxyzfil, current, radius, theta=0.0, skip="near" + ), + cfsem.flux_density_triangle_mesh_hierarchical( + obs_rows, nodes, triangles, stream_function, theta=0.0, skip="near" + ), + cfsem.vector_potential_triangle_mesh_hierarchical( + obs_rows, nodes, triangles, stream_function, theta=0.0, skip="near" + ), + ) + for result in results: + _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, + extra_diagnostics=True, + ) + 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)) + + interaction_map = full.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]), From 3b64c03070237eafd4ba6b1afbada3eb4e03079c Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:37:02 -0400 Subject: [PATCH 07/20] [nearfield-separation] Phase 7: document and benchmark workflow --- CHANGELOG.md | 6 +++ benches/linear_filament.rs | 74 ++++++++++++++++++++++++++++++++- docs/python/boundary_element.md | 4 ++ docs/python/dipole.md | 5 +++ docs/python/filament.md | 63 ++++++++++++++++++++++++++++ test/test_electromagnetics.py | 11 +++++ 6 files changed, 162 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1a1e2b2..478e8cd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,9 @@ ## 12.0.0 2026-08-11 * Rust + * Add evaluator-level near/far filtering for hierarchical solves through required `Skip::{Near, Far}` filtered APIs + * Add canonical `(nsrc, ntgt)` CSC near-field traversal diagnostics using original source indices + * Add serial and parallel sparse CSC linear-filament inductance evaluation with three-point target-segment quadrature * Fix BEM triangle inductance for touching and near high-aspect-ratio triangles with an exact uniform-triangle source potential * Near-field defers to fixed D5 target integration at subdivision depth 2 * Bound near/self work at 28 exact-potential observations per directed pair @@ -11,6 +14,9 @@ * Evaluate direct triangle B fields from the analytic potential gradient and define as zero directly on the surface * Protects zero self-force per Newton's third law * Python + * Add optional `skip="near"` and `skip="far"` filtering to all hierarchical field solvers + * Add `HierarchicalDiagnostics.near_field_interaction_map` when extra diagnostics are requested + * Add `inductance_linear_filaments_sparse` with exact CSC sparsity preservation * Add an absolute aspect-67 annular stored-energy regression test * !Remove `quad` from triangle-mesh direct, mapping, and hierarchical B-field and vector potential functions diff --git a/benches/linear_filament.rs b/benches/linear_filament.rs index 21931cf..cdee3b0 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::*; @@ -366,6 +367,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 +444,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/docs/python/boundary_element.md b/docs/python/boundary_element.md index f5cbbf2..dd56b57 100644 --- a/docs/python/boundary_element.md +++ b/docs/python/boundary_element.md @@ -1,5 +1,9 @@ # 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. Extra diagnostics include the resulting +direct-interaction pattern as a canonical `(ntri, ntgt)` SciPy CSC matrix. + ## 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..99474e8 100644 --- a/docs/python/dipole.md +++ b/docs/python/dipole.md @@ -1,5 +1,10 @@ # Dipole +Hierarchical dipole solvers accept an optional interaction filter. `skip="near"` returns the +far-only result, `skip="far"` returns the direct near-only result, 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. + ## Fields ::: cfsem.flux_density_dipole diff --git a/docs/python/filament.md b/docs/python/filament.md index a9a03db..c2dfe68 100644 --- a/docs/python/filament.md +++ b/docs/python/filament.md @@ -1,5 +1,11 @@ # Linear Filament +Hierarchical field solvers accept `skip=None`, `skip="near"`, or `skip="far"`. 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. A full solve is therefore the sum of +those two filtered 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 +16,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/test/test_electromagnetics.py b/test/test_electromagnetics.py index 264e33a..0fbac4e 100644 --- a/test/test_electromagnetics.py +++ b/test/test_electromagnetics.py @@ -1113,6 +1113,17 @@ def test_inductance_linear_filaments_sparse_preserves_csc_pattern(par): 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)) From 754953d35838e27fd9f201bf567086d5724a2406 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 15:58:55 -0400 Subject: [PATCH 08/20] [nearfield-separation] Phase 8: add diagnostics-only skip mode Add Skip::Both zero-field fast paths across Rust and Python while retaining a single traversal when diagnostics are requested. Plan: /Users/jlogan/git/cfsem-py/nearfield_separation_plan.md Co-Authored-By: Codex --- CHANGELOG.md | 4 +- cfsem/cfsem.pyi | 24 ++-- docs/python/boundary_element.md | 6 +- docs/python/dipole.md | 7 +- docs/python/filament.md | 12 +- src/physics/hierarchical/convenience.rs | 147 +++++++++++--------- src/physics/hierarchical/evaluator.rs | 20 ++- src/physics/hierarchical/kernel.rs | 2 + src/physics/hierarchical/tests.rs | 176 ++++++++++++++++++++++++ src/python.rs | 3 +- test/test_hierarchical.py | 47 +++++-- 11 files changed, 345 insertions(+), 103 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 478e8cd..fdf6f93 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,7 +3,7 @@ ## 12.0.0 2026-08-11 * Rust - * Add evaluator-level near/far filtering for hierarchical solves through required `Skip::{Near, Far}` filtered APIs + * Add evaluator-level near/far filtering for hierarchical solves through required `Skip::{Near, Far, Both}` filtered APIs * Add canonical `(nsrc, ntgt)` CSC near-field traversal diagnostics using original source indices * Add serial and parallel sparse CSC linear-filament inductance evaluation with three-point target-segment quadrature * Fix BEM triangle inductance for touching and near high-aspect-ratio triangles with an exact uniform-triangle source potential @@ -14,7 +14,7 @@ * Evaluate direct triangle B fields from the analytic potential gradient and define as zero directly on the surface * Protects zero self-force per Newton's third law * Python - * Add optional `skip="near"` and `skip="far"` filtering to all hierarchical field solvers + * Add optional `skip="near"`, `skip="far"`, and diagnostics-only `skip="both"` filtering to all hierarchical field solvers * Add `HierarchicalDiagnostics.near_field_interaction_map` when extra diagnostics are requested * Add `inductance_linear_filaments_sparse` with exact CSC sparsity preservation * Add an absolute aspect-67 annular stored-energy regression test diff --git a/cfsem/cfsem.pyi b/cfsem/cfsem.pyi index 7677199..cd222a6 100644 --- a/cfsem/cfsem.pyi +++ b/cfsem/cfsem.pyi @@ -284,7 +284,7 @@ def flux_density_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical magnetic flux density of dipoles in Cartesian coordinates. @@ -306,7 +306,7 @@ def flux_density_dipole_hierarchical( 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, and `None` returns both. + 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 @@ -324,7 +324,7 @@ def vector_potential_dipole_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical magnetic vector potential of dipoles in Cartesian coordinates. @@ -346,7 +346,7 @@ def vector_potential_dipole_hierarchical( 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, and `None` returns both. + 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 @@ -373,7 +373,7 @@ def flux_density_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for many linear filament segments. @@ -396,7 +396,7 @@ def flux_density_linear_filament_hierarchical( 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, and `None` returns both. + 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 @@ -438,7 +438,7 @@ def vector_potential_linear_filament_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for many linear filament segments. @@ -461,7 +461,7 @@ def vector_potential_linear_filament_hierarchical( 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, and `None` returns both. + 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 @@ -542,7 +542,7 @@ def flux_density_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical B-field calculation for a triangle mesh with nodal stream-function values. @@ -573,7 +573,7 @@ def flux_density_triangle_mesh_hierarchical( 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, and `None` returns both. + 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 @@ -591,7 +591,7 @@ def vector_potential_triangle_mesh_hierarchical( par: bool = True, out: ArrayTriple | None = None, extra_diagnostics: bool = False, - skip: Literal["near", "far"] | None = None, + skip: Literal["near", "far", "both"] | None = None, ) -> SolveResult: """Hierarchical A-field calculation for a triangle mesh with nodal stream-function values. @@ -622,7 +622,7 @@ def vector_potential_triangle_mesh_hierarchical( 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, and `None` returns both. + 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 dd56b57..2ab9f49 100644 --- a/docs/python/boundary_element.md +++ b/docs/python/boundary_element.md @@ -1,8 +1,10 @@ # 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. Extra diagnostics include the resulting -direct-interaction pattern as a canonical `(ntri, ntgt)` SciPy CSC matrix. +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 diff --git a/docs/python/dipole.md b/docs/python/dipole.md index 99474e8..a7b543e 100644 --- a/docs/python/dipole.md +++ b/docs/python/dipole.md @@ -1,9 +1,10 @@ # Dipole Hierarchical dipole solvers accept an optional interaction filter. `skip="near"` returns the -far-only result, `skip="far"` returns the direct near-only result, 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. +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 diff --git a/docs/python/filament.md b/docs/python/filament.md index c2dfe68..65d3362 100644 --- a/docs/python/filament.md +++ b/docs/python/filament.md @@ -1,10 +1,12 @@ # Linear Filament -Hierarchical field solvers accept `skip=None`, `skip="near"`, or `skip="far"`. 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. A full solve is therefore the sum of -those two filtered solves, up to floating-point roundoff. Filtering happens in the common evaluator, -so skipped kernel calculations are not performed. +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 diff --git a/src/physics/hierarchical/convenience.rs b/src/physics/hierarchical/convenience.rs index 77cea19..5989701 100644 --- a/src/physics/hierarchical/convenience.rs +++ b/src/physics/hierarchical/convenience.rs @@ -584,79 +584,90 @@ 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, skip) { - (true, Some(skip)) => eval_par_with_skip( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - sources, - targets, - moments, - theta, - skip, - out_components, - &mut scratch, - ), - (true, None) => eval_par( + let source_summaries = if skip == Some(Skip::Both) { + 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, - ), - (false, Some(skip)) => eval_with_skip( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - sources, - targets, - moments, - theta, - skip, - out_components, - &mut scratch, - ), - (false, None) => eval( - &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 let Some(source_summaries) = 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]; + let err = match (par, skip) { + (true, Some(skip)) => eval_par_with_skip( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + skip, + out_components, + &mut scratch, + ), + (true, None) => eval_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + out_components, + &mut scratch, + ), + (false, Some(skip)) => eval_with_skip( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + skip, + out_components, + &mut scratch, + ), + (false, None) => eval( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + sources, + targets, + moments, + theta, + 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(); diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index 0cf2f5c..fc685a8 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -130,7 +130,8 @@ where /// /// The source tree is still traversed normally so acceptance decisions do not /// change. [`Skip::Near`] returns accepted far-summary contributions only; -/// [`Skip::Far`] returns direct leaf contributions only. +/// [`Skip::Far`] returns direct leaf contributions only. [`Skip::Both`] zeroes the output without +/// target summarization or source-tree traversal. #[inline] pub fn eval_with_skip( kernel: &K, @@ -241,6 +242,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() || scratch.contribution.is_empty() { return HierarchicalError::ScratchTooSmall; } @@ -391,7 +398,10 @@ where ) } -/// Evaluate vector-valued targets in parallel while omitting one interaction class. +/// Evaluate vector-valued targets in parallel while omitting interaction classes. +/// +/// [`Skip::Both`] zeroes the output without target summarization, source-tree traversal, or +/// parallel scratch use. #[inline] pub fn eval_par_with_skip( kernel: &K, @@ -466,6 +476,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; } diff --git a/src/physics/hierarchical/kernel.rs b/src/physics/hierarchical/kernel.rs index bc33d90..2c60aa6 100644 --- a/src/physics/hierarchical/kernel.rs +++ b/src/physics/hierarchical/kernel.rs @@ -7,6 +7,8 @@ pub enum Skip { 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. diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index aff77fe..386a2a1 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -44,6 +44,8 @@ struct TargetSummary { struct MockKernel { _marker: core::marker::PhantomData, + target_summary_calls: AtomicUsize, + accept_calls: AtomicUsize, near_calls: AtomicUsize, far_calls: AtomicUsize, } @@ -52,6 +54,8 @@ 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), } @@ -124,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; @@ -140,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, @@ -172,6 +188,87 @@ 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(); @@ -277,6 +374,85 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( 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); + + kernel.target_summary_calls.store(0, Ordering::Relaxed); + kernel.accept_calls.store(0, Ordering::Relaxed); + kernel.near_calls.store(0, Ordering::Relaxed); + kernel.far_calls.store(0, Ordering::Relaxed); + let mut skipped_both = [f64::NAN]; + let mut empty_contribution = []; + let mut empty_scratch = EvaluationScratch { + contribution: &mut empty_contribution, + }; + assert_eq!( + super::eval_with_skip( + &kernel, + source_tree.as_view(), + &[], + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Skip::Both, + [&mut skipped_both], + &mut empty_scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(skipped_both, [0.0]); + assert_eq!(kernel.target_summary_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.accept_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 0); + + skipped_both.fill(f64::NAN); + assert_eq!( + super::eval_par_with_skip( + &kernel, + source_tree.as_view(), + &[], + sources.as_slice(), + targets.as_slice(), + &moments, + 0.5, + Skip::Both, + [&mut skipped_both], + &mut empty_scratch, + ), + HierarchicalError::Ok + ); + assert_eq!(skipped_both, [0.0]); + assert_eq!(kernel.target_summary_calls.load(Ordering::Relaxed), 0); + assert_eq!(kernel.accept_calls.load(Ordering::Relaxed), 0); +} + +#[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), + (&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); } #[test] diff --git a/src/python.rs b/src/python.rs index 28fb88f..17704fe 100644 --- a/src/python.rs +++ b/src/python.rs @@ -640,9 +640,10 @@ fn parse_hierarchical_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', or None." + "Unsupported hierarchical skip value: {value}. Expected 'near', 'far', 'both', or None." ), } .into()), diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 97e9f6f..1d7e62b 100644 --- a/test/test_hierarchical.py +++ b/test/test_hierarchical.py @@ -171,11 +171,22 @@ def test_hierarchical_skip_decomposes_near_and_far_fields(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) - for filtered in (far_only, near_only): + _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, @@ -185,6 +196,12 @@ def test_hierarchical_skip_decomposes_near_and_far_fields(par): 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])) @@ -223,7 +240,7 @@ def test_near_field_interaction_map_uses_original_source_rows(): np.testing.assert_array_equal(interaction_map.data, np.ones(2, dtype=bool)) -def test_all_hierarchical_methods_accept_skip(): +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])) @@ -237,25 +254,39 @@ def test_all_hierarchical_methods_accept_skip(): obs_rows = np.array([[0.2, 0.2, 0.5]]) results = ( - cfsem.flux_density_dipole_hierarchical(loc, moment, obs, radius, theta=0.0, skip="near"), - cfsem.vector_potential_dipole_hierarchical(loc, moment, obs, radius, theta=0.0, skip="near"), + 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="near" + 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="near" + 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="near" + 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="near" + 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 From ead0fd864b12d97f89cd4a25dbaf80459ad9bcd1 Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 16:14:11 -0400 Subject: [PATCH 09/20] [nearfield-separation] Test diagnostics-only inductance workflow --- test/test_hierarchical.py | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/test/test_hierarchical.py b/test/test_hierarchical.py index 1d7e62b..7fe905d 100644 --- a/test/test_hierarchical.py +++ b/test/test_hierarchical.py @@ -309,7 +309,17 @@ def test_near_field_map_drives_sparse_filament_inductance(par): 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, @@ -332,8 +342,9 @@ def test_near_field_map_drives_sparse_filament_inductance(par): skip="far", ) _assert_vec_close(full, _add_vec3(far_only.field, near_only.field)) + _assert_vec_zero(diagnostics_only) - interaction_map = full.diagnostics.near_field_interaction_map + interaction_map = diagnostics_only.diagnostics.near_field_interaction_map assert 0 < interaction_map.nnz < nsegment * nsegment sparse_inductance = cfsem.inductance_linear_filaments_sparse( xyzfil, From ec57350f1aec9ff99b36fc4e4e092c322312b3fa Mon Sep 17 00:00:00 2001 From: James Logan Date: Thu, 13 Aug 2026 16:18:17 -0400 Subject: [PATCH 10/20] update version and changelog --- CHANGELOG.md | 9 +++++++++ Cargo.lock | 2 +- Cargo.toml | 2 +- 3 files changed, 11 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fdf6f93..7c1d204 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Changelog +## 12.1.0 2026-08-13 + +* 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 + * Add sparse inductance matrix for linear filaments using CSC interaction map +* Python + * Plumb in bindings to new `skip` option and sparse inductance matrix + ## 12.0.0 2026-08-11 * Rust diff --git a/Cargo.lock b/Cargo.lock index cd09805..98b8c94 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -128,7 +128,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "cfsem" -version = "12.0.0" +version = "12.1.0" dependencies = [ "criterion", "faer", diff --git a/Cargo.toml b/Cargo.toml index 6dd1e19..d52f2f7 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "12.0.0" +version = "12.1.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" From 1cfa73531d491b6d45f216b54c4edfd81b3aa333 Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 14:20:00 -0400 Subject: [PATCH 11/20] update version --- Cargo.lock | 2 +- Cargo.toml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 69b5481..357d934 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -128,7 +128,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "cfsem" -version = "12.1.0" +version = "12.2.0" dependencies = [ "criterion", "faer", diff --git a/Cargo.toml b/Cargo.toml index 2a77c1d..87cce64 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "12.1.0" +version = "12.2.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" From 0a31dd85d1de3a677d824ecd06303c54d416a6bc Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 15:57:04 -0400 Subject: [PATCH 12/20] update changelog --- CHANGELOG.md | 6 ------ 1 file changed, 6 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 98bcfd1..ad1053f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -19,9 +19,6 @@ ## 12.0.0 2026-08-11 * Rust - * Add evaluator-level near/far filtering for hierarchical solves through required `Skip::{Near, Far, Both}` filtered APIs - * Add canonical `(nsrc, ntgt)` CSC near-field traversal diagnostics using original source indices - * Add serial and parallel sparse CSC linear-filament inductance evaluation with three-point target-segment quadrature * Fix BEM triangle inductance for touching and near high-aspect-ratio triangles with an exact uniform-triangle source potential * Near-field defers to fixed D5 target integration at subdivision depth 2 * Bound near/self work at 28 exact-potential observations per directed pair @@ -30,9 +27,6 @@ * Evaluate direct triangle B fields from the analytic potential gradient and define as zero directly on the surface * Protects zero self-force per Newton's third law * Python - * Add optional `skip="near"`, `skip="far"`, and diagnostics-only `skip="both"` filtering to all hierarchical field solvers - * Add `HierarchicalDiagnostics.near_field_interaction_map` when extra diagnostics are requested - * Add `inductance_linear_filaments_sparse` with exact CSC sparsity preservation * Add an absolute aspect-67 annular stored-energy regression test * !Remove `quad` from triangle-mesh direct, mapping, and hierarchical B-field and vector potential functions From e6b0fc27652d51a844a3e66aa7e8d7d161b749f8 Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 15:57:30 -0400 Subject: [PATCH 13/20] check CSC index values for 32 bit compat --- src/python.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/python.rs b/src/python.rs index 9b27c28..d263aee 100644 --- a/src/python.rs +++ b/src/python.rs @@ -3438,8 +3438,8 @@ fn inductance_linear_filaments_sparse_csc( PyReadonlyArray1, ), wire_radius_src: PyReadonlyArray1, - row_indices: PyReadonlyArray1, - column_pointers: PyReadonlyArray1, + row_indices: PyReadonlyArray1, + column_pointers: PyReadonlyArray1, par: bool, ) -> PyResult>> { _3tup_slice_ro!(xyzfil_tgt); @@ -3447,8 +3447,8 @@ fn inductance_linear_filaments_sparse_csc( _3tup_slice_ro!(xyzfil_src); _3tup_slice_ro!(dlxyzfil_src); let wire_radius_src = wire_radius_src.as_slice()?; - let row_indices = row_indices.as_slice()?; - let column_pointers = column_pointers.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 { @@ -3461,8 +3461,8 @@ fn inductance_linear_filaments_sparse_csc( xyzfil_src, dlxyzfil_src, wire_radius_src, - row_indices, - column_pointers, + &row_indices, + &column_pointers, &mut out, ) .map_err(|msg| PyInteropError::DimensionalityError { From c8b3f096a103ca6f83830e2678e03abefe4fac0a Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 16:15:53 -0400 Subject: [PATCH 14/20] consolidate traversal implementations and parallelize diagnostics --- src/physics/hierarchical/evaluator.rs | 481 +++++++++++++++++--------- src/physics/hierarchical/mod.rs | 2 +- src/physics/hierarchical/tests.rs | 24 +- src/python.rs | 34 +- 4 files changed, 361 insertions(+), 180 deletions(-) diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index fc685a8..2dfd84c 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -2,6 +2,7 @@ use super::{ BoundedGeometry, ClusterTreeView, HierarchicalError, HierarchicalKernel, Scalar, Skip, SourceCollection, SourceMomentCollection, TargetCollection, }; +use rayon::prelude::*; use std::sync::atomic::{AtomicU32, Ordering}; /// CPU-owned source summary storage. @@ -285,10 +286,119 @@ 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> +where + K: HierarchicalKernel, +{ + kernel: &'a K, + sources: S, + target: K::TargetGeometry, + moments: M, + skip: Option, + out: &'a mut K::Output, + contribution: &'a mut K::Output, + target_summary: &'a K::TargetSummary, +} + +impl TraversalVisitor for EvaluationTraversalVisitor<'_, K, S, M> +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 self.skip != Some(Skip::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 self.skip == Some(Skip::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( kernel: &K, @@ -302,7 +412,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 @@ -318,41 +428,25 @@ 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) { - if skip != Some(Skip::Far) { - 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 { - if skip == Some(Skip::Near) { - continue; - } - 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, + skip, + out, + contribution, + target_summary, + }; + traverse_source_tree( + kernel, + source_tree, + source_summaries, + &target, + theta, + active, + &mut visitor, + ); HierarchicalError::Ok } @@ -655,21 +749,126 @@ pub struct TraversalDiagnostics { pub near_field_interaction_map: NearFieldInteractionMap, } -/// Collect accepted levels and the direct near-field CSC pattern in one diagnostic walk. -/// -/// 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( +/// 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, -) -> Result, 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 { @@ -681,65 +880,32 @@ where if source_summaries.len() < source_tree.n_nodes() { return Err(HierarchicalError::ScratchTooSmall); } + Ok(()) +} - let target_count = targets.len(); - let source_count = source_tree.node_range_count[0] as usize; +/// 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::new(); + let mut row_indices = Vec::with_capacity(entry_count); let mut column_pointers = Vec::with_capacity(target_count + 1); - let mut active = Vec::new(); column_pointers.push(0); - for target_id in 0..target_count { - let target = targets.target(target_id); - let mut weighted_level = K::Scalar::ZERO; - let mut represented_sources = K::Scalar::ZERO; - let column_start = row_indices.len(); - - 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_at_node = 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_at_node; - represented_sources = represented_sources + source_count_at_node; - continue; - } - - 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_at_node; - represented_sources = represented_sources + source_count_at_node; - - let start = source_tree.leaf_start[source_node_index] as usize; - let end = start + leaf_count as usize; - row_indices.extend_from_slice(&source_tree.sorted_indices[start..end]); - } 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); } - - row_indices[column_start..].sort_unstable(); - column_pointers.push(row_indices.len()); - accepted_levels.push(if represented_sources > K::Scalar::ZERO { - weighted_level / represented_sources - } else { - crate::math::cast::(f64::NAN) - }); } - Ok(TraversalDiagnostics { + 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, @@ -747,84 +913,85 @@ where source_count, target_count, }, - }) + } } -/// Compute the source-tree level represented at each target by the terminal traversal nodes. +/// Collect accepted levels and the direct near-field CSC pattern. /// -/// 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( +/// 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, - out: &mut [K::Scalar], -) -> HierarchicalError +) -> Result, HierarchicalError> where K: HierarchicalKernel, K::TargetGeometry: Copy, C: TargetCollection, { - let err = validate_source_tree_layout(source_tree); - if err != HierarchicalError::Ok { - return err; - } - if targets.len() != out.len() || !targets.valid_lengths() { - return HierarchicalError::LengthMismatch; - } - if source_summaries.len() < source_tree.n_nodes() { - return HierarchicalError::ScratchTooSmall; - } - - 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; - } - - 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)); - } - } + 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, + )) +} - out[target_id] = if represented_sources > K::Scalar::ZERO { - weighted_level / represented_sources - } else { - crate::math::cast::(f64::NAN) - }; +/// 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, + )); } - HierarchicalError::Ok + 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/mod.rs b/src/physics/hierarchical/mod.rs index d72907a..1700b2e 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -40,7 +40,7 @@ pub(crate) use evaluator::{ }; pub use evaluator::{ NearFieldInteractionMap, TraversalDiagnostics, eval, eval_par, eval_par_with_skip, - eval_with_skip, traversal_diagnostics, + eval_with_skip, traversal_diagnostics, traversal_diagnostics_par, }; pub use kernel::Skip; pub(crate) use kernel::{ diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index 386a2a1..4702709 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -482,25 +482,21 @@ fn traversal_diagnostics_returns_canonical_near_field_csc_pattern() { 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]); - - let mut accepted_levels = vec![0.0; targets.len()]; - assert_eq!( - super::evaluator::accepted_levels( - &kernel, - source_tree.as_view(), - &summaries.node_summaries, - targets.as_slice(), - 0.5, - &mut accepted_levels, - ), - HierarchicalError::Ok - ); - assert_eq!(diagnostics.accepted_levels, accepted_levels); } #[test] diff --git a/src/python.rs b/src/python.rs index d263aee..a28aa20 100644 --- a/src/python.rs +++ b/src/python.rs @@ -662,6 +662,7 @@ fn traversal_diagnostics_for_python( targets: C, moments: M, theta: f64, + par: bool, ) -> PyResult> where K: physics::hierarchical::kernel::HierarchicalKernel + Sync, @@ -683,14 +684,23 @@ where return Err(py_hierarchical_error("source summary update", err)); } - physics::hierarchical::traversal_diagnostics( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - targets, - theta, - ) - .map_err(|err| py_hierarchical_error("hierarchical traversal diagnostic", err)) + let diagnostics = match par { + true => physics::hierarchical::traversal_diagnostics_par( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + targets, + theta, + ), + false => physics::hierarchical::traversal_diagnostics( + &kernel, + source_tree.as_view(), + &source_summaries.node_summaries, + targets, + theta, + ), + }; + diagnostics.map_err(|err| py_hierarchical_error("hierarchical traversal diagnostic", err)) } /// Convert an owned CSC interaction pattern into a SciPy sparse matrix. @@ -723,6 +733,7 @@ struct HierarchicalDiagnosticRequest<'a, K, S, C, M> { targets: C, moments: M, theta: f64, + par: bool, } /// Return diagnostics only when requested by the Python caller. @@ -749,6 +760,7 @@ where request.targets, request.moments, request.theta, + request.par, )?; Ok(( Some(source_tree_diagnostics_object(py, request.source_tree)?), @@ -851,6 +863,7 @@ fn flux_density_dipole_hierarchical( targets, moments, theta, + par, }, )?; solve_result_from_field( @@ -959,6 +972,7 @@ fn vector_potential_dipole_hierarchical( targets, moments, theta, + par, }, )?; solve_result_from_field( @@ -1066,6 +1080,7 @@ fn flux_density_linear_filament_hierarchical( targets, moments: ifil.as_slice(), theta, + par, }, )?; solve_result_from_field( @@ -1176,6 +1191,7 @@ fn vector_potential_linear_filament_hierarchical( targets, moments: ifil.as_slice(), theta, + par, }, )?; solve_result_from_field( @@ -1248,6 +1264,7 @@ fn flux_density_triangle_mesh_hierarchical( targets, moments, theta, + par, }, )?; solve_result_from_field( @@ -1323,6 +1340,7 @@ fn vector_potential_triangle_mesh_hierarchical( targets, moments, theta, + par, }, )?; solve_result_from_field( From 1cbc37cdfa5581bbdf2e794ffacd1e9d07cfb07d Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 16:22:47 -0400 Subject: [PATCH 15/20] consolidate sparse near-field inductance matrix impl --- src/physics/linear_filament.rs | 124 ++++++++++++++++++--------------- 1 file changed, 69 insertions(+), 55 deletions(-) diff --git a/src/physics/linear_filament.rs b/src/physics/linear_filament.rs index fd70e57..e317f71 100644 --- a/src/physics/linear_filament.rs +++ b/src/physics/linear_filament.rs @@ -1,9 +1,7 @@ //! Magnetics calculations for piecewise-linear current filaments. use rayon::{ - iter::{ - IndexedParallelIterator, IntoParallelIterator, IntoParallelRefMutIterator, ParallelIterator, - }, + iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator}, slice::{ParallelSlice, ParallelSliceMut}, }; @@ -399,6 +397,42 @@ fn validate_sparse_inductance_inputs( 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 @@ -428,32 +462,15 @@ pub fn inductance_linear_filaments_sparse_csc( )?; for target in 0..ntgt { - let tgt_start = ( - xyzfil_tgt.0[target], - xyzfil_tgt.1[target], - xyzfil_tgt.2[target], - ); - let tgt_delta = ( - dlxyzfil_tgt.0[target], - dlxyzfil_tgt.1[target], - dlxyzfil_tgt.2[target], - ); for entry in column_pointers[target]..column_pointers[target + 1] { - let source = row_indices[entry]; - out[entry] = 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], - tgt_start, - tgt_delta, + out[entry] = inductance_linear_filaments_sparse_entry( + xyzfil_tgt, + dlxyzfil_tgt, + xyzfil_src, + dlxyzfil_src, + wire_radius_src, + target, + row_indices[entry], ); } } @@ -483,33 +500,30 @@ pub fn inductance_linear_filaments_sparse_csc_par( out, )?; - out.par_iter_mut().enumerate().for_each(|(entry, value)| { - let target = column_pointers.partition_point(|&pointer| pointer <= entry) - 1; - let source = row_indices[entry]; - *value = 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], - ), - ); - }); + // 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(()) } From e71d97cd9b582b987304ac6424311476991584e2 Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 16:39:25 -0400 Subject: [PATCH 16/20] consolidate Skip options --- src/physics/hierarchical/convenience.rs | 14 +- src/physics/hierarchical/evaluator.rs | 349 ++++++++++-------------- src/physics/hierarchical/mod.rs | 4 +- src/physics/hierarchical/tests.rs | 65 +---- 4 files changed, 158 insertions(+), 274 deletions(-) diff --git a/src/physics/hierarchical/convenience.rs b/src/physics/hierarchical/convenience.rs index 5989701..8aa9709 100644 --- a/src/physics/hierarchical/convenience.rs +++ b/src/physics/hierarchical/convenience.rs @@ -17,7 +17,7 @@ use super::kernels::{ use super::{ BuildMethod, ClusterTree, EvaluationScratch, HierarchicalError, HierarchicalKernel, Scalar, Skip, SourceCollection, SourceMomentCollection, SourceNodeSummaries, TargetCollection, eval, - eval_par, eval_par_with_skip, eval_with_skip, scratch_len, scratch_len_par, update_summaries, + eval_par, scratch_len, scratch_len_par, update_summaries, }; /// Diagnostic information returned by stateless hierarchical solves. @@ -613,8 +613,10 @@ where 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_with_skip( + (true, Some(skip)) => eval_par( &kernel, source_tree.as_view(), &source_summaries.node_summaries, @@ -622,7 +624,7 @@ where targets, moments, theta, - skip, + Some(skip), out_components, &mut scratch, ), @@ -634,10 +636,11 @@ where targets, moments, theta, + None, out_components, &mut scratch, ), - (false, Some(skip)) => eval_with_skip( + (false, Some(skip)) => eval( &kernel, source_tree.as_view(), &source_summaries.node_summaries, @@ -645,7 +648,7 @@ where targets, moments, theta, - skip, + Some(skip), out_components, &mut scratch, ), @@ -657,6 +660,7 @@ where targets, moments, theta, + None, out_components, &mut scratch, ), diff --git a/src/physics/hierarchical/evaluator.rs b/src/physics/hierarchical/evaluator.rs index 2dfd84c..9acb52d 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -91,85 +91,12 @@ 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. The output slice count must match the kernel output +/// dimension `D`. #[inline] pub fn eval( - kernel: &K, - source_tree: ClusterTreeView<'_, T>, - source_summaries: &[K::SourceSummary], - sources: S, - targets: C, - moments: M, - theta: T, - out: [&mut [T]; D], - scratch: &mut EvaluationScratch<'_, [T; D]>, -) -> HierarchicalError -where - K: HierarchicalKernel, - T: Scalar, - K::TargetGeometry: Copy, - S: SourceCollection, - M: SourceMomentCollection, - C: TargetCollection, -{ - eval_optional_skip( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - None, - out, - scratch, - ) -} - -/// Evaluate vector-valued targets while omitting one interaction class. -/// -/// The source tree is still traversed normally so acceptance decisions do not -/// change. [`Skip::Near`] returns accepted far-summary contributions only; -/// [`Skip::Far`] returns direct leaf contributions only. [`Skip::Both`] zeroes the output without -/// target summarization or source-tree traversal. -#[inline] -pub fn eval_with_skip( - kernel: &K, - source_tree: ClusterTreeView<'_, T>, - source_summaries: &[K::SourceSummary], - sources: S, - targets: C, - moments: M, - theta: T, - skip: Skip, - out: [&mut [T]; D], - scratch: &mut EvaluationScratch<'_, [T; D]>, -) -> HierarchicalError -where - K: HierarchicalKernel, - T: Scalar, - K::TargetGeometry: Copy, - S: SourceCollection, - M: SourceMomentCollection, - C: TargetCollection, -{ - eval_optional_skip( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - Some(skip), - out, - scratch, - ) -} - -#[inline] -fn eval_optional_skip( kernel: &K, source_tree: ClusterTreeView<'_, T>, source_summaries: &[K::SourceSummary], @@ -193,23 +120,66 @@ where if err != HierarchicalError::Ok { return err; } - eval_validated( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - skip, - 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], @@ -217,7 +187,6 @@ fn eval_validated( targets: C, moments: M, theta: T, - skip: Option, out: [&mut [T]; D], scratch: &mut EvaluationScratch<'_, [T; D]>, ) -> HierarchicalError @@ -243,12 +212,6 @@ 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() || scratch.contribution.is_empty() { return HierarchicalError::ScratchTooSmall; } @@ -260,7 +223,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, @@ -268,7 +231,6 @@ where target, moments, theta, - skip, &mut target_out, &mut scratch.contribution[0], &mut target_summary, @@ -345,7 +307,7 @@ fn traverse_source_tree( } /// Terminal-node visitor that evaluates far summaries and direct leaf sources. -struct EvaluationTraversalVisitor<'a, K, S, M> +struct EvaluationTraversalVisitor<'a, K, S, M, const EVALUATE_NEAR: bool, const EVALUATE_FAR: bool> where K: HierarchicalKernel, { @@ -353,13 +315,13 @@ where sources: S, target: K::TargetGeometry, moments: M, - skip: Option, out: &'a mut K::Output, contribution: &'a mut K::Output, target_summary: &'a K::TargetSummary, } -impl TraversalVisitor for EvaluationTraversalVisitor<'_, K, S, M> +impl TraversalVisitor + for EvaluationTraversalVisitor<'_, K, S, M, EVALUATE_NEAR, EVALUATE_FAR> where K: HierarchicalKernel, S: SourceCollection, @@ -372,7 +334,7 @@ where _source_level: u32, source_summary: &K::SourceSummary, ) { - if self.skip != Some(Skip::Far) { + if EVALUATE_FAR { self.kernel .eval_far(self.target_summary, source_summary, self.contribution); self.kernel.accumulate(self.out, self.contribution); @@ -381,7 +343,7 @@ where #[inline] fn on_near_leaf(&mut self, _source_node_index: usize, _source_level: u32, source_ids: &[u32]) { - if self.skip == Some(Skip::Near) { + if !EVALUATE_NEAR { return; } for &source_id in source_ids { @@ -400,7 +362,7 @@ where /// 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], @@ -408,7 +370,6 @@ fn eval_scalar( target: K::TargetGeometry, moments: M, theta: K::Scalar, - skip: Option, out: &mut K::Output, contribution: &mut K::Output, target_summary: &mut K::TargetSummary, @@ -428,12 +389,11 @@ where return err; } - let mut visitor = EvaluationTraversalVisitor { + let mut visitor = EvaluationTraversalVisitor:: { kernel, sources, target, moments, - skip, out, contribution, target_summary, @@ -456,83 +416,11 @@ 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. #[inline] pub fn eval_par( - kernel: &K, - source_tree: ClusterTreeView<'_, T>, - source_summaries: &[K::SourceSummary], - sources: S, - targets: C, - moments: M, - theta: T, - out: [&mut [T]; D], - scratch: &mut EvaluationScratch<'_, [T; D]>, -) -> HierarchicalError -where - K: HierarchicalKernel + Sync, - T: Scalar, - K::TargetGeometry: Copy, - S: SourceCollection, - M: SourceMomentCollection, - C: TargetCollection, -{ - eval_par_optional_skip( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - None, - out, - scratch, - ) -} - -/// Evaluate vector-valued targets in parallel while omitting interaction classes. -/// -/// [`Skip::Both`] zeroes the output without target summarization, source-tree traversal, or -/// parallel scratch use. -#[inline] -pub fn eval_par_with_skip( - kernel: &K, - source_tree: ClusterTreeView<'_, T>, - source_summaries: &[K::SourceSummary], - sources: S, - targets: C, - moments: M, - theta: T, - skip: Skip, - out: [&mut [T]; D], - scratch: &mut EvaluationScratch<'_, [T; D]>, -) -> HierarchicalError -where - K: HierarchicalKernel + Sync, - T: Scalar, - K::TargetGeometry: Copy, - S: SourceCollection, - M: SourceMomentCollection, - C: TargetCollection, -{ - eval_par_optional_skip( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - Some(skip), - out, - scratch, - ) -} - -#[inline] -fn eval_par_optional_skip( kernel: &K, source_tree: ClusterTreeView<'_, T>, source_summaries: &[K::SourceSummary], @@ -570,12 +458,6 @@ 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; } @@ -590,27 +472,76 @@ where } let error_code = AtomicU32::new(HierarchicalError::Ok as u32); - eval_par_chunks( - kernel, - source_tree, - source_summaries, - sources, - targets, - moments, - theta, - skip, - 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) => eval_par_chunks::( + kernel, + source_tree, + source_summaries, + sources, + targets, + moments, + theta, + out, + &mut scratch.contribution[..chunk_count], + chunk_size, + &error_code, + ), + } 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], @@ -618,7 +549,6 @@ fn eval_par_chunks( targets: C, moments: M, theta: T, - skip: Option, out: [&mut [T]; D], scratch_contributions: &mut [[T; D]], chunk_size: usize, @@ -640,7 +570,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, @@ -648,7 +578,6 @@ fn eval_par_chunks( targets, moments, theta, - skip, out, &mut chunk_scratch, ); @@ -673,7 +602,7 @@ fn eval_par_chunks( rayon::join( || { - eval_par_chunks( + eval_par_chunks::( kernel, source_tree, source_summaries, @@ -681,7 +610,6 @@ fn eval_par_chunks( left_targets, moments, theta, - skip, left_out, left_scratch, chunk_size, @@ -689,7 +617,7 @@ fn eval_par_chunks( ); }, || { - eval_par_chunks( + eval_par_chunks::( kernel, source_tree, source_summaries, @@ -697,7 +625,6 @@ fn eval_par_chunks( right_targets, moments, theta, - skip, right_out, right_scratch, chunk_size, diff --git a/src/physics/hierarchical/mod.rs b/src/physics/hierarchical/mod.rs index 1700b2e..0731b14 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -39,8 +39,8 @@ pub(crate) use evaluator::{ EvaluationScratch, SourceNodeSummaries, scratch_len, scratch_len_par, update_summaries, }; pub use evaluator::{ - NearFieldInteractionMap, TraversalDiagnostics, eval, eval_par, eval_par_with_skip, - eval_with_skip, traversal_diagnostics, traversal_diagnostics_par, + NearFieldInteractionMap, TraversalDiagnostics, eval, eval_par, traversal_diagnostics, + traversal_diagnostics_par, }; pub use kernel::Skip; pub(crate) use kernel::{ diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index 4702709..cbb456e 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -302,6 +302,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( targets.as_slice(), &moments, 0.5, + None, [&mut full], &mut scratch, ), @@ -314,7 +315,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( kernel.far_calls.store(0, Ordering::Relaxed); let mut far_only = [0.0]; assert_eq!( - super::eval_with_skip( + super::eval( &kernel, source_tree.as_view(), &summaries.node_summaries, @@ -322,7 +323,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( targets.as_slice(), &moments, 0.5, - Skip::Near, + Some(Skip::Near), [&mut far_only], &mut scratch, ), @@ -335,7 +336,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( kernel.far_calls.store(0, Ordering::Relaxed); let mut near_only = [0.0]; assert_eq!( - super::eval_with_skip( + super::eval( &kernel, source_tree.as_view(), &summaries.node_summaries, @@ -343,7 +344,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( targets.as_slice(), &moments, 0.5, - Skip::Far, + Some(Skip::Far), [&mut near_only], &mut scratch, ), @@ -357,7 +358,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( kernel.far_calls.store(0, Ordering::Relaxed); let mut far_only_par = [0.0]; assert_eq!( - super::eval_par_with_skip( + super::eval_par( &kernel, source_tree.as_view(), &summaries.node_summaries, @@ -365,7 +366,7 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( targets.as_slice(), &moments, 0.5, - Skip::Near, + Some(Skip::Near), [&mut far_only_par], &mut scratch, ), @@ -374,56 +375,6 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( 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); - - kernel.target_summary_calls.store(0, Ordering::Relaxed); - kernel.accept_calls.store(0, Ordering::Relaxed); - kernel.near_calls.store(0, Ordering::Relaxed); - kernel.far_calls.store(0, Ordering::Relaxed); - let mut skipped_both = [f64::NAN]; - let mut empty_contribution = []; - let mut empty_scratch = EvaluationScratch { - contribution: &mut empty_contribution, - }; - assert_eq!( - super::eval_with_skip( - &kernel, - source_tree.as_view(), - &[], - sources.as_slice(), - targets.as_slice(), - &moments, - 0.5, - Skip::Both, - [&mut skipped_both], - &mut empty_scratch, - ), - HierarchicalError::Ok - ); - assert_eq!(skipped_both, [0.0]); - assert_eq!(kernel.target_summary_calls.load(Ordering::Relaxed), 0); - assert_eq!(kernel.accept_calls.load(Ordering::Relaxed), 0); - assert_eq!(kernel.near_calls.load(Ordering::Relaxed), 0); - assert_eq!(kernel.far_calls.load(Ordering::Relaxed), 0); - - skipped_both.fill(f64::NAN); - assert_eq!( - super::eval_par_with_skip( - &kernel, - source_tree.as_view(), - &[], - sources.as_slice(), - targets.as_slice(), - &moments, - 0.5, - Skip::Both, - [&mut skipped_both], - &mut empty_scratch, - ), - HierarchicalError::Ok - ); - assert_eq!(skipped_both, [0.0]); - assert_eq!(kernel.target_summary_calls.load(Ordering::Relaxed), 0); - assert_eq!(kernel.accept_calls.load(Ordering::Relaxed), 0); } #[test] @@ -597,6 +548,7 @@ where targets, moments, theta, + None, column_slices, scratch, ); @@ -644,6 +596,7 @@ where targets, moments, theta, + None, column_slices, scratch, ); From 0e93eff78be38636f5f2aebdd9867c0d548b90ee Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 16:44:05 -0400 Subject: [PATCH 17/20] roll major version --- CHANGELOG.md | 4 ++-- Cargo.toml | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ad1053f..918da1e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,10 +1,10 @@ # Changelog -## 12.2.0 2026-08-13 +## 13.0.0 2026-08-13 * 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 + * !Add option for hierarchical evaluator to skip evaluating kernels for near-field, far-field, or both * Add sparse inductance matrix for linear filaments using CSC interaction map * Python * Plumb in bindings to new `skip` option and sparse inductance matrix diff --git a/Cargo.toml b/Cargo.toml index 87cce64..8fff106 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "12.2.0" +version = "13.0.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" From b823e5160bddfb173db4e2f2a0bc151144513b2c Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 16:45:52 -0400 Subject: [PATCH 18/20] accept scipy csc array --- Cargo.lock | 2 +- cfsem/bindings.py | 10 +++++----- test/test_electromagnetics.py | 7 ++++--- 3 files changed, 10 insertions(+), 9 deletions(-) diff --git a/Cargo.lock b/Cargo.lock index 357d934..1639f71 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -128,7 +128,7 @@ checksum = "9330f8b2ff13f34540b44e946ef35111825727b38d33286ef986142615121801" [[package]] name = "cfsem" -version = "12.2.0" +version = "13.0.0" dependencies = [ "criterion", "faer", diff --git a/cfsem/bindings.py b/cfsem/bindings.py index 01318e6..424a82c 100644 --- a/cfsem/bindings.py +++ b/cfsem/bindings.py @@ -12,7 +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_matrix +from scipy.sparse import csc_array, csc_matrix from cfsem.types import Array3xN @@ -1112,7 +1112,7 @@ def inductance_linear_filaments_sparse( dlxyzfil_tgt: Array3xN, xyzfil_src: Array3xN, dlxyzfil_src: Array3xN, - interaction_map: csc_matrix, + interaction_map: csc_matrix | csc_array, wire_radius_src: float | NDArray[float64] = 0.0, par: bool = True, ) -> csc_matrix: @@ -1136,12 +1136,12 @@ def inductance_linear_filaments_sparse( [H] CSC inductance matrix with the supplied sparsity pattern Raises: - TypeError: If ``interaction_map`` is not a SciPy ``csc_matrix``. + 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): - raise TypeError("interaction_map must be a scipy.sparse.csc_matrix") + 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") diff --git a/test/test_electromagnetics.py b/test/test_electromagnetics.py index 0fbac4e..3db2506 100644 --- a/test/test_electromagnetics.py +++ b/test/test_electromagnetics.py @@ -1055,7 +1055,8 @@ def test_inductance_linear_filaments_matrix_contracts_to_vector(par): @mark.parametrize("par", [True, False]) -def test_inductance_linear_filaments_sparse_preserves_csc_pattern(par): +@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]), @@ -1074,7 +1075,7 @@ def test_inductance_linear_filaments_sparse_preserves_csc_pattern(par): ) 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 = sparse.csc_matrix( + interaction_map = csc_type( (np.full(row_indices.size, np.nan), row_indices, column_pointers), shape=(3, 4), ) @@ -1129,7 +1130,7 @@ 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"): + 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") From 4f390cbc776b8ac0f0068446433e5080b93f16d8 Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 17:39:13 -0400 Subject: [PATCH 19/20] changelog entry date --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 918da1e..4a1bb73 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,6 @@ # Changelog -## 13.0.0 2026-08-13 +## 13.0.0 2026-08-19 * Rust * Add near-field interaction map to hierarchical diagnostics outputs From 754d508c7c88bb864bfe9d9feaba7959cd465240 Mon Sep 17 00:00:00 2001 From: James Logan Date: Wed, 19 Aug 2026 18:07:40 -0400 Subject: [PATCH 20/20] update rust hierarchical convenience functions for parity with python --- CHANGELOG.md | 1 + benches/linear_filament.rs | 4 + benches/point_source.rs | 4 + src/physics/hierarchical/convenience.rs | 111 +++++++++- src/physics/hierarchical/evaluator.rs | 34 +-- src/physics/hierarchical/mod.rs | 2 +- src/physics/hierarchical/tests.rs | 98 +++++++++ src/python.rs | 261 ++++-------------------- 8 files changed, 261 insertions(+), 254 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4a1bb73..265496b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ * 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 diff --git a/benches/linear_filament.rs b/benches/linear_filament.rs index cdee3b0..5a7d505 100644 --- a/benches/linear_filament.rs +++ b/benches/linear_filament.rs @@ -187,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(), @@ -217,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(), @@ -323,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(), @@ -354,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(), 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/src/physics/hierarchical/convenience.rs b/src/physics/hierarchical/convenience.rs index 8aa9709..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, - Skip, 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); @@ -91,6 +123,7 @@ pub fn flux_density_dipole_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -114,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, @@ -130,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); @@ -144,6 +179,7 @@ pub fn vector_potential_dipole_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -168,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, @@ -185,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); @@ -198,6 +236,7 @@ pub fn flux_density_linear_filament_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -222,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, @@ -239,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); @@ -252,6 +293,7 @@ pub fn vector_potential_linear_filament_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -282,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, @@ -297,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) @@ -313,6 +357,7 @@ pub fn flux_density_triangle_mesh_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -343,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, @@ -358,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) @@ -374,6 +421,7 @@ pub fn vector_potential_triangle_mesh_hierarchical( theta, par, None, + extra_diagnostics, out, ) } @@ -388,6 +436,7 @@ pub fn flux_density_dipole_hierarchical_with_skip( 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); @@ -402,6 +451,7 @@ pub fn flux_density_dipole_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -416,6 +466,7 @@ pub fn vector_potential_dipole_hierarchical_with_skip( 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); @@ -430,6 +481,7 @@ pub fn vector_potential_dipole_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -445,6 +497,7 @@ pub fn flux_density_linear_filament_hierarchical_with_skip( 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); @@ -458,6 +511,7 @@ pub fn flux_density_linear_filament_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -473,6 +527,7 @@ pub fn vector_potential_linear_filament_hierarchical_with_skip( 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); @@ -486,6 +541,7 @@ pub fn vector_potential_linear_filament_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -499,6 +555,7 @@ pub fn flux_density_triangle_mesh_hierarchical_with_skip( theta: f64, par: bool, skip: Skip, + extra_diagnostics: bool, out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result>, HierarchicalError> { mesh.validate_nodal_values(s) @@ -515,6 +572,7 @@ pub fn flux_density_triangle_mesh_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -528,6 +586,7 @@ pub fn vector_potential_triangle_mesh_hierarchical_with_skip( theta: f64, par: bool, skip: Skip, + extra_diagnostics: bool, out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result>, HierarchicalError> { mesh.validate_nodal_values(s) @@ -544,6 +603,7 @@ pub fn vector_potential_triangle_mesh_hierarchical_with_skip( theta, par, Some(skip), + extra_diagnostics, out, ) } @@ -558,6 +618,7 @@ pub(crate) fn one_shot_vec3( theta: T, par: bool, skip: Option, + extra_diagnostics: bool, out: (&mut [T], &mut [T], &mut [T]), ) -> Result, HierarchicalError> where @@ -584,7 +645,7 @@ where BuildMethod::LongestAxis => ClusterTree::build(sources)?, BuildMethod::MortonLbvh => ClusterTree::build_morton_lbvh(sources)?, }; - let source_summaries = if skip == Some(Skip::Both) { + let source_summaries = if skip == Some(Skip::Both) && !extra_diagnostics { None } else { let mut source_summaries = SourceNodeSummaries::::new(source_tree.as_view()); @@ -603,7 +664,10 @@ where let construction_seconds = construction_start.elapsed().as_secs_f64(); let evaluation_start = Instant::now(); - if let Some(source_summaries) = source_summaries { + 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(), @@ -675,11 +739,36 @@ where } 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 9acb52d..a20e3c3 100644 --- a/src/physics/hierarchical/evaluator.rs +++ b/src/physics/hierarchical/evaluator.rs @@ -93,8 +93,9 @@ where /// as a single target leaf, walked against the source tree, and written directly /// 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. The output slice count must match the kernel output -/// dimension `D`. +/// 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, @@ -212,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; } @@ -418,7 +425,8 @@ where /// 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. /// [`Skip::Near`] retains accepted far-summary contributions only, [`Skip::Far`] retains direct -/// leaf contributions only, and `None` evaluates both interaction classes. +/// 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, @@ -458,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; } @@ -512,19 +526,7 @@ where chunk_size, &error_code, ), - Some(Skip::Both) => 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)) diff --git a/src/physics/hierarchical/mod.rs b/src/physics/hierarchical/mod.rs index 0731b14..0f286fb 100644 --- a/src/physics/hierarchical/mod.rs +++ b/src/physics/hierarchical/mod.rs @@ -21,7 +21,7 @@ pub mod kernels; pub mod tree; pub use convenience::{ - flux_density_dipole_hierarchical, flux_density_dipole_hierarchical_with_skip, + 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, diff --git a/src/physics/hierarchical/tests.rs b/src/physics/hierarchical/tests.rs index cbb456e..58047a7 100644 --- a/src/physics/hierarchical/tests.rs +++ b/src/physics/hierarchical/tests.rs @@ -377,6 +377,62 @@ fn filtered_evaluation_skips_required_kernel_calls_and_reconstructs_full_output( 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]]); @@ -395,6 +451,7 @@ fn one_shot_skip_both_bypasses_field_traversal_and_zeroes_outputs() { 0.5, true, Some(Skip::Both), + false, (&mut out0, &mut out1, &mut out2), ) .unwrap(); @@ -404,6 +461,47 @@ fn one_shot_skip_both_bypasses_field_traversal_and_zeroes_outputs() { 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] diff --git a/src/python.rs b/src/python.rs index a28aa20..70258eb 100644 --- a/src/python.rs +++ b/src/python.rs @@ -654,56 +654,7 @@ fn parse_hierarchical_skip( } } -/// Compute traversal diagnostics after rebuilding source summaries for the returned tree. -fn traversal_diagnostics_for_python( - kernel: K, - source_tree: &physics::hierarchical::tree::ClusterTree, - sources: S, - targets: C, - moments: M, - theta: f64, - par: bool, -) -> 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 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 diagnostics = match par { - true => physics::hierarchical::traversal_diagnostics_par( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - targets, - theta, - ), - false => physics::hierarchical::traversal_diagnostics( - &kernel, - source_tree.as_view(), - &source_summaries.node_summaries, - targets, - theta, - ), - }; - diagnostics.map_err(|err| py_hierarchical_error("hierarchical traversal diagnostic", err)) -} - -/// Convert an owned CSC interaction pattern into a SciPy sparse matrix. +/// Move a CSC interaction pattern into a SciPy sparse matrix. fn near_field_interaction_map_object( py: Python<'_>, map: physics::hierarchical::NearFieldInteractionMap, @@ -726,48 +677,26 @@ type OptionalDiagnosticsPy = ( 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, - par: bool, -} - -/// 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 { + let Some(traversal) = diagnostics.take_traversal_diagnostics() else { return Ok((None, None, None)); - } - - let diagnostics = traversal_diagnostics_for_python( - request.kernel, - request.source_tree, - request.sources, - request.targets, - request.moments, - request.theta, - request.par, - )?; + }; Ok(( - Some(source_tree_diagnostics_object(py, request.source_tree)?), - Some(PyArray1::from_vec(py, diagnostics.accepted_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, - diagnostics.near_field_interaction_map, + traversal.near_field_interaction_map, )?), )) } @@ -809,7 +738,7 @@ fn flux_density_dipole_hierarchical( let obs = read_xyz_tuple(py, &obs, "obs")?; let construction_method = parse_build_method(construction_method)?; let skip = parse_hierarchical_skip(skip)?; - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.len(), "flux_density", |out| { match skip { Some(skip) => physics::hierarchical::flux_density_dipole_hierarchical_with_skip( @@ -821,6 +750,7 @@ fn flux_density_dipole_hierarchical( theta, par, skip, + extra_diagnostics, out, ), None => physics::hierarchical::flux_density_dipole_hierarchical( @@ -831,41 +761,14 @@ fn flux_density_dipole_hierarchical( 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, near_field_interaction_map) = - optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::DipoleFluxDensityKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -916,7 +819,7 @@ fn vector_potential_dipole_hierarchical( let obs = read_xyz_tuple(py, &obs, "obs")?; let construction_method = parse_build_method(construction_method)?; let skip = parse_hierarchical_skip(skip)?; - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, obs.len(), "vector_potential", |out| { match skip { Some(skip) => { @@ -929,6 +832,7 @@ fn vector_potential_dipole_hierarchical( theta, par, skip, + extra_diagnostics, out, ) } @@ -940,41 +844,14 @@ fn vector_potential_dipole_hierarchical( 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, near_field_interaction_map) = - optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: physics::hierarchical::kernels::DipoleVectorPotentialKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1027,7 +904,7 @@ fn flux_density_linear_filament_hierarchical( let wire_radius = read_float_input_array1(py, &wire_radius, "wire_radius")?; let construction_method = parse_build_method(construction_method)?; let skip = parse_hierarchical_skip(skip)?; - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, xyzp.len(), "flux_density", |out| { match skip { Some(skip) => { @@ -1041,6 +918,7 @@ fn flux_density_linear_filament_hierarchical( theta, par, skip, + extra_diagnostics, out, ) } @@ -1053,36 +931,14 @@ fn flux_density_linear_filament_hierarchical( 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, near_field_interaction_map) = - 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, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1135,7 +991,7 @@ fn vector_potential_linear_filament_hierarchical( let wire_radius = read_float_input_array1(py, &wire_radius, "wire_radius")?; let construction_method = parse_build_method(construction_method)?; let skip = parse_hierarchical_skip(skip)?; - let (field, diagnostics) = + let (field, mut diagnostics) = evaluate_hierarchical_vec3(py, out, xyzp.len(), "vector_potential", |out| { match skip { Some(skip) => { @@ -1149,6 +1005,7 @@ fn vector_potential_linear_filament_hierarchical( theta, par, skip, + extra_diagnostics, out, ) } @@ -1161,6 +1018,7 @@ fn vector_potential_linear_filament_hierarchical( construction_method, theta, par, + extra_diagnostics, out, ), } @@ -1168,32 +1026,8 @@ fn vector_potential_linear_filament_hierarchical( 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, near_field_interaction_map) = - 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, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1237,7 +1071,7 @@ fn flux_density_triangle_mesh_hierarchical( 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(), @@ -1248,25 +1082,13 @@ fn flux_density_triangle_mesh_hierarchical( theta, par, skip, + extra_diagnostics, out, ) .map_err(|err| py_hierarchical_error("hierarchical triangle-mesh flux density", err)) })?; let (source_tree, accepted_levels, near_field_interaction_map) = - optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: - physics::hierarchical::kernels::BoundaryElementFluxDensityKernel::::new(), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field, @@ -1310,7 +1132,7 @@ fn vector_potential_triangle_mesh_hierarchical( 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(), @@ -1321,6 +1143,7 @@ fn vector_potential_triangle_mesh_hierarchical( theta, par, skip, + extra_diagnostics, out, ) .map_err(|err| { @@ -1328,21 +1151,7 @@ fn vector_potential_triangle_mesh_hierarchical( }) })?; let (source_tree, accepted_levels, near_field_interaction_map) = - optional_hierarchical_diagnostics( - py, - extra_diagnostics, - HierarchicalDiagnosticRequest { - kernel: - physics::hierarchical::kernels::BoundaryElementVectorPotentialKernel::::new( - ), - source_tree: &diagnostics.source_tree, - sources, - targets, - moments, - theta, - par, - }, - )?; + hierarchical_diagnostics_objects(py, &mut diagnostics)?; solve_result_from_field( py, field,