Performance improvements (#22)

## 2.6.0 2025-10-15

Substantial performance improvement for `flux_density_linear_filament`
Biot-Savart methods.
This also improves performance in calculations that use these methods,
such as linear filament
body force density calcs.

### Added

* Rust
    * Add `dot3f` and `cross3f` 32-bit float variants

### Changed

* Rust
* Use mixed-precision method for `flux_density_linear_filament_scalar`
* High-dynamic-range part of the calc is still done using 64-bit floats
    * Low-dynamic-range part of the calc is now done using 32-bit floats
        * _All_ addition operations in 32-bit section are done using 
        fused multiply-add operations, usually chained to defer
        roundoff to final operation. As a result, total roundoff error
        accumulated in this section is minimal.
* Return is upcast back to 64-bit float to support precise summation
downstream
    * 1.4-2x speedup without any meaningful loss of precision
* No change to unit test tolerances needed; unlike an all-32-bit
implementation,
this mixed-precision method passes all the same tests as the 64-bit-only
method
* Python
    * Update dep versions
* Use latest rust backend version, which includes 1.4-2x speedup for
flux_density_linear_filament Biot-Savart calcs

Author: jlogan-cfs

CHANGELOG.md

@@ -0,0 +1,32 @@
+# Changelog
+
+See archived changelogs for versions prior to 2.6.0.
+
+## 2.6.0 2025-10-15
+
+Substantial performance improvement for `flux_density_linear_filament` Biot-Savart methods.
+This also improves performance in calculations that use these methods, such as linear filament
+body force density calcs.
+
+### Added
+
+* Rust
+    * Add `dot3f` and `cross3f` 32-bit float variants
+
+### Changed
+
+* Rust
+    * Use mixed-precision method for `flux_density_linear_filament_scalar`
+    * High-dynamic-range part of the calc is still done using 64-bit floats
+    * Low-dynamic-range part of the calc is now done using 32-bit floats
+        * _All_ addition operations in 32-bit section are done using 
+        fused multiply-add operations, usually chained to defer
+        roundoff to final operation. As a result, total roundoff error
+        accumulated in this section is minimal.
+    * Return is upcast back to 64-bit float to support precise summation downstream
+    * 1.4-2x speedup without any meaningful loss of precision
+        * No change to unit test tolerances needed; unlike an all-32-bit implementation,
+        this mixed-precision method passes all the same tests as the 64-bit-only method
+* Python
+    * Update dep versions
+    * Use latest rust backend version, which includes 1.4-2x speedup for flux_density_linear_filament Biot-Savart calcs

CHANGELOG_PYTHON.md CHANGELOG_PYTHON_ARCHIVE.mdCHANGELOG_PYTHON.md renamed to CHANGELOG_PYTHON_ARCHIVE.md

@@ -1,6 +1,6 @@
 [package]
 name = "cfsem"
-version = "2.5.0"
+version = "2.6.0"
 edition = "2024"
 authors = ["Commonwealth Fusion Systems <jlogan@cfs.energy>"]
 license = "MIT"
@@ -15,11 +15,11 @@ name = "cfsem"
 crate-type = ["cdylib", "rlib"]
 
 [dependencies]
-pyo3 = { version="0.25.1", features=["extension-module"], optional=true }
-numpy = { version="0.25.0", optional=true }  # This must match pyo3 version!
+pyo3 = { version="0.26.0", features=["extension-module"], optional=true }
+numpy = { version="0.26.0", optional=true }  # This must match pyo3 version!
 
-nalgebra = "^0.33.2"
-rayon = "^1.10.0"
+nalgebra = "^0.34.1"
+rayon = "^1.11.0"
 libm = "^0.2"
 num-traits = { version = "0.2.19", features = ["libm"] }
 

CHANGELOG_RUST.md CHANGELOG_RUST_ARCHIVE.mdCHANGELOG_RUST.md renamed to CHANGELOG_RUST_ARCHIVE.md

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "cfsem"
-version = "2.5.0"
+version = "2.6.0"
 description = "Quasi-steady electromagnetics including filamentized approximations, Biot-Savart, and Grad-Shafranov."
 authors = [{name = "Commonwealth Fusion Systems", email = "jlogan@cfs.energy"}]
 requires-python = ">=3.9, <3.14"
@@ -15,7 +15,7 @@ classifiers = [
 ]
 dependencies = [
     "numpy >= 2",
-    "interpn >= 0.2.5",
+    "interpn >= 0.6.1",
     "findiff >= 0.12.1",
     "pydantic >= 2",
     "pydantic-numpy >= 6"

Cargo.lock

@@ -1,5 +1,5 @@
 #![doc=include_str!("../README.md")]
-#![allow(non_snake_case)]
+#![allow(clippy::doc_overindented_list_items)]
 #![allow(clippy::needless_range_loop)]
 #![allow(clippy::needless_late_init)]
 #![allow(non_snake_case)]

Cargo.toml

@@ -45,6 +45,9 @@ pub fn ellipk(m: f64) -> f64 {
     let mut ellip: f64 = 0.0;
     let c: f64 = 1.0 - m;
     let logterm = c.powi(-1).ln();
+
+    // NOTE: This loop is unrolled at compile-time automatically,
+    // and the repeated calls to `powi` are de-duplicated by the compiler.
     for i in 0..5 {
         ellip = logterm
             .mul_add(ELLIPK_B[i], ELLIPK_A[i])
@@ -69,6 +72,9 @@ pub fn ellipe(m: f64) -> f64 {
     let mut ellip: f64 = 0.0;
     let c: f64 = 1.0 - m;
     let logterm = c.powi(-1).ln();
+
+    // NOTE: This loop is unrolled at compile-time automatically,
+    // and the repeated calls to `powi` are de-duplicated by the compiler.
     for i in 0..5 {
         ellip = logterm
             .mul_add(ELLIPE_B[i], ELLIPE_A[i])
@@ -99,12 +105,35 @@ pub fn cross3(x0: f64, y0: f64, z0: f64, x1: f64, y1: f64, z1: f64) -> (f64, f64
     (cx, cy, cz)
 }
 
+/// Evaluate the cross products for each axis component
+/// separately using `mul_add` which would not be assumed usable
+/// in a more general implementation.
+/// 32-bit float variant.
+#[inline]
+pub fn cross3f(x0: f32, y0: f32, z0: f32, x1: f32, y1: f32, z1: f32) -> (f32, f32, f32) {
+    let xy = -x1 * y0;
+    let yz = -y1 * z0;
+    let zx = -z1 * x0;
+    let cx = y0.mul_add(z1, yz);
+    let cy = z0.mul_add(x1, zx);
+    let cz = x0.mul_add(y1, xy);
+
+    (cx, cy, cz)
+}
+
 /// Scalar dot product using `mul_add`.
 #[inline]
 pub fn dot3(x0: f64, y0: f64, z0: f64, x1: f64, y1: f64, z1: f64) -> f64 {
     x0.mul_add(x1, y0.mul_add(y1, z0 * z1))
 }
 
+/// Scalar dot product using `mul_add`.
+/// 32-bit float variant.
+#[inline]
+pub fn dot3f(x0: f32, y0: f32, z0: f32, x1: f32, y1: f32, z1: f32) -> f32 {
+    x0.mul_add(x1, y0.mul_add(y1, z0 * z1))
+}
+
 /// Convert a point from cartesian to cylindrical coordinates.
 #[inline]
 pub fn cartesian_to_cylindrical(x: f64, y: f64, z: f64) -> (f64, f64, f64) {

pyproject.toml

@@ -1268,7 +1268,7 @@ mod test {
         let zfil: Vec<f64> = (0..NFIL)
             .map(|i| (i as f64) - (NFIL as f64) / 2.0)
             .collect();
-        let ifil: Vec<f64> = (0..NFIL).map(|i| (i as f64)).collect();
+        let ifil: Vec<f64> = (0..NFIL).map(|i| i as f64).collect();
 
         // Build a scattering of observation locations
         let rprime: Vec<f64> = (0..NOBS).map(|i| 2.0 * (i as f64).sin() + 2.1).collect();