diff --git a/.cargo/config.toml b/.cargo/config.toml index ab3794ad..51372ae0 100644 --- a/.cargo/config.toml +++ b/.cargo/config.toml @@ -1,2 +1,2 @@ -[target.'cfg(any(target_arch = "x86", target_arch = "x86_64", target_arch = "x64"))'] -rustflags = ["-C", "target-feature=+sse3,+ssse3,+sse4.1,+sse4.2,+popcnt,+cmpxchg16b,+avx,+avx2,+fma,+bmi1,+bmi2,+lzcnt,+pclmulqdq,+movbe"] +[target.'cfg(any(target_arch = "x86_64", target_arch = "x64"))'] +rustflags = ["-C", "target-cpu=x86-64-v3"] diff --git a/.github/workflows/release_python.yml b/.github/workflows/release_python.yml index 58bc5f56..cb7b24fa 100644 --- a/.github/workflows/release_python.yml +++ b/.github/workflows/release_python.yml @@ -24,28 +24,112 @@ jobs: platform: - runner: ubuntu-22.04 target: x86_64 + maturin_args: --release --out dist --find-interpreter - runner: ubuntu-22.04 target: x86 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: aarch64 + maturin_args: --release --out dist --find-interpreter - runner: ubuntu-22.04 target: armv7 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: s390x + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: ppc64le + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python steps: - uses: actions/checkout@v4 + - name: Mark workspace as safe for git + run: git config --global --add safe.directory "$GITHUB_WORKSPACE" - uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Build wheels uses: PyO3/maturin-action@v1 + env: + RUST_LOG: cc=debug + CARGO_TERM_VERBOSE: "true" with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter --features=python + args: ${{ matrix.platform.maturin_args }} sccache: 'true' - manylinux: auto + manylinux: 2_28 + before-script-linux: | + git config --global --add safe.directory /github/workspace + git config --global --add safe.directory /home/runner/work/cfsem-py/cfsem-py + + if [ "${{ matrix.platform.target }}" = "aarch64" ]; then + export CC=clang + export CXX=clang++ + + # discover GCC version/triple inside the container + GCC_VER="$(gcc -dumpversion)" + GCC_TRIPLE="$(gcc -dumpmachine)" + + export CXXFLAGS="--gcc-toolchain=/usr -stdlib=libstdc++ \ + -isystem /usr/include/c++/${GCC_VER} \ + -isystem /usr/include/c++/${GCC_VER}/${GCC_TRIPLE} \ + -isystem /usr/include/c++/${GCC_VER}/backward" + fi + + if [ "${{ matrix.platform.target }}" = "x86_64" ] || [ "${{ matrix.platform.target }}" = "aarch64" ]; then + echo "PATH=$PATH"; + ls -la /usr/bin | grep -E 'dnf|microdnf|yum|apt|apk' || true; + + # aarch64 runner is newer + if command -v dnf >/dev/null 2>&1; then + echo "Found dnf package manager"; + dnf -y makecache; + dnf -y update; + dnf -y install --skip-broken clang lld gcc gcc-c++ libstdc++-devel openblas-devel || true; + fi + + if command -v apt-get >/dev/null 2>&1; then + echo "Found apt-get package manager"; + apt-get update; + apt-get install -y clang lld g++ libstdc++-12-dev libopenblas-dev || true; + fi + + if command -v microdnf >/dev/null 2>&1; then + echo "Found microdnf package manager"; + microdnf -y update; + microdnf -y install clang lld gcc gcc-c++ libstdc++-devel openblas-devel || true; + fi + + # manylinux2014 runner is based on centos7 + if command -v yum >/dev/null 2>&1; then + echo "Found yum package manager"; + yum -y makecache; + yum -y update; + yum -y --setopt=skip_missing_names_on_install=True install --skip-broken clang lld gcc gcc-c++ libstdc++-devel openblas-devel || true; + yum -y --setopt=skip_missing_names_on_install=True install --skip-broken centos-release-scl || true; + yum -y --setopt=skip_missing_names_on_install=True install --skip-broken llvm-toolset-7 clang llvm-toolset-7-llvm gcc gcc-c++ libstdc++-devel openblas-devel || true; + fi + + if [ -f /opt/rh/llvm-toolset-7/enable ]; then + source /opt/rh/llvm-toolset-7/enable + fi + + if [ "${{ matrix.platform.target }}" = "aarch64" ]; then + # Clang on manylinux doesn't always discover libstdc++ headers without help. + export CC=clang + export CXX=clang++ + export CXXFLAGS="--gcc-toolchain=/usr -stdlib=libstdc++" + fi + + if ! command -v ld.lld >/dev/null 2>&1; then + echo "warning: ld.lld not found after install; continuing without it" >&2 + fi + else + echo "Skipping linux before-script for target ${{ matrix.platform.target }}" + fi + - name: Install auditwheel + run: python -m pip install auditwheel + - name: Repair wheels + run: bindings/repair_wheel.sh dist - name: Upload wheels uses: actions/upload-artifact@v4 with: @@ -60,14 +144,20 @@ jobs: platform: - runner: ubuntu-22.04 target: x86_64 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: x86 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: aarch64 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python - runner: ubuntu-22.04 target: armv7 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python steps: - uses: actions/checkout@v4 + - name: Mark workspace as safe for git + run: git config --global --add safe.directory "$GITHUB_WORKSPACE" - uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} @@ -75,9 +165,13 @@ jobs: uses: PyO3/maturin-action@v1 with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter --features=python + args: ${{ matrix.platform.maturin_args }} sccache: 'true' manylinux: musllinux_1_2 + - name: Install auditwheel + run: python -m pip install auditwheel + - name: Repair wheels + run: bindings/repair_wheel.sh dist - name: Upload wheels uses: actions/upload-artifact@v4 with: @@ -92,20 +186,55 @@ jobs: platform: - runner: windows-latest target: x64 + maturin_args: --release --out dist --find-interpreter - runner: windows-latest target: x86 + maturin_args: --release --out dist --find-interpreter --no-default-features --features=python steps: - uses: actions/checkout@v4 + - name: Mark workspace as safe for git + run: git config --global --add safe.directory "$GITHUB_WORKSPACE" + - name: Initialize submodules + run: git submodule update --init --recursive --depth 1 - uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} architecture: ${{ matrix.platform.target }} + - name: Install OpenBLAS (windows) + shell: pwsh + run: | + $triplet = if ("${{ matrix.platform.target }}" -eq "x86") { "x86-windows" } else { "x64-windows" } + & "C:\vcpkg\vcpkg.exe" install "openblas:$triplet" + "CMAKE_PREFIX_PATH=C:\vcpkg\installed\$triplet" | Out-File -FilePath $env:GITHUB_ENV -Append - name: Build wheels uses: PyO3/maturin-action@v1 + env: + RUST_LOG: cc=debug + CARGO_TERM_VERBOSE: "true" with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter --features=python + args: ${{ matrix.platform.maturin_args }} sccache: 'true' + - name: Dump Boost b2 logs (on failure) + if: failure() + shell: pwsh + run: | + Get-ChildItem -Recurse -Filter b2.log -Path target | ForEach-Object { + Write-Host "=== $($_.FullName) ===" + Write-Host "=== $($_.FullName) (head 200) ===" + Get-Content $_.FullName -TotalCount 200 + Get-Content $_.FullName -Tail 5000 + Write-Host "=== $($_.FullName) (error scan) ===" + Get-Content $_.FullName | Select-String -CaseSensitive -Pattern 'error','failed','LNK','fatal' -Context 2,2 | Select-Object -Last 200 + } + Get-ChildItem -Recurse -Path target -Include *.err,*.log | ForEach-Object { + Write-Host "=== $($_.FullName) ===" + Get-Content $_.FullName -Tail 200 + } + - name: Install delvewheel + run: python -m pip install delvewheel + - name: Repair wheels + run: bindings/repair_wheel.sh dist - name: Upload wheels uses: actions/upload-artifact@v4 with: @@ -121,8 +250,10 @@ jobs: - runner: macos-15-intel target: x86_64 # macos-15 runner is arm64 + maturin_args: --release --out dist --find-interpreter - runner: macos-15 target: aarch64 + maturin_args: --release --out dist --find-interpreter steps: - uses: actions/checkout@v4 - uses: actions/setup-python@v5 @@ -130,10 +261,17 @@ jobs: python-version: ${{ matrix.python-version }} - name: Build wheels uses: PyO3/maturin-action@v1 + env: + RUST_LOG: cc=debug + CARGO_TERM_VERBOSE: "true" with: target: ${{ matrix.platform.target }} - args: --release --out dist --find-interpreter --features=python + args: ${{ matrix.platform.maturin_args }} sccache: 'true' + - name: Install delocate + run: python -m pip install delocate + - name: Repair wheels + run: bindings/repair_wheel.sh dist - name: Upload wheels uses: actions/upload-artifact@v4 with: @@ -144,6 +282,10 @@ jobs: runs-on: ubuntu-latest steps: - uses: actions/checkout@v4 + - name: Mark workspace as safe for git + run: git config --global --add safe.directory "$GITHUB_WORKSPACE" + - name: Initialize submodules + run: git submodule update --init --recursive --depth 1 - name: Build sdist uses: PyO3/maturin-action@v1 with: diff --git a/.github/workflows/test_python.yml b/.github/workflows/test_python.yml index 015000cd..e8ffdb8f 100644 --- a/.github/workflows/test_python.yml +++ b/.github/workflows/test_python.yml @@ -10,13 +10,12 @@ on: push: branches: [ "main", "develop", "release" ] workflow_dispatch: - tag: "Manual Run" workflow_call: jobs: build: runs-on: ubuntu-22.04 - timeout-minutes: 15 + timeout-minutes: 45 strategy: matrix: python-version: ['3.10', '3.11', '3.12', '3.13'] @@ -24,10 +23,18 @@ jobs: steps: - uses: actions/checkout@v4 + - name: Mark workspace as safe for git + run: git config --global --add safe.directory "$GITHUB_WORKSPACE" + + - name: Initialize submodules + run: git submodule update --init --recursive --depth 1 - name: Install Rust toolchain uses: dtolnay/rust-toolchain@stable + - name: Install Linux build deps + run: sudo apt-get update && sudo apt-get install -y clang lld libstdc++-12-dev libopenblas-dev + - name: Install uv uses: astral-sh/setup-uv@v5 with: @@ -47,7 +54,7 @@ jobs: uv run --locked ruff format ./cfsem --check uv run --locked pyright ./cfsem --pythonversion 3.10 uv run --locked coverage run --source=./cfsem -m pytest ./test/ - uv run --locked coverage report + uv run --locked coverage report -m # Build docs sh build_docs.sh diff --git a/.github/workflows/test_rust.yml b/.github/workflows/test_rust.yml index 03165b17..802ad1d3 100644 --- a/.github/workflows/test_rust.yml +++ b/.github/workflows/test_rust.yml @@ -18,6 +18,8 @@ jobs: steps: - uses: actions/checkout@v4 + - name: Install Linux build deps + run: sudo apt-get update && sudo apt-get install -y clang lld libstdc++-12-dev - name: Install Rust toolchain uses: dtolnay/rust-toolchain@stable - name: Format diff --git a/.gitignore b/.gitignore index dde82289..f806f466 100644 --- a/.gitignore +++ b/.gitignore @@ -78,4 +78,7 @@ target site examples/*.png -docs/python/example_outputs/*.png \ No newline at end of file +docs/python/example_outputs/*.png + +# Vendored code and build artefacts +vendor/armadillo-15.2.3/ \ No newline at end of file diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 00000000..5e799856 --- /dev/null +++ b/.gitmodules @@ -0,0 +1,15 @@ +[submodule "vendor/boost-boost-1.90.0"] + path = vendor/boost-boost-1.90.0 + url = https://github.com/boostorg/boost +[submodule "vendor/jsoncpp"] + path = vendor/jsoncpp + url = https://github.com/open-source-parsers/jsoncpp.git +[submodule "vendor/tclap"] + path = vendor/tclap + url = https://github.com/mirror/tclap.git +[submodule "vendor/rat-common"] + path = vendor/rat-common + url = https://gitlab.com/Project-Rat/rat-common.git +[submodule "vendor/rat-mlfmm"] + path = vendor/rat-mlfmm + url = https://gitlab.com/Project-Rat/rat-mlfmm.git diff --git a/CHANGELOG.md b/CHANGELOG.md index 8eeecfd4..7099a769 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,25 @@ # Changelog +## 4.0.0 2025-01-16 + +### Added + +* Add `mlfmm` module with limited bindings to Project Rat's multipole methods for B-field and A-field + * Exposed as `cfsem::mlfmm::fields_linear_filament_mlfmm` when `rat-mlfmm` feature is enabled + * Add vendored dependencies and thin C FFI to support those bindings + * Available on Windows, linux, and macOS only; disabled on other platforms +* !Remove `flux_density_biot_savart` backwards-compatibility alias +* !Upgrade `flux_density_linear_filament`, `vector_potential_linear_filament`, and `body_force_density_linear_filament` functions + * Now handle finite wire length and finite wire thickness analytically + * New `wire_radius` input + * Removed mixed-precision section from flux density calculation in favor of all f64; this no longer affects overall throughput + * For both linear_filament and point_segment variants + * Old point-source segment formulations moved to `point_source::segment` module and available as `flux_density_point_segment` and `vector_potential_point_segment` functions +* Add `biot_savart.py` and `vector_potential.py` examples with plots comparing linear filament, point segment, and multipole calcs +* Improve parallelism heuristics + * Minimum chunk size of 1024 + * Use half of available parallelism as heuristic for physical cores + ## 3.1.0 2025-12-19 ### Added diff --git a/Cargo.lock b/Cargo.lock index de101506..25e74baf 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -2,6 +2,23 @@ # It is not intended for manual editing. version = 4 +[[package]] +name = "adler2" +version = "2.0.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "320119579fcad9c21884f5c4861d16174d0e06250625266f50fe6898340abefa" + +[[package]] +name = "aes" +version = "0.8.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b169f7a6d4742236a0a00c541b845991d0ac43e546831af1249753ab4c3aa3a0" +dependencies = [ + "cfg-if", + "cipher", + "cpufeatures", +] + [[package]] name = "aho-corasick" version = "1.1.3" @@ -11,6 +28,15 @@ dependencies = [ "memchr", ] +[[package]] +name = "alloca" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e5a7d05ea6aea7e9e64d25b9156ba2fee3fdd659e34e41063cd2fc7cd020d7f4" +dependencies = [ + "cc", +] + [[package]] name = "anes" version = "0.1.6" @@ -38,11 +64,20 @@ version = "1.4.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "ace50bade8e6234aa140d9a2f552bbee1db4d353f69b8217bc503490fc1a9f26" +[[package]] +name = "block-buffer" +version = "0.10.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3078c7629b62d3f0439517fa394996acacc5cbc91c5a20d8c658e77abd503a71" +dependencies = [ + "generic-array", +] + [[package]] name = "bumpalo" -version = "3.18.1" +version = "3.19.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "793db76d6187cd04dff33004d8e6c9cc4e05cd330500379d2394209271b4aeee" +checksum = "5dd9dc738b7a8311c7ade152424974d8115f2cdad61e8dab8dac9f2362298510" [[package]] name = "bytemuck" @@ -50,12 +85,33 @@ version = "1.23.1" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "5c76a5792e44e4abe34d3abf15636779261d45a7450612059293d1d2cfc63422" +[[package]] +name = "bzip2" +version = "0.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f3a53fac24f34a81bc9954b5d6cfce0c21e18ec6959f44f56e8e90e4bb7c346c" +dependencies = [ + "libbz2-rs-sys", +] + [[package]] name = "cast" version = "0.3.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "37b2a672a2cb129a2e41c10b1224bb368f9f37a2b16b612598138befd7b37eb5" +[[package]] +name = "cc" +version = "1.2.54" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6354c81bbfd62d9cfa9cb3c773c2b7b2a3a482d569de977fd0e961f6e7c00583" +dependencies = [ + "find-msvc-tools", + "jobserver", + "libc", + "shlex", +] + [[package]] name = "cfg-if" version = "1.0.1" @@ -64,8 +120,10 @@ checksum = "9555578bc9e57714c812a1f84e4fc5b4d21fcb063490c624de019f7464c91268" [[package]] name = "cfsem" -version = "3.1.0" +version = "4.0.0" dependencies = [ + "cc", + "cmake", "criterion", "itertools 0.14.0", "libm", @@ -74,6 +132,7 @@ dependencies = [ "numpy", "pyo3", "rayon", + "zip", ] [[package]] @@ -103,6 +162,16 @@ dependencies = [ "half", ] +[[package]] +name = "cipher" +version = "0.4.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "773f3b9af64447d2ce9850330c473515014aa235e6a783b02db81ff39e4a3dad" +dependencies = [ + "crypto-common", + "inout", +] + [[package]] name = "clap" version = "4.5.40" @@ -128,12 +197,61 @@ version = "0.7.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "b94f61472cee1439c0b966b47e3aca9ae07e45d070759512cd390ea2bebc6675" +[[package]] +name = "cmake" +version = "0.1.57" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "75443c44cd6b379beb8c5b45d85d0773baf31cce901fe7bb252f4eff3008ef7d" +dependencies = [ + "cc", +] + +[[package]] +name = "constant_time_eq" +version = "0.3.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7c74b8349d32d297c9134b8c88677813a227df8f779daa29bfc29c183fe3dca6" + +[[package]] +name = "cpufeatures" +version = "0.2.17" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "59ed5838eebb26a2bb2e58f6d5b5316989ae9d08bab10e0e6d103e656d1b0280" +dependencies = [ + "libc", +] + +[[package]] +name = "crc" +version = "3.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9710d3b3739c2e349eb44fe848ad0b7c8cb1e42bd87ee49371df2f7acaf3e675" +dependencies = [ + "crc-catalog", +] + +[[package]] +name = "crc-catalog" +version = "2.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "19d374276b40fb8bbdee95aef7c7fa6b5316ec764510eb64b8dd0e2ed0d7e7f5" + +[[package]] +name = "crc32fast" +version = "1.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9481c1c90cbf2ac953f07c8d4a58aa3945c425b7185c9154d67a65e4230da511" +dependencies = [ + "cfg-if", +] + [[package]] name = "criterion" -version = "0.7.0" +version = "0.8.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "e1c047a62b0cc3e145fa84415a3191f628e980b194c2755aa12300a4e6cbd928" +checksum = "4d883447757bb0ee46f233e9dc22eb84d93a9508c9b868687b274fc431d886bf" dependencies = [ + "alloca", "anes", "cast", "ciborium", @@ -142,6 +260,7 @@ dependencies = [ "itertools 0.13.0", "num-traits", "oorandom", + "page_size", "plotters", "rayon", "regex", @@ -153,9 +272,9 @@ dependencies = [ [[package]] name = "criterion-plot" -version = "0.6.0" +version = "0.8.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "9b1bcc0dc7dfae599d84ad0b1a55f80cde8af3725da8313b528da95ef783e338" +checksum = "ed943f81ea2faa8dcecbbfa50164acf95d555afec96a27871663b300e387b2e4" dependencies = [ "cast", "itertools 0.13.0", @@ -192,12 +311,95 @@ version = "0.2.3" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "43da5946c66ffcc7745f48db692ffbb10a83bfe0afd96235c5c2a4fb23994929" +[[package]] +name = "crypto-common" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "78c8292055d1c1df0cce5d180393dc8cce0abec0a7102adb6c7b1eef6016d60a" +dependencies = [ + "generic-array", + "typenum", +] + +[[package]] +name = "deflate64" +version = "0.1.10" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "26bf8fc351c5ed29b5c2f0cbbac1b209b74f60ecd62e675a998df72c49af5204" + +[[package]] +name = "deranged" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ececcb659e7ba858fb4f10388c250a7252eb0a27373f1a72b8748afdd248e587" +dependencies = [ + "powerfmt", +] + +[[package]] +name = "digest" +version = "0.10.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9ed9a281f7bc9b7576e61468ba615a66a5c8cfdff42420a70aa82701a3b1e292" +dependencies = [ + "block-buffer", + "crypto-common", + "subtle", +] + [[package]] name = "either" version = "1.15.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719" +[[package]] +name = "equivalent" +version = "1.0.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "877a4ace8713b0bcf2a4e7eec82529c029f1d0619886d18145fea96c3ffe5c0f" + +[[package]] +name = "find-msvc-tools" +version = "0.1.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8591b0bcc8a98a64310a2fae1bb3e9b8564dd10e381e6e28010fde8e8e8568db" + +[[package]] +name = "flate2" +version = "1.1.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b375d6465b98090a5f25b1c7703f3859783755aa9a80433b36e0379a3ec2f369" +dependencies = [ + "crc32fast", + "miniz_oxide", + "zlib-rs", +] + +[[package]] +name = "generic-array" +version = "0.14.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85649ca51fd72272d7821adaf274ad91c288277713d9c18820d8499a7ff69e9a" +dependencies = [ + "typenum", + "version_check", +] + +[[package]] +name = "getrandom" +version = "0.3.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "899def5c37c4fd7b2664648c28120ecec138e4d395b459e5ca34f9cce2dd77fd" +dependencies = [ + "cfg-if", + "js-sys", + "libc", + "r-efi", + "wasip2", + "wasm-bindgen", +] + [[package]] name = "glam" version = "0.14.0" @@ -304,18 +506,52 @@ dependencies = [ "crunchy", ] +[[package]] +name = "hashbrown" +version = "0.16.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "841d1cc9bed7f9236f321df977030373f4a4163ae1a7dbfe1a51a2c1a51d9100" + [[package]] name = "heck" version = "0.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" +[[package]] +name = "hmac" +version = "0.12.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "6c49c37c09c17a53d937dfbb742eb3a961d65a994e6bcdcf37e7399d0cc8ab5e" +dependencies = [ + "digest", +] + +[[package]] +name = "indexmap" +version = "2.13.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7714e70437a7dc3ac8eb7e6f8df75fd8eb422675fc7678aff7364301092b1017" +dependencies = [ + "equivalent", + "hashbrown", +] + [[package]] name = "indoc" version = "2.0.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "f4c7245a08504955605670dbf141fceab975f15ca21570696aebe9d2e71576bd" +[[package]] +name = "inout" +version = "0.1.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "879f10e63c20629ecabbb64a8010319738c66a5cd0c29b02d63d272b03751d01" +dependencies = [ + "generic-array", +] + [[package]] name = "itertools" version = "0.13.0" @@ -340,6 +576,16 @@ version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "4a5f13b858c8d314ee3e8f639011f7ccefe71f97f96e50151fb991f267928e2c" +[[package]] +name = "jobserver" +version = "0.1.34" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9afb3de4395d6b3e67a780b6de64b51c978ecf11cb9a462c66be7d4ca9039d33" +dependencies = [ + "getrandom", + "libc", +] + [[package]] name = "js-sys" version = "0.3.77" @@ -350,6 +596,12 @@ dependencies = [ "wasm-bindgen", ] +[[package]] +name = "libbz2-rs-sys" +version = "0.2.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2c4a545a15244c7d945065b5d392b2d2d7f21526fba56ce51467b06ed445e8f7" + [[package]] name = "libc" version = "0.2.173" @@ -364,9 +616,19 @@ checksum = "f9fbbcab51052fe104eb5e5d351cf728d30a5be1fe14d9be8a3b097481fb97de" [[package]] name = "log" -version = "0.4.27" +version = "0.4.29" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "13dc2df351e3202783a1fe0d44375f7295ffb4049267b0f3018346dc122a1d94" +checksum = "5e5032e24019045c762d3c0f28f5b6b8bbf38563a65908389bf7978758920897" + +[[package]] +name = "lzma-rust2" +version = "0.15.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1670343e58806300d87950e3401e820b519b9384281bbabfb15e3636689ffd69" +dependencies = [ + "crc", + "sha2", +] [[package]] name = "matrixmultiply" @@ -393,6 +655,16 @@ dependencies = [ "autocfg", ] +[[package]] +name = "miniz_oxide" +version = "0.8.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1fa76a2c86f704bdb222d66965fb3d63269ce38518b83cb0575fca855ebb6316" +dependencies = [ + "adler2", + "simd-adler32", +] + [[package]] name = "nalgebra" version = "0.34.1" @@ -470,6 +742,12 @@ dependencies = [ "num-traits", ] +[[package]] +name = "num-conv" +version = "0.1.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "51d515d32fb182ee37cda2ccdcb92950d6a3c2893aa280e540671c2cd0f3b1d9" + [[package]] name = "num-integer" version = "0.1.46" @@ -502,9 +780,9 @@ dependencies = [ [[package]] name = "numpy" -version = "0.27.0" +version = "0.27.1" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "0fa24ffc88cf9d43f7269d6b6a0d0a00010924a8cc90604a21ef9c433b66998d" +checksum = "7aac2e6a6e4468ffa092ad43c39b81c79196c2bb773b8db4085f695efe3bba17" dependencies = [ "libc", "ndarray", @@ -528,12 +806,38 @@ version = "11.1.5" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "d6790f58c7ff633d8771f42965289203411a5e5c68388703c06e14f24770b41e" +[[package]] +name = "page_size" +version = "0.6.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "30d5b2194ed13191c1999ae0704b7839fb18384fa22e49b57eeaa97d79ce40da" +dependencies = [ + "libc", + "winapi", +] + [[package]] name = "paste" version = "1.0.15" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "57c0d7b74b563b49d38dae00a0c37d4d6de9b432382b2892f0574ddcae73fd0a" +[[package]] +name = "pbkdf2" +version = "0.12.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f8ed6a7761f76e3b9f92dfb0a60a6a6477c61024b775147ff0973a02653abaf2" +dependencies = [ + "digest", + "hmac", +] + +[[package]] +name = "pkg-config" +version = "0.3.32" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7edddbd0b52d732b21ad9a5fab5c704c14cd949e5e9a1ec5929a24fded1b904c" + [[package]] name = "plotters" version = "0.3.7" @@ -577,6 +881,18 @@ dependencies = [ "portable-atomic", ] +[[package]] +name = "powerfmt" +version = "0.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "439ee305def115ba05938db6eb1644ff94165c5ab5e9420d1c1bcedbba909391" + +[[package]] +name = "ppmd-rust" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d558c559f0450f16f2a27a1f017ef38468c1090c9ce63c8e51366232d53717b4" + [[package]] name = "proc-macro2" version = "1.0.95" @@ -588,9 +904,9 @@ dependencies = [ [[package]] name = "pyo3" -version = "0.27.1" +version = "0.27.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "37a6df7eab65fc7bee654a421404947e10a0f7085b6951bf2ea395f4659fb0cf" +checksum = "ab53c047fcd1a1d2a8820fe84f05d6be69e9526be40cb03b73f86b6b03e6d87d" dependencies = [ "indoc", "libc", @@ -605,18 +921,18 @@ dependencies = [ [[package]] name = "pyo3-build-config" -version = "0.27.1" +version = "0.27.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "f77d387774f6f6eec64a004eac0ed525aab7fa1966d94b42f743797b3e395afb" +checksum = "b455933107de8642b4487ed26d912c2d899dec6114884214a0b3bb3be9261ea6" dependencies = [ "target-lexicon", ] [[package]] name = "pyo3-ffi" -version = "0.27.1" +version = "0.27.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "2dd13844a4242793e02df3e2ec093f540d948299a6a77ea9ce7afd8623f542be" +checksum = "1c85c9cbfaddf651b1221594209aed57e9e5cff63c4d11d1feead529b872a089" dependencies = [ "libc", "pyo3-build-config", @@ -624,9 +940,9 @@ dependencies = [ [[package]] name = "pyo3-macros" -version = "0.27.1" +version = "0.27.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "eaf8f9f1108270b90d3676b8679586385430e5c0bb78bb5f043f95499c821a71" +checksum = "0a5b10c9bf9888125d917fb4d2ca2d25c8df94c7ab5a52e13313a07e050a3b02" dependencies = [ "proc-macro2", "pyo3-macros-backend", @@ -636,9 +952,9 @@ dependencies = [ [[package]] name = "pyo3-macros-backend" -version = "0.27.1" +version = "0.27.2" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "70a3b2274450ba5288bc9b8c1b69ff569d1d61189d4bff38f8d22e03d17f932b" +checksum = "03b51720d314836e53327f5871d4c0cfb4fb37cc2c4a11cc71907a86342c40f9" dependencies = [ "heck", "proc-macro2", @@ -656,6 +972,12 @@ dependencies = [ "proc-macro2", ] +[[package]] +name = "r-efi" +version = "5.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "69cdb34c158ceb288df11e18b4bd39de994f6657d83847bdffdbd7f346754b0f" + [[package]] name = "rawpointer" version = "0.2.1" @@ -749,18 +1071,28 @@ dependencies = [ [[package]] name = "serde" -version = "1.0.219" +version = "1.0.228" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5f0e2c6ed6606019b4e29e69dbaba95b11854410e5347d525002456dbbb786b6" +checksum = "9a8e94ea7f378bd32cbbd37198a4a91436180c5bb472411e48b5ec2e2124ae9e" +dependencies = [ + "serde_core", + "serde_derive", +] + +[[package]] +name = "serde_core" +version = "1.0.228" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "41d385c7d4ca58e59fc732af25c3983b67ac852c1a25000afe1175de458b67ad" dependencies = [ "serde_derive", ] [[package]] name = "serde_derive" -version = "1.0.219" +version = "1.0.228" source = "registry+https://github.com/rust-lang/crates.io-index" -checksum = "5b0276cf7f2c73365f7157c8123c21cd9a50fbbd844757af28ca1f5925fc2a00" +checksum = "d540f220d3187173da220f885ab66608367b6574e925011a9353e4badda91d79" dependencies = [ "proc-macro2", "quote", @@ -779,6 +1111,34 @@ dependencies = [ "serde", ] +[[package]] +name = "sha1" +version = "0.10.6" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e3bf829a2d51ab4a5ddf1352d8470c140cadc8301b2ae1789db023f01cedd6ba" +dependencies = [ + "cfg-if", + "cpufeatures", + "digest", +] + +[[package]] +name = "sha2" +version = "0.10.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "a7507d819769d01a365ab707794a4084392c824f54a7a6a7862f8c3d0892b283" +dependencies = [ + "cfg-if", + "cpufeatures", + "digest", +] + +[[package]] +name = "shlex" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0fda2ff0d084019ba4d7c6f371c95d8fd75ce3524c3cb8fb653a3023f6323e64" + [[package]] name = "simba" version = "0.9.0" @@ -792,6 +1152,18 @@ dependencies = [ "wide", ] +[[package]] +name = "simd-adler32" +version = "0.3.8" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e320a6c5ad31d271ad523dcf3ad13e2767ad8b1cb8f047f75a8aeaf8da139da2" + +[[package]] +name = "subtle" +version = "2.6.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "13c2bddecc57b384dee18652358fb23172facb8a2c51ccc10d74c157bdea3292" + [[package]] name = "syn" version = "2.0.103" @@ -809,6 +1181,26 @@ version = "0.13.2" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "e502f78cdbb8ba4718f566c418c52bc729126ffd16baee5baa718cf25dd5a69a" +[[package]] +name = "time" +version = "0.3.45" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f9e442fc33d7fdb45aa9bfeb312c095964abdf596f7567261062b2a7107aaabd" +dependencies = [ + "deranged", + "js-sys", + "num-conv", + "powerfmt", + "serde_core", + "time-core", +] + +[[package]] +name = "time-core" +version = "0.1.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8b36ee98fd31ec7426d599183e8fe26932a8dc1fb76ddb6214d05493377d34ca" + [[package]] name = "tinytemplate" version = "1.2.1" @@ -819,6 +1211,12 @@ dependencies = [ "serde_json", ] +[[package]] +name = "typed-path" +version = "0.12.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "7922f2cdc51280d47b491af9eafc41eb0cdab85eabcb390c854412fcbf26dbe8" + [[package]] name = "typenum" version = "1.18.0" @@ -837,6 +1235,12 @@ version = "0.2.4" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "7264e107f553ccae879d21fbea1d6724ac785e8c3bfc762137959b5802826ef3" +[[package]] +name = "version_check" +version = "0.9.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0b928f33d975fc6ad9f86c8f283853ad26bdd5b10b7f1542aa2fa15e2289105a" + [[package]] name = "walkdir" version = "2.5.0" @@ -847,6 +1251,15 @@ dependencies = [ "winapi-util", ] +[[package]] +name = "wasip2" +version = "1.0.2+wasi-0.2.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9517f9239f02c069db75e65f174b3da828fe5f5b945c4dd26bd25d89c03ebcf5" +dependencies = [ + "wit-bindgen", +] + [[package]] name = "wasm-bindgen" version = "0.2.100" @@ -925,6 +1338,22 @@ dependencies = [ "safe_arch", ] +[[package]] +name = "winapi" +version = "0.3.9" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5c839a674fcd7a98952e593242ea400abe93992746761e38641405d28b00f419" +dependencies = [ + "winapi-i686-pc-windows-gnu", + "winapi-x86_64-pc-windows-gnu", +] + +[[package]] +name = "winapi-i686-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ac3b87c63620426dd9b991e5ce0329eff545bccbbb34f3be09ff6fb6ab51b7b6" + [[package]] name = "winapi-util" version = "0.1.9" @@ -934,6 +1363,12 @@ dependencies = [ "windows-sys", ] +[[package]] +name = "winapi-x86_64-pc-windows-gnu" +version = "0.4.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "712e227841d057c1ee1cd2fb22fa7e5a5461ae8e48fa2ca79ec42cfc1931183f" + [[package]] name = "windows-sys" version = "0.59.0" @@ -1006,3 +1441,103 @@ name = "windows_x86_64_msvc" version = "0.52.6" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "589f6da84c646204747d1270a2a5661ea66ed1cced2631d546fdfb155959f9ec" + +[[package]] +name = "wit-bindgen" +version = "0.51.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d7249219f66ced02969388cf2bb044a09756a083d0fab1e566056b04d9fbcaa5" + +[[package]] +name = "zeroize" +version = "1.8.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b97154e67e32c85465826e8bcc1c59429aaaf107c1e4a9e53c8d8ccd5eff88d0" +dependencies = [ + "zeroize_derive", +] + +[[package]] +name = "zeroize_derive" +version = "1.4.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "85a5b4158499876c763cb03bc4e49185d3cccbabb15b33c627f7884f43db852e" +dependencies = [ + "proc-macro2", + "quote", + "syn", +] + +[[package]] +name = "zip" +version = "7.2.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c42e33efc22a0650c311c2ef19115ce232583abbe80850bc8b66509ebef02de0" +dependencies = [ + "aes", + "bzip2", + "constant_time_eq", + "crc32fast", + "deflate64", + "flate2", + "generic-array", + "getrandom", + "hmac", + "indexmap", + "lzma-rust2", + "memchr", + "pbkdf2", + "ppmd-rust", + "sha1", + "time", + "typed-path", + "zeroize", + "zopfli", + "zstd", +] + +[[package]] +name = "zlib-rs" +version = "0.5.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "40990edd51aae2c2b6907af74ffb635029d5788228222c4bb811e9351c0caad3" + +[[package]] +name = "zopfli" +version = "0.8.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "f05cd8797d63865425ff89b5c4a48804f35ba0ce8d125800027ad6017d2b5249" +dependencies = [ + "bumpalo", + "crc32fast", + "log", + "simd-adler32", +] + +[[package]] +name = "zstd" +version = "0.13.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "e91ee311a569c327171651566e07972200e76fcfe2242a4fa446149a3881c08a" +dependencies = [ + "zstd-safe", +] + +[[package]] +name = "zstd-safe" +version = "7.2.4" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "8f49c4d5f0abb602a93fb8736af2a4f4dd9512e36f7f570d66e65ff867ed3b9d" +dependencies = [ + "zstd-sys", +] + +[[package]] +name = "zstd-sys" +version = "2.0.16+zstd.1.5.7" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "91e19ebc2adc8f83e43039e79776e3fda8ca919132d68a1fed6a5faca2683748" +dependencies = [ + "cc", + "pkg-config", +] diff --git a/Cargo.toml b/Cargo.toml index 9c9de289..354f84a1 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "cfsem" -version = "3.1.0" +version = "4.0.0" edition = "2024" authors = ["Commonwealth Fusion Systems "] license = "MIT" @@ -8,6 +8,7 @@ repository = "https://github.com/cfs-energy/cfsem-py" homepage = "https://github.com/cfs-energy/cfsem-py" description = "Quasi-steady electromagnetics including filamentized approximations, Biot-Savart, and Grad-Shafranov." readme = "README.md" +build = "build.rs" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html [lib] @@ -15,8 +16,8 @@ name = "cfsem" crate-type = ["cdylib", "rlib"] [dependencies] -pyo3 = { version="0.27.1", features=["extension-module"], optional=true } -numpy = { version="0.27.0", optional=true } # This must match pyo3 version! +pyo3 = { version="0.27.2", features=["extension-module"], optional=true } +numpy = { version="0.27.1", optional=true } # This must match pyo3 version! nalgebra = "^0.34.1" rayon = "^1.11.0" @@ -24,9 +25,14 @@ libm = "^0.2" num-traits = { version = "0.2.19", features = ["libm"] } [dev-dependencies] -criterion = "0.7.0" +criterion = "0.8.1" itertools = "0.14.0" +[build-dependencies] +cc = "1.2.54" +cmake = "0.1.57" +zip = "7.1.0" + [profile.release] opt-level = 3 lto = true @@ -43,6 +49,7 @@ codegen-units = 16 [features] default = [] python = ["pyo3", "numpy"] +rat-mlfmm = [] [[bench]] name = "linear_filament" diff --git a/README.md b/README.md index 83b6208b..dc81c304 100644 --- a/README.md +++ b/README.md @@ -4,15 +4,33 @@ Quasi-steady electromagnetics including filamentized approximations, Biot-Savart, and Grad-Shafranov. +## Installation - System Dependencies + +When building the `mlfmm` feature, some system dependencies are needed. + +* cmake +* C/C++ toolchain (clang + make) + * Linux: clang, lld, libstdc++-12-dev + * NOTE: the g++ toolchain does not work for this project due to dependency ordering issues during linking + when building cdylib builds for the python extension library. +* git (for submodules) +* BLAS/LAPACK + * Linux & Windows: openblas + * Mac: already included by the OS (Accelerate) +* zlib (for Boost iostreams; typically provided by the OS) + +as well as some run-time dependencies: + +* zlib +* System C++ runtime +* BLAS/LAPACK (MacOS only, nominally provided by OS) + ## Installation - Python Requirements * Python 3.9-3.13 and pip -* Don't worry about this: - * This info provided for troubleshooting purposes: - * If on an x86 processor, you will need a CPU that supports SSE through 4.1, AVX, and FMA. - * This should be true on any modern machine. +* If on an x86 processor, you will need a CPU from roughly 2013 or later. ```bash pip install cfsem @@ -26,7 +44,11 @@ To include this library in a Rust project, add an entry to your Cargo.toml's `[d cfsem = "*" ``` -For Python installation, see the docs for the Python library. +If building with the `rat-mlfmm` feature, the library must be included as a git dependency + +```toml +cfsem = { git = "https://github.com/cfs-energy/cfsem-py.git", tag = "4.0.0" } +``` ## Benchmarking - Rust @@ -36,6 +58,12 @@ Benchmarks are configured in Cargo.toml, and can be run via cargo: cargo bench ``` +To build the docs with katex math rendering: + +```bash +RUSTDOCFLAGS="--html-in-header=katex-header.html" cargo rustdoc --open +``` + ## Development - Python Requirements @@ -48,12 +76,6 @@ To install in the active python environment, do uv pip install -e . --group dev ``` -To build the Rust bindings only, do - -```bash -maturin develop --release --features=python -``` - No part of installation requires root. If access issues are encountered, this can likely be resolved by using a virtual environment. Some computationally-expensive calculations are written in Rust. These calculations and their python bindings are installed from pre-built binaries when installing from pypi or compiled during local development installation, with no intervention from the user in either case. Symmetric bindings with docstrings are available in the `bindings.py` module and re-exported at the library level. diff --git a/benches/linear_filament.rs b/benches/linear_filament.rs index ddbadf5f..2efc4874 100644 --- a/benches/linear_filament.rs +++ b/benches/linear_filament.rs @@ -46,12 +46,14 @@ fn bench_flux_density_linear_filament(c: &mut Criterion) { b.iter(|| { let n = xobs.len(); let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let wire_radius = vec![0.0; ifil.len()]; black_box( flux_density_linear_filament( (&xobs[..], &yobs[..], &zobs[..]), (&xfil[..], &yfil[..], &zfil[..]), (&dlxfil[..], &dlyfil[..], &dlzfil[..]), &ifil[..], + &wire_radius, (&mut bx, &mut by, &mut bz), ) .unwrap(), @@ -72,12 +74,14 @@ fn bench_flux_density_linear_filament(c: &mut Criterion) { b.iter(|| { let n = xobs.len(); let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let wire_radius = vec![0.0; ifil.len()]; black_box( flux_density_linear_filament_par( (&xobs[..], &yobs[..], &zobs[..]), (&xfil[..], &yfil[..], &zfil[..]), (&dlxfil[..], &dlyfil[..], &dlzfil[..]), &ifil[..], + &wire_radius, (&mut bx, &mut by, &mut bz), ) .unwrap(), @@ -131,12 +135,14 @@ fn bench_vector_potential_linear_filament(c: &mut Criterion) { b.iter(|| { let n = xobs.len(); let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let wire_radius = vec![0.0; ifil.len()]; black_box( vector_potential_linear_filament( (&xobs[..], &yobs[..], &zobs[..]), (&xfil[..], &yfil[..], &zfil[..]), (&dlxfil[..], &dlyfil[..], &dlzfil[..]), &ifil[..], + &wire_radius, (&mut bx, &mut by, &mut bz), ) .unwrap(), @@ -157,12 +163,14 @@ fn bench_vector_potential_linear_filament(c: &mut Criterion) { b.iter(|| { let n = xobs.len(); let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let wire_radius = vec![0.0; ifil.len()]; black_box( vector_potential_linear_filament_par( (&xobs[..], &yobs[..], &zobs[..]), (&xfil[..], &yfil[..], &zfil[..]), (&dlxfil[..], &dlyfil[..], &dlzfil[..]), &ifil[..], + &wire_radius, (&mut bx, &mut by, &mut bz), ) .unwrap(), diff --git a/bindings/rat-mlfmm-c/CMakeLists.txt b/bindings/rat-mlfmm-c/CMakeLists.txt new file mode 100644 index 00000000..dfd77ce1 --- /dev/null +++ b/bindings/rat-mlfmm-c/CMakeLists.txt @@ -0,0 +1,280 @@ +cmake_minimum_required(VERSION 3.20) +project(RatMlfmmC LANGUAGES CXX) + +list(INSERT CMAKE_MODULE_PATH 0 ${CMAKE_CURRENT_SOURCE_DIR}/cmake) + +if(CFSEM_EXPECT_RELEASE AND NOT CMAKE_BUILD_TYPE STREQUAL "Release") + message(FATAL_ERROR "CFSEM_EXPECT_RELEASE is ON but CMAKE_BUILD_TYPE is '${CMAKE_BUILD_TYPE}'.") +endif() + +set(RAT_MLFMM_DIR "" CACHE PATH "Path to rat-mlfmm source directory") +if(NOT RAT_MLFMM_DIR) + message(FATAL_ERROR "RAT_MLFMM_DIR is required") +endif() +if(NOT EXISTS "${RAT_MLFMM_DIR}/CMakeLists.txt") + message(FATAL_ERROR "RAT_MLFMM_DIR does not point to a rat-mlfmm source tree") +endif() + +set(RAT_COMMON_DIR "" CACHE PATH "Path to rat-common source directory") +if(NOT RAT_COMMON_DIR) + message(FATAL_ERROR "RAT_COMMON_DIR is required") +endif() +if(NOT EXISTS "${RAT_COMMON_DIR}/CMakeLists.txt") + message(FATAL_ERROR "RAT_COMMON_DIR does not point to a rat-common source tree") +endif() + +set(JSONCPP_SRC_DIR "" CACHE PATH "Path to jsoncpp source directory") +if(NOT JSONCPP_SRC_DIR) + message(FATAL_ERROR "JSONCPP_SRC_DIR is required") +endif() +if(NOT EXISTS "${JSONCPP_SRC_DIR}/CMakeLists.txt") + message(FATAL_ERROR "JSONCPP_SRC_DIR does not point to a jsoncpp source tree") +endif() + +set(ARMADILLO_SRC_DIR "" CACHE PATH "Path to armadillo source directory") +if(NOT ARMADILLO_SRC_DIR) + message(FATAL_ERROR "ARMADILLO_SRC_DIR is required") +endif() +if(NOT EXISTS "${ARMADILLO_SRC_DIR}/CMakeLists.txt") + message(FATAL_ERROR "ARMADILLO_SRC_DIR does not point to an armadillo source tree") +endif() + +set(TCLAP_SRC_DIR "" CACHE PATH "Path to tclap source directory") +if(NOT TCLAP_SRC_DIR) + message(FATAL_ERROR "TCLAP_SRC_DIR is required") +endif() +if(NOT EXISTS "${TCLAP_SRC_DIR}/CMakeLists.txt") + message(FATAL_ERROR "TCLAP_SRC_DIR does not point to a tclap source tree") +endif() + +set(BOOST_SRC_DIR "" CACHE PATH "Path to boost source directory") +if(NOT BOOST_SRC_DIR) + message(FATAL_ERROR "BOOST_SRC_DIR is required") +endif() +if(NOT EXISTS "${BOOST_SRC_DIR}/CMakeLists.txt") + message(FATAL_ERROR "BOOST_SRC_DIR does not point to a boost source tree") +endif() + +set(JSONCPP_WITH_TESTS OFF CACHE BOOL "" FORCE) +set(JSONCPP_WITH_POST_BUILD_UNITTEST OFF CACHE BOOL "" FORCE) +set(JSONCPP_WITH_PKGCONFIG_SUPPORT OFF CACHE BOOL "" FORCE) +set(JSONCPP_WITH_WARNING_AS_ERROR OFF CACHE BOOL "" FORCE) +set(BUILD_SHARED_LIBS OFF CACHE BOOL "" FORCE) +set(BUILD_STATIC_LIBS ON CACHE BOOL "" FORCE) +add_subdirectory(${JSONCPP_SRC_DIR} ${CMAKE_BINARY_DIR}/jsoncpp-build EXCLUDE_FROM_ALL) +set(JsonCPP_INCLUDE_DIR ${JSONCPP_SRC_DIR}/include CACHE PATH "" FORCE) +set(JSONCPP_OUTPUT_DIR ${CMAKE_BINARY_DIR}/lib) +set(JsonCPP_LIBRARY ${JSONCPP_OUTPUT_DIR}/${CMAKE_STATIC_LIBRARY_PREFIX}jsoncpp${CMAKE_STATIC_LIBRARY_SUFFIX} CACHE FILEPATH "" FORCE) +add_custom_command( + OUTPUT ${JsonCPP_LIBRARY} + COMMAND ${CMAKE_COMMAND} -E make_directory ${JSONCPP_OUTPUT_DIR} + COMMAND ${CMAKE_COMMAND} -E copy_if_different $ ${JsonCPP_LIBRARY} + DEPENDS jsoncpp_static +) +add_custom_target(jsoncpp_copy ALL DEPENDS ${JsonCPP_LIBRARY}) + +set(BUILD_SMOKE_TEST OFF CACHE BOOL "" FORCE) +set(ALLOW_BLAS_LAPACK_MACOS ON CACHE BOOL "" FORCE) +set(STATIC_LIB ON CACHE BOOL "" FORCE) +add_subdirectory(${ARMADILLO_SRC_DIR} ${CMAKE_BINARY_DIR}/armadillo-build EXCLUDE_FROM_ALL) +set(Armadillo_DIR ${CMAKE_CURRENT_SOURCE_DIR}/cmake CACHE PATH "" FORCE) +file(WRITE ${CMAKE_BINARY_DIR}/armadillo-build/ArmadilloLibraryDepends.cmake + "# Generated by rat-mlfmm wrapper; Armadillo target provided by add_subdirectory.\\n" +) + +set(TCLAP_INCLUDE_DIR ${TCLAP_SRC_DIR}/include CACHE PATH "" FORCE) + +set(BOOST_INSTALL_DIR ${CMAKE_BINARY_DIR}/boost-install) +set(BOOST_BUILD_DIR ${CMAKE_BINARY_DIR}/boost-build) +set(BOOST_LIB_DIR ${BOOST_INSTALL_DIR}/lib) + +set(_boost_lib_suffix ${CMAKE_STATIC_LIBRARY_SUFFIX}) +set(_boost_lib_prefix ${CMAKE_STATIC_LIBRARY_PREFIX}) + +set(_boost_filesystem ${BOOST_LIB_DIR}/${_boost_lib_prefix}boost_filesystem${_boost_lib_suffix}) +set(_boost_iostreams ${BOOST_LIB_DIR}/${_boost_lib_prefix}boost_iostreams${_boost_lib_suffix}) +set(_boost_thread ${BOOST_LIB_DIR}/${_boost_lib_prefix}boost_thread${_boost_lib_suffix}) +set(_boost_chrono ${BOOST_LIB_DIR}/${_boost_lib_prefix}boost_chrono${_boost_lib_suffix}) + +if(NOT EXISTS "${_boost_filesystem}" OR NOT EXISTS "${_boost_iostreams}" OR NOT EXISTS "${_boost_thread}" OR NOT EXISTS "${_boost_chrono}") + if(CMAKE_BUILD_TYPE STREQUAL "Debug") + set(_boost_variant "debug") + else() + set(_boost_variant "release") + endif() + + if(WIN32) + set(_boost_bootstrap ${BOOST_SRC_DIR}/bootstrap.bat) + set(_boost_b2 ${BOOST_SRC_DIR}/b2.exe) + else() + set(_boost_bootstrap ${BOOST_SRC_DIR}/bootstrap.sh) + set(_boost_b2 ${BOOST_SRC_DIR}/b2) + endif() + + if(NOT EXISTS "${_boost_b2}") + execute_process( + COMMAND ${_boost_bootstrap} + WORKING_DIRECTORY ${BOOST_SRC_DIR} + RESULT_VARIABLE _boost_bootstrap_result + OUTPUT_VARIABLE _boost_bootstrap_stdout + ERROR_VARIABLE _boost_bootstrap_stderr + ) + if(NOT _boost_bootstrap_result EQUAL 0) + message(STATUS "Boost bootstrap stdout:\n${_boost_bootstrap_stdout}") + message(STATUS "Boost bootstrap stderr:\n${_boost_bootstrap_stderr}") + message(FATAL_ERROR "Boost bootstrap failed with code ${_boost_bootstrap_result}") + endif() + endif() + + if(DEFINED CFSEM_BOOST_CXXFLAGS AND NOT CFSEM_BOOST_CXXFLAGS STREQUAL "") + set(_boost_cxxflags "${CFSEM_BOOST_CXXFLAGS}") + else() + set(_boost_cxxflags "-O3") + endif() + + set(_boost_b2_args + --build-dir=${BOOST_BUILD_DIR} + --with-filesystem + --with-iostreams + --with-thread + --with-chrono + -sNO_ICU=1 + -sNO_LZMA=1 + -sNO_ZSTD=1 + -sNO_ZLIB=1 + -sNO_BZIP2=1 + define=BOOST_REGEX_NO_ICU + link=static + runtime-link=static + variant=${_boost_variant} + optimization=speed + cxxflags=${_boost_cxxflags} + --prefix=${BOOST_INSTALL_DIR} + install + ) + if(WIN32) + if(CMAKE_SIZEOF_VOID_P EQUAL 8) + set(_boost_address_model 64) + else() + set(_boost_address_model 32) + endif() + list(APPEND _boost_b2_args toolset=msvc architecture=x86) + list(APPEND _boost_b2_args address-model=${_boost_address_model}) + # Force Windows target to skip POSIX-only config probes (pthread/fallocate). + list(APPEND _boost_b2_args target-os=windows) + list(APPEND _boost_b2_args -sNO_PTHREAD=1) + list(APPEND _boost_b2_args define=BOOST_FILESYSTEM_NO_FALLOCATE) + endif() + message(STATUS "Boost b2 path: ${_boost_b2}") + message(STATUS "Boost b2 args: ${_boost_b2_args}") + execute_process( + COMMAND ${_boost_b2} ${_boost_b2_args} + WORKING_DIRECTORY ${BOOST_SRC_DIR} + RESULT_VARIABLE _boost_build_result + OUTPUT_VARIABLE _boost_build_stdout + ERROR_VARIABLE _boost_build_stderr + ) + if(NOT _boost_build_result EQUAL 0) + set(_boost_build_log "${BOOST_BUILD_DIR}/b2.log") + file(WRITE "${_boost_build_log}" "${_boost_build_stdout}\n${_boost_build_stderr}") + message(STATUS "Boost build log written to: ${_boost_build_log}") + string(LENGTH "${_boost_build_stdout}" _boost_stdout_len) + if(_boost_stdout_len GREATER 4000) + math(EXPR _boost_stdout_start "${_boost_stdout_len} - 4000") + string(SUBSTRING "${_boost_build_stdout}" ${_boost_stdout_start} 4000 _boost_stdout_tail) + else() + set(_boost_stdout_tail "${_boost_build_stdout}") + endif() + string(LENGTH "${_boost_build_stderr}" _boost_stderr_len) + if(_boost_stderr_len GREATER 4000) + math(EXPR _boost_stderr_start "${_boost_stderr_len} - 4000") + string(SUBSTRING "${_boost_build_stderr}" ${_boost_stderr_start} 4000 _boost_stderr_tail) + else() + set(_boost_stderr_tail "${_boost_build_stderr}") + endif() + message(STATUS "Boost build stdout (tail):\n${_boost_stdout_tail}") + message(STATUS "Boost build stderr (tail):\n${_boost_stderr_tail}") + if(WIN32 AND EXISTS "${_boost_filesystem}" AND EXISTS "${_boost_iostreams}" AND EXISTS "${_boost_thread}" AND EXISTS "${_boost_chrono}") + message(WARNING "Boost build reported failure (${_boost_build_result}) but required libs are present; continuing.") + else() + message(FATAL_ERROR "Boost build failed with code ${_boost_build_result}") + endif() + endif() +endif() + +configure_file( + ${CMAKE_CURRENT_SOURCE_DIR}/cmake/BoostConfig.cmake.in + ${CMAKE_CURRENT_BINARY_DIR}/BoostConfig.cmake + @ONLY +) +include(CMakePackageConfigHelpers) +write_basic_package_version_file( + ${CMAKE_CURRENT_BINARY_DIR}/BoostConfigVersion.cmake + VERSION 1.90.0 + COMPATIBILITY AnyNewerVersion +) +set(Boost_DIR ${CMAKE_CURRENT_BINARY_DIR} CACHE PATH "" FORCE) + +set(ENABLE_DOUBLE_PRECISION ON CACHE BOOL "" FORCE) +set(ENABLE_TESTING OFF CACHE BOOL "" FORCE) +set(ENABLE_EXAMPLES OFF CACHE BOOL "" FORCE) +set(ENABLE_BLAS ON CACHE BOOL "" FORCE) +set(ENABLE_SUPERLU OFF CACHE BOOL "" FORCE) +set(ENABLE_MKL_ALLOC OFF CACHE BOOL "" FORCE) +set(RAT_COMMON_BUILD_SHARED OFF CACHE BOOL "" FORCE) + +add_subdirectory(${RAT_COMMON_DIR} ${CMAKE_BINARY_DIR}/rat-common-build EXCLUDE_FROM_ALL) +set(RatCommon_DIR ${CMAKE_BINARY_DIR}/rat-common-build CACHE PATH "" FORCE) + +set(ENABLE_CUDA OFF CACHE BOOL "" FORCE) +set(ENABLE_TESTING OFF CACHE BOOL "" FORCE) +set(ENABLE_EXAMPLES OFF CACHE BOOL "" FORCE) +set(ENABLE_MATLAB OFF CACHE BOOL "" FORCE) +set(RAT_MLFMM_BUILD_SHARED OFF CACHE BOOL "" FORCE) + +add_subdirectory(${RAT_MLFMM_DIR} ${CMAKE_BINARY_DIR}/rat-mlfmm-build EXCLUDE_FROM_ALL) + +if(TARGET ratcmn AND TARGET jsoncpp_lib) + add_dependencies(ratcmn jsoncpp_lib) +endif() +if(TARGET ratcmn AND TARGET jsoncpp_static) + add_dependencies(ratcmn jsoncpp_static) +endif() +if(TARGET ratcmn AND TARGET jsoncpp_copy) + add_dependencies(ratcmn jsoncpp_copy) +endif() + +set(RAT_COMMON_BUILD_INCLUDE ${CMAKE_BINARY_DIR}/rat-common-build/include/rat/common) +file(MAKE_DIRECTORY ${RAT_COMMON_BUILD_INCLUDE}) +file(GLOB RAT_COMMON_HEADERS ${RAT_COMMON_DIR}/include/*.hh) +foreach(_hdr IN LISTS RAT_COMMON_HEADERS) + file(COPY ${_hdr} DESTINATION ${RAT_COMMON_BUILD_INCLUDE}) +endforeach() + +if(TARGET ratmlfmm) + target_include_directories(ratmlfmm PRIVATE ${CMAKE_BINARY_DIR}/rat-common-build/include) +endif() + +add_library(rat_mlfmm_c STATIC + src/rat_mlfmm_c.cpp +) + +target_compile_features(rat_mlfmm_c PRIVATE cxx_std_14) + +target_include_directories(rat_mlfmm_c + PUBLIC + ${CMAKE_CURRENT_SOURCE_DIR}/include + PRIVATE + ${RAT_MLFMM_DIR}/include +) + +target_link_libraries(rat_mlfmm_c PRIVATE ratmlfmm Rat::Common) + +set_target_properties(rat_mlfmm_c PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib + ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib + RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin + POSITION_INDEPENDENT_CODE ON +) + +target_compile_definitions(rat_mlfmm_c PRIVATE RAT_MLFMM_C_BUILD) diff --git a/bindings/rat-mlfmm-c/cmake/ArmadilloConfig.cmake b/bindings/rat-mlfmm-c/cmake/ArmadilloConfig.cmake new file mode 100644 index 00000000..85ded647 --- /dev/null +++ b/bindings/rat-mlfmm-c/cmake/ArmadilloConfig.cmake @@ -0,0 +1,4 @@ +set(ARMADILLO_INCLUDE_DIRS "${CMAKE_BINARY_DIR}/armadillo-build/tmp/include") +set(ARMADILLO_LIBRARY_DIRS "") +set(ARMADILLO_LIBRARIES armadillo) +set(Armadillo_VERSION "15.2.3") diff --git a/bindings/rat-mlfmm-c/cmake/ArmadilloConfigVersion.cmake b/bindings/rat-mlfmm-c/cmake/ArmadilloConfigVersion.cmake new file mode 100644 index 00000000..ba84befc --- /dev/null +++ b/bindings/rat-mlfmm-c/cmake/ArmadilloConfigVersion.cmake @@ -0,0 +1,10 @@ +set(PACKAGE_VERSION "15.2.3") + +if(PACKAGE_VERSION VERSION_LESS PACKAGE_FIND_VERSION) + set(PACKAGE_VERSION_COMPATIBLE FALSE) +else() + set(PACKAGE_VERSION_COMPATIBLE TRUE) + if(PACKAGE_VERSION VERSION_EQUAL PACKAGE_FIND_VERSION) + set(PACKAGE_VERSION_EXACT TRUE) + endif() +endif() diff --git a/bindings/rat-mlfmm-c/cmake/BoostConfig.cmake.in b/bindings/rat-mlfmm-c/cmake/BoostConfig.cmake.in new file mode 100644 index 00000000..68b42d16 --- /dev/null +++ b/bindings/rat-mlfmm-c/cmake/BoostConfig.cmake.in @@ -0,0 +1,57 @@ +set(Boost_VERSION 1.90.0) +set(Boost_VERSION_STRING "1.90.0") +set(Boost_FOUND TRUE) + +set(Boost_INCLUDE_DIR "@BOOST_INSTALL_DIR@/include") +set(Boost_INCLUDE_DIRS "@BOOST_INSTALL_DIR@/include") +set(_boost_lib_dir "@BOOST_INSTALL_DIR@/lib") + +find_library(Boost_FILESYSTEM_LIBRARY NAMES boost_filesystem PATHS "${_boost_lib_dir}" NO_DEFAULT_PATH) +find_library(Boost_IOSTREAMS_LIBRARY NAMES boost_iostreams PATHS "${_boost_lib_dir}" NO_DEFAULT_PATH) +find_library(Boost_THREAD_LIBRARY NAMES boost_thread PATHS "${_boost_lib_dir}" NO_DEFAULT_PATH) +find_library(Boost_CHRONO_LIBRARY NAMES boost_chrono PATHS "${_boost_lib_dir}" NO_DEFAULT_PATH) + +if(NOT Boost_FILESYSTEM_LIBRARY OR NOT Boost_IOSTREAMS_LIBRARY OR NOT Boost_THREAD_LIBRARY OR NOT Boost_CHRONO_LIBRARY) + set(Boost_FOUND FALSE) + return() +endif() + +if(NOT TARGET Boost::boost) + add_library(Boost::boost INTERFACE IMPORTED) + set_target_properties(Boost::boost PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}" + ) +endif() + +if(NOT TARGET Boost::filesystem) + add_library(Boost::filesystem STATIC IMPORTED) + set_target_properties(Boost::filesystem PROPERTIES + IMPORTED_LOCATION "${Boost_FILESYSTEM_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}" + ) +endif() + +if(NOT TARGET Boost::iostreams) + add_library(Boost::iostreams STATIC IMPORTED) + set_target_properties(Boost::iostreams PROPERTIES + IMPORTED_LOCATION "${Boost_IOSTREAMS_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}" + INTERFACE_LINK_LIBRARIES "z" + ) +endif() + +if(NOT TARGET Boost::thread) + add_library(Boost::thread STATIC IMPORTED) + set_target_properties(Boost::thread PROPERTIES + IMPORTED_LOCATION "${Boost_THREAD_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}" + ) +endif() + +if(NOT TARGET Boost::chrono) + add_library(Boost::chrono STATIC IMPORTED) + set_target_properties(Boost::chrono PROPERTIES + IMPORTED_LOCATION "${Boost_CHRONO_LIBRARY}" + INTERFACE_INCLUDE_DIRECTORIES "${Boost_INCLUDE_DIRS}" + ) +endif() diff --git a/bindings/rat-mlfmm-c/cmake/FindJsonCPP.cmake b/bindings/rat-mlfmm-c/cmake/FindJsonCPP.cmake new file mode 100644 index 00000000..156c587e --- /dev/null +++ b/bindings/rat-mlfmm-c/cmake/FindJsonCPP.cmake @@ -0,0 +1,51 @@ +if(TARGET jsoncpp_lib OR TARGET jsoncpp_static) + set(JsonCPP_FOUND TRUE) + set(JsonCPP_INCLUDE_DIR "${JSONCPP_SRC_DIR}/include") + if(NOT TARGET JsonCPP::JsonCPP) + add_library(JsonCPP::JsonCPP INTERFACE IMPORTED) + if(TARGET jsoncpp_lib) + set(_jsoncpp_target jsoncpp_lib) + else() + set(_jsoncpp_target jsoncpp_static) + endif() + set_target_properties(JsonCPP::JsonCPP PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${JSONCPP_SRC_DIR}/include" + INTERFACE_LINK_LIBRARIES ${_jsoncpp_target} + ) + unset(_jsoncpp_target) + endif() + return() +endif() + +find_package(PkgConfig) +pkg_check_modules(PC_JsonCPP QUIET JsonCPP) + +find_path(JsonCPP_INCLUDE_DIR + NAMES json/config.h + PATHS ${PC_JsonCPP_INCLUDE_DIRS} + HINTS $ENV{JSONCPP}/include /usr/local/include /opt/local/include /usr/include/jsoncpp +) + +find_library(JsonCPP_LIBRARY + NAMES jsoncpp + PATHS ${PC_JsonCPP_LIBRARY_DIR} + HINTS $ENV{JSONCPP}/lib /usr/local/lib /opt/local/lib /usr/lib +) + +set(JsonCPP_VERSION ${PC_JsonCPP_VERSION}) + +mark_as_advanced(JsonCPP_FOUND JsonCPP_INCLUDE_DIR JsonCPP_VERSION) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(JsonCPP + REQUIRED_VARS JsonCPP_INCLUDE_DIR JsonCPP_LIBRARY + VERSION_VAR JsonCPP_VERSION +) + +if(JsonCPP_FOUND AND NOT TARGET JsonCPP::JsonCPP) + add_library(JsonCPP::JsonCPP INTERFACE IMPORTED) + set_target_properties(JsonCPP::JsonCPP PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${JsonCPP_INCLUDE_DIR}" + INTERFACE_LINK_LIBRARIES "${JsonCPP_LIBRARY}" + ) +endif() diff --git a/bindings/rat-mlfmm-c/cmake/FindTCLAP.cmake b/bindings/rat-mlfmm-c/cmake/FindTCLAP.cmake new file mode 100644 index 00000000..9ee46a6a --- /dev/null +++ b/bindings/rat-mlfmm-c/cmake/FindTCLAP.cmake @@ -0,0 +1,42 @@ +if(DEFINED TCLAP_INCLUDE_DIR) + set(TCLAP_FOUND TRUE) + set(TCLAP_INCLUDE_DIRS "${TCLAP_INCLUDE_DIR}") + if(NOT TARGET TCLAP::TCLAP) + add_library(TCLAP::TCLAP INTERFACE IMPORTED) + set_target_properties(TCLAP::TCLAP PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${TCLAP_INCLUDE_DIRS}" + ) + endif() + return() +endif() + +find_package(PkgConfig) +pkg_check_modules(PC_TCLAP QUIET TCLAP) + +mark_as_advanced(TCLAP_FOUND TCLAP_INCLUDE_DIR TCLAP_VERSION) + +find_path(TCLAP_INCLUDE_DIR + NAMES tclap/MultiArg.h + PATHS ${TCLAP_PKGCONF_INCLUDE_DIR} +) + +set(TCLAP_VERSION ${PC_TCLAP_VERSION}) + +mark_as_advanced(TCLAP_FOUND TCLAP_INCLUDE_DIR TCLAP_VERSION) + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(TCLAP + REQUIRED_VARS TCLAP_INCLUDE_DIR + VERSION_VAR TCLAP_VERSION +) + +if(TCLAP_FOUND) + get_filename_component(TCLAP_INCLUDE_DIRS ${TCLAP_INCLUDE_DIR} DIRECTORY) +endif() + +if(TCLAP_FOUND AND NOT TARGET TCLAP::TCLAP) + add_library(TCLAP::TCLAP INTERFACE IMPORTED) + set_target_properties(TCLAP::TCLAP PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES ${TCLAP_INCLUDE_DIRS} + ) +endif() diff --git a/bindings/rat-mlfmm-c/include/rat_mlfmm_c.h b/bindings/rat-mlfmm-c/include/rat_mlfmm_c.h new file mode 100644 index 00000000..eb5326e1 --- /dev/null +++ b/bindings/rat-mlfmm-c/include/rat_mlfmm_c.h @@ -0,0 +1,83 @@ +#ifndef RAT_MLFMM_C_H +#define RAT_MLFMM_C_H + +#include + +#ifdef __cplusplus +extern "C" { +#endif + +#if defined(_WIN32) + #if defined(RAT_MLFMM_C_BUILD) + #define RAT_MLFMM_C_API __declspec(dllexport) + #else + #define RAT_MLFMM_C_API __declspec(dllimport) + #endif +#else + #define RAT_MLFMM_C_API __attribute__((visibility("default"))) +#endif + +typedef struct rat_mlfmm_context rat_mlfmm_context; + +typedef enum rat_mlfmm_direct_mode { + RAT_MLFMM_DIRECT_ALWAYS = 0, + RAT_MLFMM_DIRECT_THRESHOLD = 1, + RAT_MLFMM_DIRECT_NEVER = 2 +} rat_mlfmm_direct_mode; + +RAT_MLFMM_C_API rat_mlfmm_context *rat_mlfmm_context_create(void); +RAT_MLFMM_C_API void rat_mlfmm_context_destroy(rat_mlfmm_context *ctx); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_sources_linear( + rat_mlfmm_context *ctx, + const double *rs_x, + const double *rs_y, + const double *rs_z, + const double *drs_x, + const double *drs_y, + const double *drs_z, + const double *currents, + const double *eps, + size_t num_sources); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_targets( + rat_mlfmm_context *ctx, + const double *rt_x, + const double *rt_y, + const double *rt_z, + size_t num_targets); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_van_lanen( + rat_mlfmm_context *ctx, + int use_van_lanen); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_num_exp( + rat_mlfmm_context *ctx, + int num_exp); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_direct_mode( + rat_mlfmm_context *ctx, + rat_mlfmm_direct_mode mode); + +RAT_MLFMM_C_API int rat_mlfmm_context_set_direct_threshold( + rat_mlfmm_context *ctx, + double threshold); + +RAT_MLFMM_C_API int rat_mlfmm_context_compute_ba( + rat_mlfmm_context *ctx, + double *out_bx, + double *out_by, + double *out_bz, + size_t out_b_len, + double *out_ax, + double *out_ay, + double *out_az, + size_t out_a_len); + +RAT_MLFMM_C_API const char *rat_mlfmm_last_error(void); + +#ifdef __cplusplus +} +#endif + +#endif diff --git a/bindings/rat-mlfmm-c/src/rat_mlfmm_c.cpp b/bindings/rat-mlfmm-c/src/rat_mlfmm_c.cpp new file mode 100644 index 00000000..61519320 --- /dev/null +++ b/bindings/rat-mlfmm-c/src/rat_mlfmm_c.cpp @@ -0,0 +1,322 @@ +#include "rat_mlfmm_c.h" + +#include +#include +#include +#include +#include + +#include "currentsources.hh" +#include "mgntargets.hh" +#include "mlfmm.hh" +#include "settings.hh" + +namespace { +thread_local std::string g_last_error; +} + +extern "C" RAT_MLFMM_C_API void rat_mlfmm_set_last_error(const char *msg) { + if (msg) { + g_last_error = msg; + } else { + g_last_error.clear(); + } +} + +extern "C" RAT_MLFMM_C_API const char *rat_mlfmm_last_error(void) { + return g_last_error.c_str(); +} + +namespace { +struct Context { + rat::fmm::ShCurrentSourcesPr sources; + rat::fmm::ShMgnTargetsPr targets; + rat::fmm::ShSettingsPr settings; + rat::fmm::ShMlfmmPr mlfmm; + bool use_van_lanen = true; +}; + +static arma::Mat copy_mat_3xn_cols( + const double *x, + const double *y, + const double *z, + size_t n_cols) { + arma::Mat out(3, n_cols); + for (size_t j = 0; j < n_cols; ++j) { + out(0, j) = static_cast(x[j]); + out(1, j) = static_cast(y[j]); + out(2, j) = static_cast(z[j]); + } + return out; +} + +static arma::Row copy_row(const double *data, size_t n) { + arma::Row out(n); + for (size_t i = 0; i < n; ++i) { + out(i) = static_cast(data[i]); + } + return out; +} + +static int set_error_and_return(const char *msg) { + rat_mlfmm_set_last_error(msg); + return 0; +} + +static int handle_exception(const std::exception &ex) { + return set_error_and_return(ex.what()); +} +} // namespace + +extern "C" rat_mlfmm_context *rat_mlfmm_context_create(void) { + try { + auto *ctx = new Context(); + ctx->settings = rat::fmm::Settings::create(); + rat_mlfmm_set_last_error(nullptr); + return reinterpret_cast(ctx); + } catch (const std::exception &ex) { + handle_exception(ex); + } catch (...) { + set_error_and_return("unknown error in rat_mlfmm_context_create"); + } + return nullptr; +} + +extern "C" void rat_mlfmm_context_destroy(rat_mlfmm_context *ctx) { + if (!ctx) { + return; + } + auto *raw = reinterpret_cast(ctx); + delete raw; +} + +extern "C" int rat_mlfmm_context_set_sources_linear( + rat_mlfmm_context *ctx, + const double *rs_x, + const double *rs_y, + const double *rs_z, + const double *drs_x, + const double *drs_y, + const double *drs_z, + const double *currents, + const double *eps, + size_t num_sources) { + + if (!ctx || !rs_x || !rs_y || !rs_z || !drs_x || !drs_y || !drs_z || !currents || !eps) { + return set_error_and_return("null pointer in set_sources_linear"); + } + if (num_sources == 0) { + return set_error_and_return("num_sources must be positive"); + } + + try { + auto *raw = reinterpret_cast(ctx); + arma::Mat Rs = copy_mat_3xn_cols(rs_x, rs_y, rs_z, num_sources); + arma::Mat dRs = copy_mat_3xn_cols(drs_x, drs_y, drs_z, num_sources); + arma::Row Is = copy_row(currents, num_sources); + arma::Row epss = copy_row(eps, num_sources); + + // Rat kernels expect segment centers; shift from start-point inputs to midpoint. + Rs += RAT_CONST(0.5) * dRs; + raw->sources = rat::fmm::CurrentSources::create(Rs, dRs, Is, epss); + raw->sources->set_van_Lanen(raw->use_van_lanen); + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_sources_linear"); + } +} + +extern "C" int rat_mlfmm_context_set_targets( + rat_mlfmm_context *ctx, + const double *rt_x, + const double *rt_y, + const double *rt_z, + size_t num_targets) { + + if (!ctx || !rt_x || !rt_y || !rt_z) { + return set_error_and_return("null pointer in set_targets"); + } + if (num_targets == 0) { + return set_error_and_return("num_targets must be positive"); + } + + try { + auto *raw = reinterpret_cast(ctx); + arma::Mat Rt = copy_mat_3xn_cols(rt_x, rt_y, rt_z, num_targets); + raw->targets = rat::fmm::MgnTargets::create(Rt); + raw->targets->set_field_type('B', 3); + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_targets"); + } +} + +extern "C" int rat_mlfmm_context_set_van_lanen( + rat_mlfmm_context *ctx, + int use_van_lanen) { + + if (!ctx) { + return set_error_and_return("null pointer in set_van_lanen"); + } + + try { + auto *raw = reinterpret_cast(ctx); + raw->use_van_lanen = (use_van_lanen != 0); + if (raw->sources) { + raw->sources->set_van_Lanen(raw->use_van_lanen); + } + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_van_lanen"); + } +} + +extern "C" int rat_mlfmm_context_set_num_exp( + rat_mlfmm_context *ctx, + int num_exp) { + + if (!ctx) { + return set_error_and_return("null pointer in set_num_exp"); + } + if (num_exp <= 0) { + return set_error_and_return("num_exp must be positive"); + } + + try { + auto *raw = reinterpret_cast(ctx); + raw->settings->set_num_exp(num_exp); + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_num_exp"); + } +} + +extern "C" int rat_mlfmm_context_set_direct_mode( + rat_mlfmm_context *ctx, + rat_mlfmm_direct_mode mode) { + + if (!ctx) { + return set_error_and_return("null pointer in set_direct_mode"); + } + + try { + auto *raw = reinterpret_cast(ctx); + rat::fmm::DirectMode dm = rat::fmm::DirectMode::NEVER; + switch (mode) { + case RAT_MLFMM_DIRECT_ALWAYS: + dm = rat::fmm::DirectMode::ALWAYS; + break; + case RAT_MLFMM_DIRECT_THRESHOLD: + dm = rat::fmm::DirectMode::TRESHOLD; + break; + case RAT_MLFMM_DIRECT_NEVER: + dm = rat::fmm::DirectMode::NEVER; + break; + default: + return set_error_and_return("invalid direct mode"); + } + raw->settings->set_direct(dm); + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_direct_mode"); + } +} + +extern "C" int rat_mlfmm_context_set_direct_threshold( + rat_mlfmm_context *ctx, + double threshold) { + + if (!ctx) { + return set_error_and_return("null pointer in set_direct_threshold"); + } + if (threshold <= 0.0) { + return set_error_and_return("direct threshold must be positive"); + } + + try { + auto *raw = reinterpret_cast(ctx); + raw->settings->set_direct_tresh(static_cast(threshold)); + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in set_direct_threshold"); + } +} + +extern "C" int rat_mlfmm_context_compute_ba( + rat_mlfmm_context *ctx, + double *out_bx, + double *out_by, + double *out_bz, + size_t out_b_len, + double *out_ax, + double *out_ay, + double *out_az, + size_t out_a_len) { + + if (!ctx || !out_bx || !out_by || !out_bz || !out_ax || !out_ay || !out_az) { + return set_error_and_return("null pointer in compute_ba"); + } + + try { + auto *raw = reinterpret_cast(ctx); + if (!raw->sources) { + return set_error_and_return("sources not set"); + } + if (!raw->targets) { + return set_error_and_return("targets not set"); + } + + const size_t num_targets = raw->targets->num_targets(); + if (out_b_len < num_targets || out_a_len < num_targets) { + return set_error_and_return("output buffer too small"); + } + + raw->targets->set_field_type("BA", arma::Row{3, 3}); + raw->mlfmm = rat::fmm::Mlfmm::create(raw->sources, raw->targets, raw->settings); + raw->mlfmm->setup(); + raw->mlfmm->calculate(); + + const arma::Mat B = raw->targets->get_field('B'); + const arma::Mat A = raw->targets->get_field('A'); + if (B.n_rows != 3 || B.n_cols != num_targets) { + return set_error_and_return("unexpected B-field shape"); + } + if (A.n_rows != 3 || A.n_cols != num_targets) { + return set_error_and_return("unexpected A-field shape"); + } + + for (size_t j = 0; j < num_targets; ++j) { + out_bx[j] = static_cast(B(0, j)); + out_by[j] = static_cast(B(1, j)); + out_bz[j] = static_cast(B(2, j)); + out_ax[j] = static_cast(A(0, j)); + out_ay[j] = static_cast(A(1, j)); + out_az[j] = static_cast(A(2, j)); + } + + rat_mlfmm_set_last_error(nullptr); + return 1; + } catch (const std::exception &ex) { + return handle_exception(ex); + } catch (...) { + return set_error_and_return("unknown error in compute_ba"); + } +} diff --git a/bindings/repair_wheel.sh b/bindings/repair_wheel.sh new file mode 100755 index 00000000..2a8f5c42 --- /dev/null +++ b/bindings/repair_wheel.sh @@ -0,0 +1,78 @@ +#!/usr/bin/env bash + +# Cross-platform packaging of dynamically linked libraries +# for wheels. + +set -euo pipefail + +usage() { + cat <<'USAGE' +Usage: repair_wheel.sh [wheel_dir] [output_dir] + +Examples: + repair_wheel.sh dist + repair_wheel.sh dist repaired + repair_wheel.sh + +Defaults: + wheel_dir = dist + output_dir = wheel_dir +USAGE +} + +wheel_dir=${1:-dist} +output_dir=${2:-${wheel_dir}} + +if [ "${wheel_dir}" = "-h" ] || [ "${wheel_dir}" = "--help" ]; then + usage + exit 0 +fi + +shopt -s nullglob +wheels=("${wheel_dir}"/*.whl) +shopt -u nullglob + +if [ ${#wheels[@]} -eq 0 ]; then + echo "no wheels found in ${wheel_dir}" >&2 + exit 1 +fi + +case "$(uname -s)" in + Darwin) + if ! command -v delocate-wheel >/dev/null 2>&1; then + echo "delocate-wheel not found; install with: python -m pip install delocate" >&2 + exit 1 + fi + for wheel in "${wheels[@]}"; do + delocate-wheel -w "${output_dir}" "${wheel}" + done + ;; + Linux) + if ! command -v auditwheel >/dev/null 2>&1; then + echo "auditwheel not found; install with: python -m pip install auditwheel" >&2 + exit 1 + fi + for wheel in "${wheels[@]}"; do + if [[ "${wheel}" == *musllinux* ]]; then + # auditwheel targets glibc/manylinux and cannot parse musllinux wheels. + # musllinux wheels are expected to work as-built without auditwheel repair. + echo "Skipping auditwheel for musllinux wheel: ${wheel}" + continue + fi + auditwheel repair "${wheel}" -w "${output_dir}" + done + ;; + MINGW*|MSYS*|CYGWIN*) + if ! command -v delvewheel >/dev/null 2>&1; then + echo "delvewheel not found; install with: python -m pip install delvewheel" >&2 + exit 1 + fi + for wheel in "${wheels[@]}"; do + delvewheel repair "${wheel}" -w "${output_dir}" + done + ;; + *) + echo "Unsupported OS: $(uname -s)" >&2 + exit 1 + ;; +esac diff --git a/bindings/vendor-patches/rat-common-static.patch b/bindings/vendor-patches/rat-common-static.patch new file mode 100644 index 00000000..17986803 --- /dev/null +++ b/bindings/vendor-patches/rat-common-static.patch @@ -0,0 +1,26 @@ +diff --git a/CMakeLists.txt b/CMakeLists.txt +index 212cbb7..d800930 100644 +--- a/CMakeLists.txt ++++ b/CMakeLists.txt +@@ -12,6 +12,7 @@ option(ENABLE_EXAMPLES "build examples" ON) + option(ENABLE_BLAS "build with optimized BLAS library" ON) + option(ENABLE_SUPERLU "link with superlu" ON) + option(ENABLE_MKL_ALLOC "enable memory allocation with mkl" OFF) ++option(RAT_COMMON_BUILD_SHARED "Build ratcmn as shared library" OFF) + + # compile commands for clangd + set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +@@ -179,7 +180,12 @@ set(source_list + ) + + # add the common library +-add_library(ratcmn SHARED ${source_list}) ++if(RAT_COMMON_BUILD_SHARED) ++ set(_rat_common_lib_type SHARED) ++else() ++ set(_rat_common_lib_type STATIC) ++endif() ++add_library(ratcmn ${_rat_common_lib_type} ${source_list}) + + # Add an alias so that library can be used inside the build tree, e.g. when testing + add_library(Rat::Common ALIAS ratcmn) diff --git a/bindings/vendor-patches/rat-mlfmm-static.patch b/bindings/vendor-patches/rat-mlfmm-static.patch new file mode 100644 index 00000000..42d30bf6 --- /dev/null +++ b/bindings/vendor-patches/rat-mlfmm-static.patch @@ -0,0 +1,26 @@ +diff --git a/CMakeLists.txt b/CMakeLists.txt +index d6951e8..6563147 100644 +--- a/CMakeLists.txt ++++ b/CMakeLists.txt +@@ -12,6 +12,7 @@ option(ENABLE_CUDA_FAST_MATH "use fast math for cuda routines" OFF) + option(ENABLE_TESTING "build unit/system tests" ON) + option(ENABLE_EXAMPLES "build examples" ON) + option(ENABLE_MATLAB "build matlab mex interface" OFF) ++option(RAT_MLFMM_BUILD_SHARED "Build ratmlfmm as shared library" OFF) + + # where are the custom CMake modules located + list(INSERT CMAKE_MODULE_PATH 0 ${CMAKE_CURRENT_SOURCE_DIR}/cmake) +@@ -123,7 +124,12 @@ if(CMAKE_CUDA_COMPILER) + endif() + + # add the library +-add_library(ratmlfmm SHARED ${source_list}) ++if(RAT_MLFMM_BUILD_SHARED) ++ set(_rat_mlfmm_lib_type SHARED) ++else() ++ set(_rat_mlfmm_lib_type STATIC) ++endif() ++add_library(ratmlfmm ${_rat_mlfmm_lib_type} ${source_list}) + + # Add an alias so that library can be used inside the build tree, e.g. when testing + add_library(Rat::MLFMM ALIAS ratmlfmm) diff --git a/build.rs b/build.rs new file mode 100644 index 00000000..11d40e59 --- /dev/null +++ b/build.rs @@ -0,0 +1,613 @@ +use std::env; +use std::io; +use std::path::{Path, PathBuf}; +use std::process::Command; + +fn main() { + if env::var("CARGO_FEATURE_RAT_MLFMM").is_err() { + return; + } + + let python_feature = env::var_os("CARGO_FEATURE_PYTHON").is_some(); + + let manifest_dir = PathBuf::from(env::var("CARGO_MANIFEST_DIR").unwrap()); + let patches_dir = manifest_dir.join("bindings").join("vendor-patches"); + let wrapper_dir = manifest_dir.join("bindings").join("rat-mlfmm-c"); + let rat_mlfmm_dir = manifest_dir.join("vendor").join("rat-mlfmm"); + let rat_common_dir = manifest_dir.join("vendor").join("rat-common"); + let jsoncpp_dir = manifest_dir.join("vendor").join("jsoncpp"); + let armadillo_dir = manifest_dir.join("vendor").join("armadillo-15.2.3"); + let armadillo_zip = manifest_dir.join("vendor").join("armadillo-15.2.3.zip"); + let tclap_dir = manifest_dir.join("vendor").join("tclap"); + let boost_dir = manifest_dir.join("vendor").join("boost-boost-1.90.0"); + + ensure_submodules( + &manifest_dir, + &[ + wrapper_dir.join("CMakeLists.txt"), + rat_mlfmm_dir.join("CMakeLists.txt"), + rat_common_dir.join("CMakeLists.txt"), + jsoncpp_dir.join("CMakeLists.txt"), + tclap_dir.join("CMakeLists.txt"), + boost_dir.join("tools/build/src/engine/build.sh"), + ], + ); + ensure_boost_headers(&boost_dir); + ensure_armadillo_extracted(&armadillo_dir, &armadillo_zip); + apply_vendor_patches(&manifest_dir, &patches_dir, &rat_common_dir, &rat_mlfmm_dir); + + let out_dir = PathBuf::from(env::var("OUT_DIR").unwrap()); + let target_arch = env::var("CARGO_CFG_TARGET_ARCH").unwrap_or_default(); + let target_os = env::var("CARGO_CFG_TARGET_OS").unwrap_or_default(); + let cpu_flag = resolve_cpu_flag(&target_arch); + let extra_c_flags = env::var("CFLAGS").ok(); + let extra_cxx_flags = env::var("CXXFLAGS").ok(); + let (mut c_flags_release, mut c_flags_debug) = + compose_c_flags(&cpu_flag, extra_c_flags.as_deref()); + let (mut cxx_flags_release, mut cxx_flags_debug) = + compose_cxx_flags(&cpu_flag, extra_cxx_flags.as_deref()); + let build_profile = env::var("PROFILE").unwrap_or_else(|_| "release".to_string()); + let is_release = build_profile == "release"; + if target_os == "linux" { + ensure_tool("clang"); + ensure_tool("clang++"); + ensure_tool("ld.lld"); + } + + let rat_mlfmm_c_src = wrapper_dir.join("src").join("rat_mlfmm_c.cpp"); + let rat_mlfmm_c_include = wrapper_dir.join("include"); + let rat_mlfmm_include = rat_mlfmm_dir.join("include"); + let rat_common_include = rat_common_dir.join("include"); + let boost_include = boost_dir.clone(); + let rat_common_shim = out_dir.join("rat-common-include"); + let rat_common_shim_dest = rat_common_shim.join("rat").join("common"); + std::fs::create_dir_all(&rat_common_shim_dest) + .unwrap_or_else(|err| panic!("failed to create {rat_common_shim_dest:?}: {err}")); + if let Ok(entries) = std::fs::read_dir(&rat_common_include) { + for entry in entries.flatten() { + let path = entry.path(); + if path.extension().and_then(|ext| ext.to_str()) != Some("hh") { + continue; + } + if let Some(name) = path.file_name() { + let dest = rat_common_shim_dest.join(name); + let _ = std::fs::copy(&path, dest); + } + } + } + let armadillo_include = armadillo_dir.join("include"); + let jsoncpp_include = jsoncpp_dir.join("include"); + let mut cc_build = cc::Build::new(); + cc_build.cargo_metadata(false); // We emit rat_mlfmm_c link directives ourselves to avoid duplicate archives. + cc_build.cpp(true); + if target_os == "linux" { + cc_build.compiler("clang++"); + } + cc_build.file(&rat_mlfmm_c_src); + cc_build.include(&rat_mlfmm_c_include); + cc_build.include(&rat_mlfmm_include); + cc_build.include(&rat_common_shim); + cc_build.include(&armadillo_include); + cc_build.include(&jsoncpp_include); + cc_build.include(&boost_include); + cc_build.define("RAT_MLFMM_C_BUILD", None); + cc_build.define("RAT_DOUBLE_PRECISION", None); + cc_build.flag_if_supported("-std=c++14"); + cc_build.flag_if_supported("-fPIC"); + if is_release { + for flag in cxx_flags_release.split_whitespace() { + cc_build.flag(flag); + } + } else { + for flag in cxx_flags_debug.split_whitespace() { + cc_build.flag(flag); + } + } + cc_build.out_dir(&out_dir); + cc_build.compile("rat_mlfmm_c"); + + let mut cfg = cmake::Config::new(&wrapper_dir); + let profile = env::var("PROFILE").unwrap_or_else(|_| "release".to_string()); + let build_type = if profile == "release" { + "Release" + } else { + "Debug" + }; + let boost_cxxflags = compose_boost_cxxflags(&cpu_flag, extra_cxx_flags.as_deref()); + cfg.profile(build_type); + cfg.define("CMAKE_BUILD_TYPE", build_type); + cfg.define("CMAKE_C_FLAGS_RELEASE", &c_flags_release); + cfg.define("CMAKE_CXX_FLAGS_RELEASE", &cxx_flags_release); + cfg.define("CMAKE_C_FLAGS_DEBUG", &c_flags_debug); + cfg.define("CMAKE_CXX_FLAGS_DEBUG", &cxx_flags_debug); + if target_os == "linux" { + cfg.define("CMAKE_C_COMPILER", "clang"); + cfg.define("CMAKE_CXX_COMPILER", "clang++"); + cfg.define("CMAKE_EXE_LINKER_FLAGS", "-fuse-ld=lld"); + cfg.define("CMAKE_SHARED_LINKER_FLAGS", "-fuse-ld=lld"); + cfg.define("CMAKE_MODULE_LINKER_FLAGS", "-fuse-ld=lld"); + } + cfg.define("CMAKE_POSITION_INDEPENDENT_CODE", "ON"); + cfg.define("CFSEM_BOOST_CXXFLAGS", &boost_cxxflags); + cfg.define( + "CFSEM_EXPECT_RELEASE", + if build_type == "Release" { "ON" } else { "OFF" }, + ); + if build_type == "Release" { + cfg.define("CMAKE_INTERPROCEDURAL_OPTIMIZATION", "ON"); + } + if target_os == "macos" { + let deploy = env::var("MACOSX_DEPLOYMENT_TARGET").unwrap_or_else(|_| "11.0".to_string()); + cfg.define("CMAKE_OSX_DEPLOYMENT_TARGET", &deploy); + } else { + cfg.define("BLA_STATIC", "ON"); + if target_os == "windows" { + cfg.define("CMAKE_FIND_LIBRARY_SUFFIXES", ".lib;.a"); + } else if target_os == "linux" { + cfg.define("CMAKE_FIND_LIBRARY_SUFFIXES", ".a"); + } + } + cfg.define("RAT_MLFMM_DIR", rat_mlfmm_dir.to_str().unwrap()); + cfg.define("RAT_COMMON_DIR", rat_common_dir.to_str().unwrap()); + cfg.define("JSONCPP_SRC_DIR", jsoncpp_dir.to_str().unwrap()); + cfg.define("ARMADILLO_SRC_DIR", armadillo_dir.to_str().unwrap()); + cfg.define("TCLAP_SRC_DIR", tclap_dir.to_str().unwrap()); + cfg.define("BOOST_SRC_DIR", boost_dir.to_str().unwrap()); + if let Ok(prefix) = env::var("CMAKE_PREFIX_PATH") { + cfg.define("CMAKE_PREFIX_PATH", &prefix); + } + if let Ok(jsoncpp_dir) = env::var("JSONCPP_DIR") { + cfg.define("JsonCPP_DIR", &jsoncpp_dir); + } + if let Ok(tclap_dir) = env::var("TCLAP_DIR") { + cfg.define("TCLAP_DIR", &tclap_dir); + } + if let Ok(armadillo_dir) = env::var("ARMADILLO_DIR") { + cfg.define("Armadillo_DIR", &armadillo_dir); + } + + let dst = cfg.build(); + + let lib_dir = dst.join("build").join("lib"); + let bin_dir = dst.join("build").join("bin"); + let rat_common_lib_dir = dst.join("build").join("rat-common-build").join("lib"); + let rat_mlfmm_lib_dir = dst.join("build").join("rat-mlfmm-build").join("lib"); + let rat_common_lib = rat_common_lib_dir.join("libratcmn.a"); + let rat_mlfmm_lib = rat_mlfmm_lib_dir.join("libratmlfmm.a"); + let rat_mlfmm_c_lib = out_dir.join("librat_mlfmm_c.a"); + let boost_lib_dir = dst.join("build").join("boost-install").join("lib"); + let armadillo_lib_dir = dst.join("build").join("armadillo-build"); + println!("cargo:rustc-link-search=native={}", out_dir.display()); + println!("cargo:rustc-link-search=native={}", lib_dir.display()); + println!( + "cargo:rustc-link-search=native={}", + rat_common_lib_dir.display() + ); + println!( + "cargo:rustc-link-search=native={}", + rat_mlfmm_lib_dir.display() + ); + println!("cargo:rustc-link-search=native={}", boost_lib_dir.display()); + println!( + "cargo:rustc-link-search=native={}", + armadillo_lib_dir.display() + ); + if bin_dir.exists() { + println!("cargo:rustc-link-search=native={}", bin_dir.display()); + } + let mut static_libs = vec![ + "rat_mlfmm_c", + "boost_filesystem", + "boost_iostreams", + "boost_thread", + "boost_chrono", + "jsoncpp", + "armadillo", + ]; + if python_feature && target_os == "linux" { + static_libs.retain(|lib| *lib != "rat_mlfmm_c"); + } + if target_os == "linux" && !python_feature { + static_libs.insert(1, "ratmlfmm"); + static_libs.insert(2, "ratcmn"); + } + if target_os == "linux" && !python_feature { + emit_link_group(&static_libs); + } else { + for lib in &static_libs { + println!("cargo:rustc-link-lib=static={lib}"); + } + } + if !python_feature { + println!("cargo:rustc-link-lib=static=ratmlfmm"); + println!("cargo:rustc-link-lib=static=ratcmn"); + } else if target_os == "linux" { + println!("cargo:rustc-link-arg-cdylib=-Wl,--no-as-needed"); + println!("cargo:rustc-link-arg-cdylib=-Wl,--whole-archive"); + println!("cargo:rustc-link-arg-cdylib={}", rat_mlfmm_c_lib.display()); + println!("cargo:rustc-link-arg-cdylib={}", rat_mlfmm_lib.display()); + println!("cargo:rustc-link-arg-cdylib={}", rat_common_lib.display()); + println!("cargo:rustc-link-arg-cdylib=-Wl,--no-whole-archive"); + } + println!("cargo:rustc-link-lib=z"); + if target_os == "linux" { + println!("cargo:rustc-link-lib=stdc++"); + println!("cargo:rustc-link-lib=openblas"); + } + if target_os == "macos" { + println!( + "cargo:rustc-link-arg-cdylib=-Wl,-force_load,{}", + rat_mlfmm_c_lib.display() + ); + println!( + "cargo:rustc-link-arg-cdylib=-Wl,-force_load,{}", + rat_mlfmm_lib.display() + ); + println!( + "cargo:rustc-link-arg-cdylib=-Wl,-force_load,{}", + rat_common_lib.display() + ); + println!("cargo:rustc-link-lib=c++"); + println!("cargo:rustc-link-lib=c++abi"); + println!("cargo:rustc-link-lib=framework=Accelerate"); + } + let _ = lib_dir; + + rerun_if_changed(&wrapper_dir.join("CMakeLists.txt")); + rerun_if_changed(&wrapper_dir.join("src/rat_mlfmm_c.cpp")); + rerun_if_changed(&wrapper_dir.join("include/rat_mlfmm_c.h")); + rerun_if_changed(&wrapper_dir.join("cmake/FindJsonCPP.cmake")); + rerun_if_changed(&wrapper_dir.join("cmake/FindTCLAP.cmake")); + rerun_if_changed(&wrapper_dir.join("cmake/ArmadilloConfig.cmake")); + rerun_if_changed(&rat_common_dir.join("CMakeLists.txt")); + rerun_if_changed(&rat_mlfmm_dir.join("CMakeLists.txt")); + rerun_if_changed(&jsoncpp_dir.join("CMakeLists.txt")); + rerun_if_changed(&armadillo_dir.join("CMakeLists.txt")); + rerun_if_changed(&armadillo_zip); + rerun_if_changed(&tclap_dir.join("CMakeLists.txt")); + rerun_if_changed(&boost_dir.join("CMakeLists.txt")); + rerun_if_changed(&wrapper_dir.join("cmake/BoostConfig.cmake.in")); + rerun_if_changed(&wrapper_dir.join("cmake/ArmadilloConfig.cmake")); + rerun_if_changed(&patches_dir.join("rat-common-static.patch")); + rerun_if_changed(&patches_dir.join("rat-mlfmm-static.patch")); +} + +fn rerun_if_changed(path: &Path) { + println!("cargo:rerun-if-changed={}", path.display()); +} + +fn resolve_cpu_flag(target_arch: &str) -> Option { + let mut target_cpu = extract_target_cpu_from_rustflags(); + if target_cpu.is_none() && matches!(target_arch, "x86_64" | "x64") { + target_cpu = Some("x86-64-v3".to_string()); + } + + target_cpu.map(|cpu| { + if matches!(target_arch, "x86_64" | "x64") { + format!("-march={cpu}") + } else if matches!(target_arch, "aarch64") { + format!("-mcpu={cpu}") + } else { + format!("-march={cpu}") + } + }) +} + +fn extract_target_cpu_from_rustflags() -> Option { + let encoded = env::var("CARGO_ENCODED_RUSTFLAGS").ok()?; + let parts: Vec<&str> = encoded.split('\u{1f}').collect(); + for i in 0..parts.len() { + let part = parts[i]; + if let Some(value) = part.strip_prefix("-Ctarget-cpu=") { + return Some(value.to_string()); + } + if part == "-C" { + if let Some(next) = parts.get(i + 1) { + if let Some(value) = next.strip_prefix("target-cpu=") { + return Some(value.to_string()); + } + } + } + } + None +} + +fn compose_c_flags(cpu_flag: &Option, extra_flags: Option<&str>) -> (String, String) { + if cfg!(target_os = "windows") { + return ("/O2".to_string(), "/O2 /Zi".to_string()); + } + let mut release = vec!["-O3"]; + let mut debug = vec!["-O3", "-g"]; + if let Some(flag) = cpu_flag { + release.push(flag); + debug.push(flag); + } + if cfg!(not(target_os = "windows")) { + release.push("-fPIC"); + debug.push("-fPIC"); + } + if let Some(extra) = extra_flags { + release.push(extra); + debug.push(extra); + } + (release.join(" "), debug.join(" ")) +} + +fn compose_cxx_flags(cpu_flag: &Option, extra_flags: Option<&str>) -> (String, String) { + compose_c_flags(cpu_flag, extra_flags) +} + +fn compose_boost_cxxflags(cpu_flag: &Option, extra_flags: Option<&str>) -> String { + if cfg!(target_os = "windows") { + return "/O2".to_string(); + } + let mut flags = vec!["-O3"]; + if let Some(flag) = cpu_flag { + flags.push(flag); + } + if cfg!(not(target_os = "windows")) { + flags.push("-fPIC"); + } + if let Some(extra) = extra_flags { + flags.push(extra); + } + flags.join(" ") +} + +fn emit_link_group(static_libs: &[&str]) { + println!("cargo:rustc-link-arg=-Wl,--start-group"); + for lib in static_libs { + println!("cargo:rustc-link-lib=static={lib}"); + } + println!("cargo:rustc-link-arg=-Wl,--end-group"); +} + +fn ensure_submodules(manifest_dir: &Path, required_paths: &[PathBuf]) { + let missing: Vec<_> = required_paths + .iter() + .filter(|path| !path.exists()) + .collect(); + + if missing.is_empty() { + return; + } + + if env::var_os("CI").is_some() || env::var_os("GITHUB_ACTIONS").is_some() { + let _ = Command::new("git") + .args([ + "config", + "--global", + "--add", + "safe.directory", + manifest_dir.to_str().unwrap(), + ]) + .status(); + let vendor_root = manifest_dir.join("vendor"); + let _ = Command::new("git") + .args(["config", "--global", "--add", "safe.directory", "*"]) + .status(); + if let Ok(entries) = std::fs::read_dir(&vendor_root) { + for entry in entries.flatten() { + let path = entry.path(); + if path.is_dir() { + if let Some(dir) = path.to_str() { + let _ = Command::new("git") + .args(["config", "--global", "--add", "safe.directory", dir]) + .status(); + } + } + } + } + } + + let status = Command::new("git") + .args([ + "submodule", + "update", + "--init", + "--recursive", + "--depth", + "1", + ]) + .current_dir(manifest_dir) + .status() + .expect("failed to run git submodule update"); + + if !status.success() { + panic!("git submodule update failed with status {status}"); + } + + for path in required_paths { + if !path.exists() { + panic!( + "required path missing after submodule update: {}", + path.display() + ); + } + } +} + +fn apply_vendor_patches( + manifest_dir: &Path, + patches_dir: &Path, + rat_common_dir: &Path, + rat_mlfmm_dir: &Path, +) { + let patches = [ + ("rat-common-static.patch", rat_common_dir), + ("rat-mlfmm-static.patch", rat_mlfmm_dir), + ]; + + for (patch_name, repo_dir) in patches { + let patch = patches_dir.join(patch_name); + if !patch.exists() { + panic!("required patch missing: {}", patch.display()); + } + if !repo_dir.exists() { + panic!("required vendor dir missing: {}", repo_dir.display()); + } + + let check = Command::new("git") + .args(["apply", "--check", patch.to_str().unwrap()]) + .current_dir(repo_dir) + .status() + .expect("failed to run git apply --check"); + + if check.success() { + let status = Command::new("git") + .args(["apply", patch.to_str().unwrap()]) + .current_dir(repo_dir) + .status() + .expect("failed to run git apply"); + if !status.success() { + panic!( + "failed to apply patch {} in {}", + patch.display(), + repo_dir.display() + ); + } + continue; + } + + let reverse_check = Command::new("git") + .args(["apply", "--reverse", "--check", patch.to_str().unwrap()]) + .current_dir(repo_dir) + .status() + .expect("failed to run git apply --reverse --check"); + + if !reverse_check.success() { + panic!( + "patch {} does not apply cleanly in {}", + patch.display(), + repo_dir.display() + ); + } + } + + let _ = manifest_dir; +} + +fn ensure_boost_headers(boost_dir: &Path) { + let boost_headers = boost_dir.join("boost"); + if boost_headers.is_dir() { + return; + } + // Boost's git superproject does not include the generated `boost/` header tree. + // Running `bootstrap.sh` + `b2 headers` creates it so includes like work. + let target_os = env::var("CARGO_CFG_TARGET_OS").unwrap_or_default(); + // Windows uses bootstrap.bat and b2.exe; Unix uses bootstrap.sh and b2. + let bootstrap = if target_os == "windows" { + boost_dir.join("bootstrap.bat") + } else { + boost_dir.join("bootstrap.sh") + }; + if !bootstrap.exists() { + panic!( + "Boost headers missing and bootstrap.sh not found at {}", + bootstrap.display() + ); + } + let status = if target_os == "windows" { + Command::new("cmd") + .current_dir(boost_dir) + .args(["/C", bootstrap.to_str().unwrap()]) + .status() + } else { + Command::new(&bootstrap).current_dir(boost_dir).status() + } + .unwrap_or_else(|err| panic!("failed to run {}: {err}", bootstrap.display())); + if !status.success() { + panic!("Boost bootstrap failed with status {}", status); + } + // b2 headers generates the consolidated boost/ headers directory from the repo layout. + let b2 = if target_os == "windows" { + boost_dir.join("b2.exe") + } else { + boost_dir.join("b2") + }; + // Generate the consolidated headers into `boost/` for downstream compilers. + let status = if target_os == "windows" { + Command::new("cmd") + .current_dir(boost_dir) + .args(["/C", b2.to_str().unwrap(), "headers"]) + .status() + } else { + Command::new(&b2) + .current_dir(boost_dir) + .arg("headers") + .status() + } + .unwrap_or_else(|err| panic!("failed to run {}: {err}", b2.display())); + if !status.success() { + panic!("Boost header generation failed with status {}", status); + } + if !boost_headers.is_dir() { + panic!( + "Boost header generation completed but {} is still missing", + boost_headers.display() + ); + } +} + +fn ensure_armadillo_extracted(armadillo_dir: &Path, armadillo_zip: &Path) { + if armadillo_dir.join("CMakeLists.txt").exists() { + return; + } + if !armadillo_zip.exists() { + panic!("armadillo zip missing: {}", armadillo_zip.display()); + } + + std::fs::create_dir_all(armadillo_dir).expect("failed to create armadillo directory"); + + let file = std::fs::File::open(armadillo_zip) + .unwrap_or_else(|err| panic!("failed to open armadillo zip: {err}")); + let mut archive = zip::ZipArchive::new(file) + .unwrap_or_else(|err| panic!("failed to read armadillo zip: {err}")); + + for i in 0..archive.len() { + let mut entry = archive + .by_index(i) + .unwrap_or_else(|err| panic!("failed to read armadillo zip entry: {err}")); + let name = entry.name().to_string(); + let stripped = name.splitn(2, '/').nth(1).unwrap_or(""); + if stripped.is_empty() { + continue; + } + let outpath = armadillo_dir.join(stripped); + if entry.is_dir() { + std::fs::create_dir_all(&outpath) + .unwrap_or_else(|err| panic!("failed to create dir {outpath:?}: {err}")); + continue; + } + + if let Some(parent) = outpath.parent() { + std::fs::create_dir_all(parent) + .unwrap_or_else(|err| panic!("failed to create dir {parent:?}: {err}")); + } + + let mut outfile = std::fs::File::create(&outpath) + .unwrap_or_else(|err| panic!("failed to create file {outpath:?}: {err}")); + io::copy(&mut entry, &mut outfile) + .unwrap_or_else(|err| panic!("failed to extract {outpath:?}: {err}")); + + #[cfg(unix)] + { + use std::os::unix::fs::PermissionsExt; + if let Some(mode) = entry.unix_mode() { + let mut perms = outfile + .metadata() + .unwrap_or_else(|err| panic!("failed to stat {outpath:?}: {err}")) + .permissions(); + perms.set_mode(mode); + std::fs::set_permissions(&outpath, perms) + .unwrap_or_else(|err| panic!("failed to set perms {outpath:?}: {err}")); + } + } + } +} + +fn ensure_tool(tool: &str) { + if Command::new(tool).arg("--version").output().is_err() { + panic!("{tool} not found in PATH"); + } +} diff --git a/cfsem/__init__.py b/cfsem/__init__.py index 2046e9e4..ee645099 100644 --- a/cfsem/__init__.py +++ b/cfsem/__init__.py @@ -9,13 +9,14 @@ from cfsem.bindings import ( body_force_density_circular_filament_cartesian, body_force_density_linear_filament, + fields_linear_filament_mlfmm, filament_helix_path, flux_circular_filament, - flux_density_biot_savart, flux_density_circular_filament, flux_density_circular_filament_cartesian, flux_density_dipole, flux_density_linear_filament, + flux_density_point_segment, gs_operator_order2, gs_operator_order4, inductance_piecewise_linear_filaments, @@ -23,6 +24,7 @@ rotate_filaments_about_path, vector_potential_circular_filament, vector_potential_linear_filament, + vector_potential_point_segment, vector_potential_dipole, ) from cfsem.types import Array3xN @@ -37,7 +39,6 @@ __all__ = [ "flux_circular_filament", - "flux_density_biot_savart", "flux_density_linear_filament", "flux_density_circular_filament", "gs_operator_order2", @@ -61,9 +62,12 @@ "vector_potential_linear_filament", "vector_potential_circular_filament", "flux_density_circular_filament_cartesian", + "fields_linear_filament_mlfmm", "mutual_inductance_circular_to_linear", "flux_density_dipole", "vector_potential_dipole", + "flux_density_point_segment", + "vector_potential_point_segment", "body_force_density_circular_filament_cartesian", "body_force_density_linear_filament", ] diff --git a/cfsem/bindings.py b/cfsem/bindings.py index 5ccb8436..b0b40b72 100644 --- a/cfsem/bindings.py +++ b/cfsem/bindings.py @@ -5,7 +5,7 @@ passed as contiguous and reallocating into contiguous inputs if necessary. """ -from numpy import ascontiguousarray, float64, zeros_like +from numpy import asarray, ascontiguousarray, float64, full, zeros_like from numpy.typing import NDArray from cfsem.types import Array3xN @@ -25,6 +25,7 @@ from .cfsem import flux_density_dipole as em_flux_density_dipole from .cfsem import vector_potential_dipole as em_vector_potential_dipole from .cfsem import flux_density_linear_filament as em_flux_density_linear_filament +from .cfsem import flux_density_point_segment as em_flux_density_point_segment from .cfsem import gs_operator_order2 as em_gs_operator_order2 from .cfsem import gs_operator_order4 as em_gs_operator_order4 from .cfsem import ( @@ -40,6 +41,16 @@ from .cfsem import ( vector_potential_linear_filament as em_vector_potential_linear_filament, ) +from .cfsem import ( + vector_potential_point_segment as em_vector_potential_point_segment, +) + +try: + from .cfsem import ( + fields_linear_filament_mlfmm as em_fields_linear_filament_mlfmm, + ) +except ImportError: # pragma: no cover - depends on optional rat-mlfmm feature + em_fields_linear_filament_mlfmm = None def flux_circular_filament( @@ -176,6 +187,7 @@ def flux_density_linear_filament( xyzfil: Array3xN, dlxyzfil: Array3xN, ifil: NDArray[float64], + wire_radius: float | NDArray[float64] = 0.0, par: bool = True, ) -> Array3xN: """ @@ -184,9 +196,10 @@ def flux_density_linear_filament( Args: xyzp: [m] x,y,z coords of observation points - xyzfil: [m] x,y,z coords of current filament origins (start of segment) - dlxyzfil: [m] x,y,z length delta of current filaments + xyzfil: [m] x,y,z coords of filament segment start points + dlxyzfil: [m] x,y,z deltas from segment start to segment end ifil: [A] current in each filament segment + wire_radius: [m] filament radius, scalar or array of length `m` par: Whether to use CPU parallelism Returns: @@ -196,10 +209,38 @@ def flux_density_linear_filament( xyzfil = _3tup_contig(xyzfil) dlxyzfil = _3tup_contig(dlxyzfil) ifil = ascontiguousarray(ifil).ravel() - return em_flux_density_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, par) + if asarray(wire_radius).ndim == 0: + wire_radius = full(ifil.size, float(wire_radius)) + wire_radius = ascontiguousarray(wire_radius).ravel() + return em_flux_density_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par) + + +def flux_density_point_segment( + xyzp: Array3xN, + xyzfil: Array3xN, + dlxyzfil: Array3xN, + ifil: NDArray[float64], + par: bool = True, +) -> Array3xN: + """ + Biot-Savart law calculation for B-field contributions from many filament segments + to many observation points, treating each segment as a point source. + Args: + xyzp: [m] x,y,z coords of observation points + xyzfil: [m] x,y,z coords of filament segment start points + dlxyzfil: [m] x,y,z deltas from segment start to segment end + ifil: [A] current in each filament segment + par: Whether to use CPU parallelism -flux_density_biot_savart = flux_density_linear_filament # For backwards-compatibility + Returns: + [T] (Bx, By, Bz) magnetic flux density at observation points + """ + xyzp = _3tup_contig(xyzp) + xyzfil = _3tup_contig(xyzfil) + dlxyzfil = _3tup_contig(dlxyzfil) + ifil = ascontiguousarray(ifil).ravel() + return em_flux_density_point_segment(xyzp, xyzfil, dlxyzfil, ifil, par) def vector_potential_linear_filament( @@ -207,6 +248,7 @@ def vector_potential_linear_filament( xyzfil: Array3xN, dlxyzfil: Array3xN, ifil: NDArray[float64], + wire_radius: float | NDArray[float64] = 0.0, par: bool = True, ) -> Array3xN: """ @@ -215,8 +257,40 @@ def vector_potential_linear_filament( Args: xyzp: [m] x,y,z coords of observation points - xyzfil: [m] x,y,z coords of current filament origins (start of segment) - dlxyzfil: [m] x,y,z length delta of current filaments + xyzfil: [m] x,y,z coords of filament segment start points + dlxyzfil: [m] x,y,z deltas from segment start to segment end + ifil: [A] current in each filament segment + wire_radius: [m] filament radius, scalar or array of length `m` + par: Whether to use CPU parallelism + + Returns: + [Wb/m] or [V-s/m] (Ax, Ay, Az) magnetic vector potential at observation points + """ + xyzp = _3tup_contig(xyzp) + xyzfil = _3tup_contig(xyzfil) + dlxyzfil = _3tup_contig(dlxyzfil) + ifil = ascontiguousarray(ifil).ravel() + if asarray(wire_radius).ndim == 0: + wire_radius = full(ifil.size, float(wire_radius)) + wire_radius = ascontiguousarray(wire_radius).ravel() + return em_vector_potential_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par) + + +def vector_potential_point_segment( + xyzp: Array3xN, + xyzfil: Array3xN, + dlxyzfil: Array3xN, + ifil: NDArray[float64], + par: bool = True, +) -> Array3xN: + """ + Vector potential calculation for A-field contribution from many filament + segments to many observation points, treating each segment as a point source. + + Args: + xyzp: [m] x,y,z coords of observation points + xyzfil: [m] x,y,z coords of filament segment start points + dlxyzfil: [m] x,y,z deltas from segment start to segment end ifil: [A] current in each filament segment par: Whether to use CPU parallelism @@ -227,7 +301,54 @@ def vector_potential_linear_filament( xyzfil = _3tup_contig(xyzfil) dlxyzfil = _3tup_contig(dlxyzfil) ifil = ascontiguousarray(ifil).ravel() - return em_vector_potential_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, par) + return em_vector_potential_point_segment(xyzp, xyzfil, dlxyzfil, ifil, par) + + +def fields_linear_filament_mlfmm( + xyzp: Array3xN, + xyzfil: Array3xN, + dlxyzfil: Array3xN, + ifil: NDArray[float64], + eps: NDArray[float64], + *, + use_linear_filament: bool = True, + direct_threshold: int = 10_000_000, + order: int | None = None, +) -> tuple[Array3xN, Array3xN]: + """ + MLFMM calculation for B- and A-field contributions from many filament segments + to many observation points. + + Args: + xyzp: [m] x,y,z coords of observation points + xyzfil: [m] x,y,z coords of filament segment start points + dlxyzfil: [m] x,y,z deltas from segment start to segment end + ifil: [A] current in each filament segment + eps: [m] van Lanen softening parameter + use_linear_filament: Whether to use linear-filament kernel + direct_threshold: Interaction count threshold for direct evaluation + order: Multipole expansion order (None uses library default) + + Returns: + (B, A) tuples of [T] and [Wb/m] field components at observation points + """ + if em_fields_linear_filament_mlfmm is None: + raise RuntimeError("rat-mlfmm feature is not enabled in this build") # pragma: no cover + xyzp = _3tup_contig(xyzp) + xyzfil = _3tup_contig(xyzfil) + dlxyzfil = _3tup_contig(dlxyzfil) + ifil = ascontiguousarray(ifil).ravel() + eps = ascontiguousarray(eps).ravel() + return em_fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + eps, + use_linear_filament, + direct_threshold, + 0 if order is None else int(order), + ) def inductance_piecewise_linear_filaments( @@ -570,6 +691,7 @@ def body_force_density_linear_filament( ifil: NDArray[float64], obs: Array3xN, j: Array3xN, + wire_radius: float | NDArray[float64] = 0.0, par: bool = True, ) -> Array3xN: """ @@ -582,6 +704,7 @@ def body_force_density_linear_filament( ifil: [A] filament current obs: [m] x,y,z coords of observation locations j: [A/m^2] current density vector at observation locations + wire_radius: [m] filament radius, scalar or array of length `m` par: Whether to use CPU parallelism Returns: @@ -592,7 +715,12 @@ def body_force_density_linear_filament( ifil = ascontiguousarray(ifil).ravel() obs = _3tup_contig(obs) j = _3tup_contig(j) - jxbx, jxby, jxbz = em_body_force_density_linear_filament(xyzfil, dlxyzfil, ifil, obs, j, par) # [N/m^3] + if asarray(wire_radius).ndim == 0: + wire_radius = full(ifil.size, float(wire_radius)) + wire_radius = ascontiguousarray(wire_radius).ravel() + jxbx, jxby, jxbz = em_body_force_density_linear_filament( + xyzfil, dlxyzfil, ifil, obs, j, wire_radius, par + ) # [N/m^3] return jxbx, jxby, jxbz # type: ignore diff --git a/examples/biot_savart.py b/examples/biot_savart.py new file mode 100644 index 00000000..38c6700c --- /dev/null +++ b/examples/biot_savart.py @@ -0,0 +1,266 @@ +from __future__ import annotations + +import os +import time + +import numpy as np + +if os.getenv("CFSEM_TESTING"): + import matplotlib + + matplotlib.use("Agg") + +import matplotlib.pyplot as plt + +import cfsem + + +def main() -> None: + # Single filament from z=-0.5 to z=0.5 along the z-axis. + xyzfil = (np.array([0.0]), np.array([0.0]), np.array([-0.5])) + dlxyzfil = (np.array([0.0]), np.array([0.0]), np.array([1.0])) + ifil = np.array([1.0]) + + # Sample plane: x-z plane at y=0 to show end effects. + n = 21 if os.getenv("CFSEM_TESTING") else 2001 + x = np.linspace(-1.0, 1.0, n) + z = np.linspace(-1.0, 1.0, n) + xx, zz = np.meshgrid(x, z, indexing="xy") + yy = np.zeros_like(xx) + + xyzp = (xx.ravel(), yy.ravel(), zz.ravel()) + + # Use a small wire radius to avoid singularities on-axis. + wire_radius = 0.01 + t0 = time.perf_counter() + bx, by, bz = cfsem.flux_density_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius=wire_radius, par=True + ) + t_linear = time.perf_counter() - t0 + + bmag = np.sqrt(bx * bx + by * by + bz * bz).reshape(xx.shape) + bmag_log10 = np.log10(bmag + 1e-30) + + # Discretize into point-segment sources for comparison. + nseg = 1000 + dz = 1.0 / nseg + xfil_ps = np.zeros(nseg) + yfil_ps = np.zeros(nseg) + zfil_ps = np.linspace(-0.5, 0.5 - dz, nseg) + dlx_ps = np.zeros(nseg) + dly_ps = np.zeros(nseg) + dlz_ps = np.full(nseg, dz) + ifil_ps = np.full(nseg, ifil[0]) + xyzfil_ps = (xfil_ps, yfil_ps, zfil_ps) + dlxyzfil_ps = (dlx_ps, dly_ps, dlz_ps) + + t0 = time.perf_counter() + bx_ps, by_ps, bz_ps = cfsem.flux_density_point_segment( + xyzp, xyzfil_ps, dlxyzfil_ps, ifil_ps, par=True + ) + t_point = time.perf_counter() - t0 + bmag_ps = np.sqrt(bx_ps * bx_ps + by_ps * by_ps + bz_ps * bz_ps).reshape(xx.shape) + + fig, axs = plt.subplots( + 4, + 3, + figsize=(12, 8), + dpi=120, + gridspec_kw={"width_ratios": [1.1, 1.0, 1.0]}, + ) + ax_map, ax_line_x, ax_line_z = axs[0] + ax_err_map, ax_err_x, ax_err_z = axs[1] + ax_mlfmm_map, ax_mlfmm_x, ax_mlfmm_z = axs[2] + ax_lf_map, ax_lf_x, ax_lf_z = axs[3] + im = ax_map.imshow( + bmag_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="magma", + aspect="equal", + ) + ax_map.contour( + xx, + zz, + bmag_log10, + levels=50, + colors="k", + linewidths=0.6, + alpha=0.6, + ) + + ax_map.plot([0.0, 0.0], [-0.5, 0.5], color="white", linewidth=2.0) + ax_map.set_xlabel("x [m]") + ax_map.set_ylabel("z [m]") + ax_map.set_title("Linear Filament\nB-field magnitude (log10)") + cbar = fig.colorbar(im, ax=ax_map) + cbar.set_label("log10(|B|) [T]") + + mid_idx = n // 2 + ax_line_x.plot(x, bmag[mid_idx, :], color="black", label="linear") + ax_line_x.plot( + x, bmag_ps[mid_idx, :], color="cyan", linestyle="--", label="point segment" + ) + ax_line_x.set_xlabel("x [m]") + ax_line_x.set_ylabel("|B| [T]") + ax_line_x.set_title("Slice along x (z = 0)") + ax_line_x.set_ylim(0.0, np.max(bmag[mid_idx, :])) + ax_line_x.grid(True, alpha=0.3) + ax_line_x.legend(frameon=False) + + ax_line_z.plot(z, bmag[:, mid_idx], color="black", label="linear") + ax_line_z.plot( + z, bmag_ps[:, mid_idx], color="cyan", linestyle="--", label="point segment" + ) + ax_line_z.set_xlabel("z [m]") + ax_line_z.set_ylabel("|B| [T]") + ax_line_z.set_title("Slice along z (x = 0)") + ax_line_z.grid(True, alpha=0.3) + ax_line_z.legend(frameon=False) + + err = np.abs(bmag - bmag_ps) + mask = np.abs(xx) < wire_radius + err = np.where(mask, np.nan, err) + err_log10 = np.log10(err) + + im_err = ax_err_map.imshow( + err_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_err_map.set_xlabel("x [m]") + ax_err_map.set_ylabel("z [m]") + ax_err_map.set_title("Point-segment discretization\nerror magnitude (log10)") + cbar_err = fig.colorbar(im_err, ax=ax_err_map) + cbar_err.set_label("log10(|ΔB|) [T]") + + rel_err_x = err[mid_idx, :] / (bmag[mid_idx, :] + 1e-30) + ax_err_x.plot(x, rel_err_x, color="black") + ax_err_x.set_xlabel("x [m]") + ax_err_x.set_ylabel("relative error") + ax_err_x.set_title("Point-segment discretization\nrelative error slice along x (z = 0)") + ax_err_x.grid(True, alpha=0.3) + + ax_err_z.set_axis_off() + + try: + t0 = time.perf_counter() + b_mlfmm, _ = cfsem.fields_linear_filament_mlfmm( + xyzp, + xyzfil_ps, + dlxyzfil_ps, + ifil_ps, + np.full(nseg, 1e-6), + use_linear_filament=False, + direct_threshold=1, # Always MLFMM to test multipole expansion + order=None, + ) + t_mlfmm = time.perf_counter() - t0 + bx_m, by_m, bz_m = b_mlfmm + bmag_m = np.sqrt(bx_m * bx_m + by_m * by_m + bz_m * bz_m).reshape(xx.shape) + err_m = np.abs(bmag_m - bmag_ps) + err_m = np.where(mask, np.nan, err_m) + err_m_log10 = np.log10(err_m) + + im_m = ax_mlfmm_map.imshow( + err_m_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_mlfmm_map.set_xlabel("x [m]") + ax_mlfmm_map.set_ylabel("z [m]") + ax_mlfmm_map.set_title("MLFMM (point-segment expansion)\nerror vs point segment (log10)") + cbar_m = fig.colorbar(im_m, ax=ax_mlfmm_map) + cbar_m.set_label("log10(|ΔB|) [T]") + + rel_err_m_x = err_m[mid_idx, :] / (bmag_ps[mid_idx, :] + 1e-30) + ax_mlfmm_x.plot(x, rel_err_m_x, color="black") + ax_mlfmm_x.set_xlabel("x [m]") + ax_mlfmm_x.set_ylabel("relative error") + ax_mlfmm_x.set_title("MLFMM (point-segment expansion)\nrelative error slice along x (z = 0)") + ax_mlfmm_x.grid(True, alpha=0.3) + + ax_line_z.plot( + z, + bmag_m[:, mid_idx], + color="lime", + linestyle=":", + label="mlfmm", + ) + ax_line_z.legend(frameon=False) + + ax_mlfmm_z.set_axis_off() + + t0 = time.perf_counter() + b_mlfmm_lf, _ = cfsem.fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + np.full(ifil.size, wire_radius), + use_linear_filament=True, + direct_threshold=1, # Still direct method in this case + order=None, + ) + t_mlfmm_lf = time.perf_counter() - t0 + bx_lf, by_lf, bz_lf = b_mlfmm_lf + bmag_lf = np.sqrt(bx_lf * bx_lf + by_lf * by_lf + bz_lf * bz_lf).reshape(xx.shape) + err_lf = np.abs(bmag_lf - bmag) + # err_lf = np.where(mask, np.nan, err_lf) + err_lf_log10 = np.log10(err_lf + 1e-30) + + im_lf = ax_lf_map.imshow( + err_lf_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_lf_map.set_xlabel("x [m]") + ax_lf_map.set_ylabel("z [m]") + ax_lf_map.set_title("MLFMM (linear-filament direct)\nerror vs linear (log10)") + cbar_lf = fig.colorbar(im_lf, ax=ax_lf_map) + cbar_lf.set_label("log10(|ΔB|) [T]") + + rel_err_lf_x = err_lf[mid_idx, :] / (bmag[mid_idx, :] + 1e-30) + ax_lf_x.plot(x, rel_err_lf_x, color="black") + ax_lf_x.set_xlabel("x [m]") + ax_lf_x.set_ylabel("relative error") + ax_lf_x.set_title("MLFMM (linear-filament direct)\nrelative error slice along x (z = 0)") + ax_lf_x.grid(True, alpha=0.3) + + ax_lf_z.plot(z, err_lf[:, mid_idx], color="black") + ax_lf_z.set_xlabel("z [m]") + ax_lf_z.set_ylabel("|ΔB| [T]") + ax_lf_z.set_title("MLFMM (linear-filament direct)\nerror slice along z (x = 0)") + ax_lf_z.grid(True, alpha=0.3) + + n_targets = xyzp[0].size + n_linear = ifil.size * n_targets + n_point = ifil_ps.size * n_targets + n_mlfmm = ifil_ps.size * n_targets + n_mlfmm_lf = ifil.size * n_targets + print(f"Linear filament (cfsem direct): {t_linear:.3f} s ({n_linear:.1e} interactions)") + print(f"Point segment (cfsem direct): {t_point:.3f} s ({n_point:.1e} interactions)") + print(f"MLFMM (point-segment exp.): {t_mlfmm:.3f} s ({n_mlfmm:.1e} interactions)") + print(f"MLFMM (linear-filament direct): {t_mlfmm_lf:.3f} s ({n_mlfmm_lf:.1e} interactions)") + except RuntimeError: + n_targets = xyzp[0].size + n_linear = ifil.size * n_targets + n_point = ifil_ps.size * n_targets + print(f"Linear filament (cfsem direct): {t_linear:.3f} s ({n_linear:.1e} interactions)") + print(f"Point segment (cfsem direct): {t_point:.3f} s ({n_point:.1e} interactions)") + for ax in (ax_mlfmm_map, ax_mlfmm_x, ax_mlfmm_z, ax_lf_map, ax_lf_x, ax_lf_z): + ax.text(0.5, 0.5, "MLFMM not available", ha="center", va="center") + ax.set_axis_off() + + fig.tight_layout() + plt.show() + + +if __name__ == "__main__": + main() diff --git a/examples/helmholtz.py b/examples/helmholtz.py index 6046a9e6..92bc86c2 100644 --- a/examples/helmholtz.py +++ b/examples/helmholtz.py @@ -1,8 +1,15 @@ """Calculate the B-field from a Helmholtz coil pair.""" import os +import time import numpy as np + +if os.getenv("CFSEM_TESTING"): + import matplotlib + + matplotlib.use("Agg") + from matplotlib import pyplot as plt from cfsem import MU_0, flux_density_circular_filament @@ -27,9 +34,10 @@ zmesh_flat = zmesh.flatten() # Calculate the B-field at every mesh point using cfsem. -Br_flat, Bz_flat = flux_density_circular_filament( - ifil, rfil, zfil, rmesh_flat, zmesh_flat -) +t0 = time.perf_counter() +Br_flat, Bz_flat = flux_density_circular_filament(ifil, rfil, zfil, rmesh_flat, zmesh_flat) +elapsed = time.perf_counter() - t0 +print(f"Computed Helmholtz coil field at {rmesh.size} locations in {elapsed:.6f} s") Br = Br_flat.reshape(rmesh.shape) diff --git a/examples/vector_potential.py b/examples/vector_potential.py new file mode 100644 index 00000000..96f7df21 --- /dev/null +++ b/examples/vector_potential.py @@ -0,0 +1,265 @@ +from __future__ import annotations + +import os +import time + +import numpy as np + +if os.getenv("CFSEM_TESTING"): + import matplotlib + + matplotlib.use("Agg") + +import matplotlib.pyplot as plt + +import cfsem + + +def main() -> None: + # Single filament from z=-0.5 to z=0.5 along the z-axis. + xyzfil = (np.array([0.0]), np.array([0.0]), np.array([-0.5])) + dlxyzfil = (np.array([0.0]), np.array([0.0]), np.array([1.0])) + ifil = np.array([1.0]) + + # Sample plane: x-z plane at y=0 to show end effects. + n = 21 if os.getenv("CFSEM_TESTING") else 2001 + x = np.linspace(-1.0, 1.0, n) + z = np.linspace(-1.0, 1.0, n) + xx, zz = np.meshgrid(x, z, indexing="xy") + yy = np.zeros_like(xx) + xyzp = (xx.ravel(), yy.ravel(), zz.ravel()) + + wire_radius = 0.01 + t0 = time.perf_counter() + ax, ay, az = cfsem.vector_potential_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius=wire_radius, par=True + ) + t_linear = time.perf_counter() - t0 + amag = np.sqrt(ax * ax + ay * ay + az * az).reshape(xx.shape) + amag_log10 = np.log10(amag + 1e-30) + + # Discretize into point-segment sources for comparison. + nseg = 1000 + dz = 1.0 / nseg + xfil_ps = np.zeros(nseg) + yfil_ps = np.zeros(nseg) + zfil_ps = np.linspace(-0.5, 0.5 - dz, nseg) + dlx_ps = np.zeros(nseg) + dly_ps = np.zeros(nseg) + dlz_ps = np.full(nseg, dz) + ifil_ps = np.full(nseg, ifil[0]) + xyzfil_ps = (xfil_ps, yfil_ps, zfil_ps) + dlxyzfil_ps = (dlx_ps, dly_ps, dlz_ps) + + t0 = time.perf_counter() + ax_ps, ay_ps, az_ps = cfsem.vector_potential_point_segment( + xyzp, xyzfil_ps, dlxyzfil_ps, ifil_ps, par=True + ) + t_point = time.perf_counter() - t0 + amag_ps = np.sqrt(ax_ps * ax_ps + ay_ps * ay_ps + az_ps * az_ps).reshape(xx.shape) + + fig, axs = plt.subplots( + 4, + 3, + figsize=(12, 8), + dpi=120, + gridspec_kw={"width_ratios": [1.1, 1.0, 1.0]}, + ) + ax_map, ax_line_x, ax_line_z = axs[0] + ax_err_map, ax_err_x, ax_err_z = axs[1] + ax_mlfmm_map, ax_mlfmm_x, ax_mlfmm_z = axs[2] + ax_lf_map, ax_lf_x, ax_lf_z = axs[3] + + im = ax_map.imshow( + amag_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="magma", + aspect="equal", + ) + ax_map.contour( + xx, + zz, + amag_log10, + levels=50, + colors="k", + linewidths=0.6, + alpha=0.6, + ) + ax_map.plot([0.0, 0.0], [-0.5, 0.5], color="white", linewidth=2.0) + ax_map.set_xlabel("x [m]") + ax_map.set_ylabel("z [m]") + ax_map.set_title("Linear Filament\nA-field magnitude (log10)") + cbar = fig.colorbar(im, ax=ax_map) + cbar.set_label("log10(|A|) [T m]") + + mid_idx = n // 2 + ax_line_x.plot(x, amag[mid_idx, :], color="black", label="linear") + ax_line_x.plot( + x, amag_ps[mid_idx, :], color="cyan", linestyle="--", label="point segment" + ) + ax_line_x.set_xlabel("x [m]") + ax_line_x.set_ylabel("|A| [T m]") + ax_line_x.set_title("Slice along x (z = 0)") + ax_line_x.set_ylim(0.0, np.max(amag[mid_idx, :])) + ax_line_x.grid(True, alpha=0.3) + ax_line_x.legend(frameon=False) + + ax_line_z.plot(z, amag[:, mid_idx], color="black", label="linear") + ax_line_z.plot( + z, amag_ps[:, mid_idx], color="cyan", linestyle="--", label="point segment" + ) + ax_line_z.set_xlabel("z [m]") + ax_line_z.set_ylabel("|A| [T m]") + ax_line_z.set_title("Slice along z (x = 0)") + ax_line_z.grid(True, alpha=0.3) + ax_line_z.legend(frameon=False) + + err = np.abs(amag - amag_ps) + mask = np.abs(xx) < wire_radius + err = np.where(mask, np.nan, err) + err_log10 = np.log10(err) + + im_err = ax_err_map.imshow( + err_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_err_map.set_xlabel("x [m]") + ax_err_map.set_ylabel("z [m]") + ax_err_map.set_title("Point-segment discretization\nerror magnitude (log10)") + cbar_err = fig.colorbar(im_err, ax=ax_err_map) + cbar_err.set_label("log10(|ΔA|) [T m]") + + rel_err_x = err[mid_idx, :] / (amag[mid_idx, :] + 1e-30) + ax_err_x.plot(x, rel_err_x, color="black") + ax_err_x.set_xlabel("x [m]") + ax_err_x.set_ylabel("relative error") + ax_err_x.set_title("Point-segment discretization\nrelative error slice along x (z = 0)") + ax_err_x.grid(True, alpha=0.3) + + ax_err_z.set_axis_off() + + try: + t0 = time.perf_counter() + _, a_mlfmm = cfsem.fields_linear_filament_mlfmm( + xyzp, + xyzfil_ps, + dlxyzfil_ps, + ifil_ps, + np.full(nseg, 1e-6), + use_linear_filament=False, + direct_threshold=1, + order=None, + ) + t_mlfmm = time.perf_counter() - t0 + ax_m, ay_m, az_m = a_mlfmm + amag_m = np.sqrt(ax_m * ax_m + ay_m * ay_m + az_m * az_m).reshape(xx.shape) + err_m = np.abs(amag_m - amag_ps) + err_m = np.where(mask, np.nan, err_m) + err_m_log10 = np.log10(err_m) + + im_m = ax_mlfmm_map.imshow( + err_m_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_mlfmm_map.set_xlabel("x [m]") + ax_mlfmm_map.set_ylabel("z [m]") + ax_mlfmm_map.set_title("MLFMM (point-segment expansion)\nerror vs point segment (log10)") + cbar_m = fig.colorbar(im_m, ax=ax_mlfmm_map) + cbar_m.set_label("log10(|ΔA|) [T m]") + + rel_err_m_x = err_m[mid_idx, :] / (amag_ps[mid_idx, :] + 1e-30) + ax_mlfmm_x.plot(x, rel_err_m_x, color="black") + ax_mlfmm_x.set_xlabel("x [m]") + ax_mlfmm_x.set_ylabel("relative error") + ax_mlfmm_x.set_title("MLFMM (point-segment expansion)\nrelative error slice along x (z = 0)") + ax_mlfmm_x.grid(True, alpha=0.3) + + ax_line_z.plot( + z, + amag_m[:, mid_idx], + color="lime", + linestyle=":", + label="mlfmm", + ) + ax_line_z.legend(frameon=False) + + ax_mlfmm_z.set_axis_off() + + t0 = time.perf_counter() + _, a_mlfmm_lf = cfsem.fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + np.full(ifil.size, wire_radius), + use_linear_filament=True, + direct_threshold=1, + order=None, + ) + t_mlfmm_lf = time.perf_counter() - t0 + ax_lf_m, ay_lf_m, az_lf_m = a_mlfmm_lf + amag_lf = np.sqrt(ax_lf_m * ax_lf_m + ay_lf_m * ay_lf_m + az_lf_m * az_lf_m).reshape( + xx.shape + ) + err_lf = np.abs(amag_lf - amag) + # err_lf = np.where(mask, np.nan, err_lf) + err_lf_log10 = np.log10(err_lf + 1e-30) + + im_lf = ax_lf_map.imshow( + err_lf_log10, + extent=(x.min(), x.max(), z.min(), z.max()), + origin="lower", + cmap="viridis", + aspect="equal", + ) + ax_lf_map.set_xlabel("x [m]") + ax_lf_map.set_ylabel("z [m]") + ax_lf_map.set_title("MLFMM (linear-filament direct)\nerror vs linear (log10)") + cbar_lf = fig.colorbar(im_lf, ax=ax_lf_map) + cbar_lf.set_label("log10(|ΔA|) [T m]") + + rel_err_lf_x = err_lf[mid_idx, :] / (amag[mid_idx, :] + 1e-30) + ax_lf_x.plot(x, rel_err_lf_x, color="black") + ax_lf_x.set_xlabel("x [m]") + ax_lf_x.set_ylabel("relative error") + ax_lf_x.set_title("MLFMM (linear-filament direct)\nrelative error slice along x (z = 0)") + ax_lf_x.grid(True, alpha=0.3) + + ax_lf_z.plot(z, err_lf[:, mid_idx], color="black") + ax_lf_z.set_xlabel("z [m]") + ax_lf_z.set_ylabel("|ΔA| [T m]") + ax_lf_z.set_title("MLFMM (linear-filament direct)\nerror slice along z (x = 0)") + ax_lf_z.grid(True, alpha=0.3) + + n_targets = xyzp[0].size + n_linear = ifil.size * n_targets + n_point = ifil_ps.size * n_targets + n_mlfmm = ifil_ps.size * n_targets + n_mlfmm_lf = ifil.size * n_targets + print(f"Linear filament (cfsem direct): {t_linear:.3f} s ({n_linear:.1e} interactions)") + print(f"Point segment (cfsem direct): {t_point:.3f} s ({n_point:.1e} interactions)") + print(f"MLFMM (point-segment exp.): {t_mlfmm:.3f} s ({n_mlfmm:.1e} interactions)") + print(f"MLFMM (linear-filament direct): {t_mlfmm_lf:.3f} s ({n_mlfmm_lf:.1e} interactions)") + except RuntimeError: + n_targets = xyzp[0].size + n_linear = ifil.size * n_targets + n_point = ifil_ps.size * n_targets + print(f"Linear filament (cfsem direct): {t_linear:.3f} s ({n_linear:.1e} interactions)") + print(f"Point segment (cfsem direct): {t_point:.3f} s ({n_point:.1e} interactions)") + for ax in (ax_mlfmm_map, ax_mlfmm_x, ax_mlfmm_z, ax_lf_map, ax_lf_x, ax_lf_z): + ax.text(0.5, 0.5, "MLFMM not available", ha="center", va="center") + ax.set_axis_off() + + fig.tight_layout() + plt.show() + + +if __name__ == "__main__": + main() diff --git a/pyproject.toml b/pyproject.toml index 57af54dc..21a09c84 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "maturin" [project] name = "cfsem" -version = "3.1.0" +dynamic = ["version"] description = "Quasi-steady electromagnetics including filamentized approximations, Biot-Savart, and Grad-Shafranov." authors = [{name = "Commonwealth Fusion Systems", email = "jlogan@cfs.energy"}] requires-python = ">=3.10, <3.14" @@ -15,7 +15,7 @@ classifiers = [ ] dependencies = [ "numpy >= 2", - "interpn[pydantic]>=0.8.2,<0.9", + "interpn[pydantic]>=0.8.2,<0.12", "findiff >= 0.12.1", "pydantic >= 2", "pydantic-numpy >= 6" @@ -39,11 +39,13 @@ dev = [ ] [tool.maturin] -features = ["pyo3/extension-module", "python"] +features = ["python", "rat-mlfmm"] module-name = "cfsem.cfsem" +version = "cargo" [tool.ruff] line-length = 110 +exclude = ["vendor"] [tool.ruff.lint] select = [ @@ -61,3 +63,6 @@ select = [ [tool.coverage.report] fail_under = 100 + +[tool.pytest.ini_options] +norecursedirs = ["vendor", "target", ".venv", "dist", "build"] diff --git a/src/lib.rs b/src/lib.rs index af584e8e..7563fb3d 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -17,6 +17,9 @@ pub mod math; pub mod mesh; pub mod physics; +#[cfg(feature = "rat-mlfmm")] +pub mod mlfmm; + #[cfg(test)] pub(crate) mod testing; @@ -33,6 +36,8 @@ pub(crate) fn chunksize(nelem: usize) -> usize { .unwrap_or(NonZeroUsize::MIN) .get(); + let ncores = ncores / 2; // Heuristic for physical cores + (nelem / ncores).max(1) } diff --git a/src/math.rs b/src/math.rs index 49283a87..7f93366e 100644 --- a/src/math.rs +++ b/src/math.rs @@ -152,7 +152,7 @@ pub fn cylindrical_to_cartesian(r: f64, phi: f64, z: f64) -> (f64, f64, f64) { /// Decompose two filament endpoints into a midpoint and a length vector #[inline] -pub fn decompose_filament( +pub(crate) fn decompose_filament( start: (f64, f64, f64), end: (f64, f64, f64), ) -> ((f64, f64, f64), (f64, f64, f64)) { @@ -167,6 +167,95 @@ pub fn decompose_filament( (midpoint, dl) } +pub(crate) struct PointLineDistance { + pub(crate) perp: f64, + pub(crate) dist_a: f64, + pub(crate) dist_b: f64, + pub(crate) frac: f64, + pub(crate) para_a: f64, + pub(crate) para_b: f64, + pub(crate) ab_norm: (f64, f64, f64), +} + +/// Minimum perpendicular distance to the infinite line defined by endpoints, +/// distances to each endpoint, clamp fraction based on `r_min`, +/// and parallel distances from each endpoint to the target. +pub(crate) fn point_line_distance_with_endpoints( + a: (f64, f64, f64), + b: (f64, f64, f64), + p: (f64, f64, f64), + r_min: f64, +) -> PointLineDistance { + // Vectors and distances between points. + let ab = (b.0 - a.0, b.1 - a.1, b.2 - a.2); + let ap = (p.0 - a.0, p.1 - a.1, p.2 - a.2); + let bp = (p.0 - b.0, p.1 - b.1, p.2 - b.2); + + let dist_a_raw = rss3(ap.0, ap.1, ap.2); + let dist_b_raw = rss3(bp.0, bp.1, bp.2); + + // Normalized segment vector. + // This might be zero, and that will be handled as late as possible to avoid disrupting + // calculations in nominal non-zero-length cases. + let ab2 = dot3(ab.0, ab.1, ab.2, ab.0, ab.1, ab.2); // (m^2) squared length. + let ab_len_inv = ab2.sqrt().recip(); + let ab_norm = (ab.0 * ab_len_inv, ab.1 * ab_len_inv, ab.2 * ab_len_inv); + + // Find the closest point on the infinite line defined by this segment to the target point. + let t = dot3(ap.0, ap.1, ap.2, ab.0, ab.1, ab.2) / ab2; + let closest = ( + t.mul_add(ab.0, a.0), + t.mul_add(ab.1, a.1), + t.mul_add(ab.2, a.2), + ); + let dp = (p.0 - closest.0, p.1 - closest.1, p.2 - closest.2); // (m) Vector from target to infinite line. + let perp_raw = rss3(dp.0, dp.1, dp.2); // (m) Un-clamped perpendicular distance. + + let r_min = r_min.max(0.0); + let r_min_frac = r_min.max(f64::MIN_POSITIVE); + + // Fraction of perpendicular distance to r_min, to be used for handling + // fields inside finite-thickness wires. + let frac = (perp_raw / r_min_frac).min(1.0); + + // Clamped parallel, perpendicular, and direct distances from segment to target. + let perp = perp_raw.max(r_min); + let dist_a = dist_a_raw.max(r_min); + let dist_b = dist_b_raw.max(r_min); + + let para_a = dot3(ap.0, ap.1, ap.2, ab.0, ab.1, ab.2) * ab_len_inv; + let para_b = dot3(bp.0, bp.1, bp.2, ab.0, ab.1, ab.2) * ab_len_inv; + + // Handle zero-length special case. + if ab2 == 0.0 { + let r_min = r_min.max(0.0); + let r_min_frac = r_min.max(f64::MIN_POSITIVE); + let frac = (dist_a_raw / r_min_frac).min(1.0); + let dist_a = dist_a_raw.max(r_min); + let dist_b = dist_b_raw.max(r_min); + let perp = dist_a; + return PointLineDistance { + perp, + dist_a, + dist_b, + frac, + para_a: 0.0, + para_b: 0.0, + ab_norm: (0.0, 0.0, 0.0), + }; + } + + PointLineDistance { + perp, + dist_a, + dist_b, + frac, + para_a, + para_b, + ab_norm, + } +} + /// Clip NaN values to the provided value. /// This is carefully organized to avoid producing any `jmp` /// instructions on modern-as-of-2025 systems. diff --git a/src/mlfmm/ffi.rs b/src/mlfmm/ffi.rs new file mode 100644 index 00000000..a5acb756 --- /dev/null +++ b/src/mlfmm/ffi.rs @@ -0,0 +1,71 @@ +use std::os::raw::{c_char, c_int}; + +#[repr(C)] +pub struct RatMlfmmContext { + _private: [u8; 0], +} + +#[repr(C)] +#[derive(Copy, Clone, Debug)] +pub enum RatMlfmmDirectMode { + Always = 0, + Threshold = 1, + Never = 2, +} + +unsafe extern "C" { + pub fn rat_mlfmm_context_create() -> *mut RatMlfmmContext; + pub fn rat_mlfmm_context_destroy(ctx: *mut RatMlfmmContext); + + pub fn rat_mlfmm_context_set_sources_linear( + ctx: *mut RatMlfmmContext, + rs_x: *const f64, + rs_y: *const f64, + rs_z: *const f64, + drs_x: *const f64, + drs_y: *const f64, + drs_z: *const f64, + currents: *const f64, + eps: *const f64, + num_sources: usize, + ) -> c_int; + + pub fn rat_mlfmm_context_set_targets( + ctx: *mut RatMlfmmContext, + rt_x: *const f64, + rt_y: *const f64, + rt_z: *const f64, + num_targets: usize, + ) -> c_int; + + pub fn rat_mlfmm_context_set_van_lanen( + ctx: *mut RatMlfmmContext, + use_van_lanen: c_int, + ) -> c_int; + + pub fn rat_mlfmm_context_set_num_exp(ctx: *mut RatMlfmmContext, num_exp: c_int) -> c_int; + + pub fn rat_mlfmm_context_set_direct_mode( + ctx: *mut RatMlfmmContext, + mode: RatMlfmmDirectMode, + ) -> c_int; + + pub fn rat_mlfmm_context_set_direct_threshold( + ctx: *mut RatMlfmmContext, + threshold: f64, + ) -> c_int; + + pub fn rat_mlfmm_context_compute_ba( + ctx: *mut RatMlfmmContext, + out_bx: *mut f64, + out_by: *mut f64, + out_bz: *mut f64, + out_b_len: usize, + out_ax: *mut f64, + out_ay: *mut f64, + out_az: *mut f64, + out_a_len: usize, + ) -> c_int; + + pub fn rat_mlfmm_last_error() -> *const c_char; +} diff --git a/src/mlfmm/mod.rs b/src/mlfmm/mod.rs new file mode 100644 index 00000000..ef74e5c8 --- /dev/null +++ b/src/mlfmm/mod.rs @@ -0,0 +1,520 @@ +//! Limited bindings to Project Rat's rat-mlfmm C++ library. + +mod ffi; + +use crate::macros::check_length_3tup; +use std::ffi::CStr; + +struct Context { + raw: *mut ffi::RatMlfmmContext, +} + +impl Context { + pub fn new() -> Result { + let raw = unsafe { ffi::rat_mlfmm_context_create() }; + if raw.is_null() { + Err(last_error()) + } else { + Ok(Self { raw }) + } + } + + pub fn from_options(opts: Option<&MlfmmOptions>) -> Result { + let mut ctx = Self::new()?; + ctx.apply_options(opts)?; + Ok(ctx) + } + + pub fn set_sources_linear( + &mut self, + rs_xyz: (&[f64], &[f64], &[f64]), + drs_xyz: (&[f64], &[f64], &[f64]), + currents: &[f64], + eps: &[f64], + ) -> Result<(), String> { + let num_sources = rs_xyz.0.len(); + check_length_3tup_result(num_sources, rs_xyz).map_err(|err| err.to_string())?; + check_length_3tup_result(num_sources, drs_xyz).map_err(|err| err.to_string())?; + if currents.len() != num_sources || eps.len() != num_sources { + return Err("currents and eps must match sources length".to_string()); + } + let ok = unsafe { + ffi::rat_mlfmm_context_set_sources_linear( + self.raw, + rs_xyz.0.as_ptr(), + rs_xyz.1.as_ptr(), + rs_xyz.2.as_ptr(), + drs_xyz.0.as_ptr(), + drs_xyz.1.as_ptr(), + drs_xyz.2.as_ptr(), + currents.as_ptr(), + eps.as_ptr(), + num_sources, + ) + }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_targets(&mut self, rt_xyz: (&[f64], &[f64], &[f64])) -> Result<(), String> { + let num_targets = rt_xyz.0.len(); + check_length_3tup_result(num_targets, rt_xyz).map_err(|err| err.to_string())?; + let ok = unsafe { + ffi::rat_mlfmm_context_set_targets( + self.raw, + rt_xyz.0.as_ptr(), + rt_xyz.1.as_ptr(), + rt_xyz.2.as_ptr(), + num_targets, + ) + }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_van_lanen(&mut self, use_van_lanen: bool) -> Result<(), String> { + let ok = unsafe { ffi::rat_mlfmm_context_set_van_lanen(self.raw, use_van_lanen as i32) }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_num_exp(&mut self, num_exp: i32) -> Result<(), String> { + let ok = unsafe { ffi::rat_mlfmm_context_set_num_exp(self.raw, num_exp) }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_direct_mode(&mut self, mode: ffi::RatMlfmmDirectMode) -> Result<(), String> { + let ok = unsafe { ffi::rat_mlfmm_context_set_direct_mode(self.raw, mode) }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_direct_threshold(&mut self, threshold: f64) -> Result<(), String> { + let ok = unsafe { ffi::rat_mlfmm_context_set_direct_threshold(self.raw, threshold) }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + pub fn set_direct_threshold_count(&mut self, threshold: u64) -> Result<(), String> { + self.set_direct_threshold(threshold as f64) + } + + pub fn compute_ba( + &mut self, + out_b_xyz: (&mut [f64], &mut [f64], &mut [f64]), + out_a_xyz: (&mut [f64], &mut [f64], &mut [f64]), + ) -> Result<(), String> { + let n = out_b_xyz.0.len(); + check_length_3tup_result(n, (&out_b_xyz.0, &out_b_xyz.1, &out_b_xyz.2)) + .map_err(|err| err.to_string())?; + check_length_3tup_result(n, (&out_a_xyz.0, &out_a_xyz.1, &out_a_xyz.2)) + .map_err(|err| err.to_string())?; + let ok = unsafe { + ffi::rat_mlfmm_context_compute_ba( + self.raw, + out_b_xyz.0.as_mut_ptr(), + out_b_xyz.1.as_mut_ptr(), + out_b_xyz.2.as_mut_ptr(), + n, + out_a_xyz.0.as_mut_ptr(), + out_a_xyz.1.as_mut_ptr(), + out_a_xyz.2.as_mut_ptr(), + n, + ) + }; + if ok == 0 { Err(last_error()) } else { Ok(()) } + } + + fn apply_options(&mut self, opts: Option<&MlfmmOptions>) -> Result<(), String> { + let default_opts; + let opts = match opts { + Some(value) => value, + None => { + default_opts = MlfmmOptions::default(); + &default_opts + } + }; + + self.set_van_lanen(opts.use_linear_filament)?; + self.set_direct_mode(ffi::RatMlfmmDirectMode::Threshold)?; + let threshold = opts + .direct_threshold + .unwrap_or(DEFAULT_DIRECT_THRESHOLD) + .max(1); + self.set_direct_threshold_count(threshold)?; + if let Some(order) = opts.order { + self.set_num_exp(order)?; + } + Ok(()) + } +} + +impl Drop for Context { + fn drop(&mut self) { + unsafe { ffi::rat_mlfmm_context_destroy(self.raw) }; + } +} + +#[derive(Clone, Debug)] +pub struct MlfmmOptions { + pub use_linear_filament: bool, + pub direct_threshold: Option, + pub order: Option, +} + +impl Default for MlfmmOptions { + fn default() -> Self { + Self { + use_linear_filament: true, + direct_threshold: None, + order: None, + } + } +} + +const DEFAULT_DIRECT_THRESHOLD: u64 = 10_000_000; + +/// MLFMM calculation of B-field and A-field from linear filament segments. +/// +/// Args: +/// xyzp: observation points (x, y, z), shape (3, N). +/// xyzfil: filament segment start points (x, y, z), shape (3, M). +/// dlxyzfil: segment deltas from start to end (x, y, z), shape (3, M). +/// ifil: filament segment currents, length M. +/// eps: Van Lanen softening parameter, length M. +/// opts: MLFMM options (direct threshold, linear filament kernel, etc.). +/// out_b_xyz: output B-field components (x, y, z), length N each. +/// out_a_xyz: output A-field components (x, y, z), length N each. +/// +/// Returns: +/// Ok(()) on success, Err(String) on failure. +pub fn fields_linear_filament_mlfmm( + xyzp: (&[f64], &[f64], &[f64]), + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + eps: &[f64], + opts: Option<&MlfmmOptions>, + out_b_xyz: (&mut [f64], &mut [f64], &mut [f64]), + out_a_xyz: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), String> { + let mut ctx = Context::from_options(opts)?; + ctx.set_sources_linear(xyzfil, dlxyzfil, ifil, eps)?; + ctx.set_targets(xyzp)?; + ctx.compute_ba(out_b_xyz, out_a_xyz) +} + +fn last_error() -> String { + unsafe { + let err = ffi::rat_mlfmm_last_error(); + if err.is_null() { + return "unknown error".to_string(); + } + CStr::from_ptr(err).to_string_lossy().into_owned() + } +} + +fn check_length_3tup_result(n: usize, tuple: (&[f64], &[f64], &[f64])) -> Result<(), &'static str> { + check_length_3tup!(n, tuple); + Ok(()) +} + +#[cfg(test)] +mod tests { + use super::{MlfmmOptions, fields_linear_filament_mlfmm}; + use crate::physics::linear_filament::flux_density_linear_filament; + use crate::physics::linear_filament::vector_potential_linear_filament; + + #[test] + fn compare_mlfmm_with_linear_filament() { + let xfil = [0.0, 0.5]; + let yfil = [0.0, 0.0]; + let zfil = [0.0, 0.0]; + let dlx = [0.5, 0.5]; + let dly = [0.0, 0.0]; + let dlz = [0.0, 0.0]; + let ifil = [10.0, 10.0]; + + let eps = [1e-9, 1e-9]; + + let targets_x = [0.25, 0.75, 0.5]; + let targets_y = [0.1, 0.2, 0.3]; + let targets_z = [0.0, 0.1, -0.2]; + + let xp = targets_x; + let yp = targets_y; + let zp = targets_z; + + let mut bx = vec![0.0; xp.len()]; + let mut by = vec![0.0; xp.len()]; + let mut bz = vec![0.0; xp.len()]; + let wire_radius = vec![0.0; ifil.len()]; + flux_density_linear_filament( + (&xp, &yp, &zp), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &wire_radius, + (&mut bx, &mut by, &mut bz), + ) + .expect("linear filament calc failed"); + + let opts = MlfmmOptions { + use_linear_filament: false, + direct_threshold: Some(1_000_000), + order: None, + }; + let mut bx_mlfmm = vec![0.0; xp.len()]; + let mut by_mlfmm = vec![0.0; xp.len()]; + let mut bz_mlfmm = vec![0.0; xp.len()]; + let mut ax_mlfmm = vec![0.0; xp.len()]; + let mut ay_mlfmm = vec![0.0; xp.len()]; + let mut az_mlfmm = vec![0.0; xp.len()]; + fields_linear_filament_mlfmm( + (&targets_x, &targets_y, &targets_z), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &eps, + Some(&opts), + (&mut bx_mlfmm, &mut by_mlfmm, &mut bz_mlfmm), + (&mut ax_mlfmm, &mut ay_mlfmm, &mut az_mlfmm), + ) + .expect("mlfmm compute failed"); + + let tol = 1e-4_f64; + for i in 0..3 { + let expect = [bx[i], by[i], bz[i]]; + let got = [bx_mlfmm[i], by_mlfmm[i], bz_mlfmm[i]]; + for j in 0..3 { + let denom = expect[j].abs().max(1.0); + let err = (got[j] - expect[j]).abs() / denom; + assert!( + err <= tol, + "component mismatch at target {i} axis {j}: got {}, expected {}, rel err {}", + got[j], + expect[j], + err + ); + } + } + } + + #[test] + fn compare_mlfmm_with_linear_filament_van_lanen_fmm() { + let xfil = [0.0, 0.5]; + let yfil = [0.0, 0.0]; + let zfil = [0.0, 0.0]; + let dlx = [0.5, 0.5]; + let dly = [0.0, 0.0]; + let dlz = [0.0, 0.0]; + let ifil = [10.0, 10.0]; + + let eps = [1e-3, 1e-3]; + + let targets_x = [0.25, 0.75, 0.5]; + let targets_y = [2.0, -2.5, 3.0]; + let targets_z = [0.5, 1.0, -1.5]; + + let xp = targets_x; + let yp = targets_y; + let zp = targets_z; + + let mut bx = vec![0.0; xp.len()]; + let mut by = vec![0.0; xp.len()]; + let mut bz = vec![0.0; xp.len()]; + let wire_radius = vec![0.0; ifil.len()]; + flux_density_linear_filament( + (&xp, &yp, &zp), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &wire_radius, + (&mut bx, &mut by, &mut bz), + ) + .expect("linear filament calc failed"); + + let opts = MlfmmOptions { + use_linear_filament: true, + direct_threshold: Some(1), + order: None, + }; + let mut bx_mlfmm = vec![0.0; xp.len()]; + let mut by_mlfmm = vec![0.0; xp.len()]; + let mut bz_mlfmm = vec![0.0; xp.len()]; + let mut ax_mlfmm = vec![0.0; xp.len()]; + let mut ay_mlfmm = vec![0.0; xp.len()]; + let mut az_mlfmm = vec![0.0; xp.len()]; + fields_linear_filament_mlfmm( + (&targets_x, &targets_y, &targets_z), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &eps, + Some(&opts), + (&mut bx_mlfmm, &mut by_mlfmm, &mut bz_mlfmm), + (&mut ax_mlfmm, &mut ay_mlfmm, &mut az_mlfmm), + ) + .expect("mlfmm compute failed"); + + let tol = 5e-3_f64; + for i in 0..3 { + let expect = [bx[i], by[i], bz[i]]; + let got = [bx_mlfmm[i], by_mlfmm[i], bz_mlfmm[i]]; + for j in 0..3 { + let denom = expect[j].abs().max(1.0); + let err = (got[j] - expect[j]).abs() / denom; + assert!( + err <= tol, + "component mismatch at target {i} axis {j}: got {}, expected {}, rel err {}", + got[j], + expect[j], + err + ); + } + } + } + + #[test] + fn compare_mlfmm_with_vector_potential() { + let xfil = [0.0, 0.5]; + let yfil = [0.0, 0.0]; + let zfil = [0.0, 0.0]; + let dlx = [0.5, 0.5]; + let dly = [0.0, 0.0]; + let dlz = [0.0, 0.0]; + let ifil = [10.0, 10.0]; + + let eps = [1e-6, 1e-6]; + + let targets_x = [0.25, 0.75, 0.5]; + let targets_y = [0.2, -0.4, 0.3]; + let targets_z = [0.0, 0.1, -0.2]; + + let xp = targets_x; + let yp = targets_y; + let zp = targets_z; + + let mut ax = vec![0.0; xp.len()]; + let mut ay = vec![0.0; xp.len()]; + let mut az = vec![0.0; xp.len()]; + vector_potential_linear_filament( + (&xp, &yp, &zp), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &vec![0.0; ifil.len()], + (&mut ax, &mut ay, &mut az), + ) + .expect("vector potential calc failed"); + + let opts = MlfmmOptions { + use_linear_filament: false, + direct_threshold: Some(1_000_000), + order: None, + }; + let mut bx_mlfmm = vec![0.0; xp.len()]; + let mut by_mlfmm = vec![0.0; xp.len()]; + let mut bz_mlfmm = vec![0.0; xp.len()]; + let mut ax_mlfmm = vec![0.0; xp.len()]; + let mut ay_mlfmm = vec![0.0; xp.len()]; + let mut az_mlfmm = vec![0.0; xp.len()]; + fields_linear_filament_mlfmm( + (&targets_x, &targets_y, &targets_z), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &eps, + Some(&opts), + (&mut bx_mlfmm, &mut by_mlfmm, &mut bz_mlfmm), + (&mut ax_mlfmm, &mut ay_mlfmm, &mut az_mlfmm), + ) + .expect("mlfmm compute failed"); + + let tol = 1e-5_f64; + for i in 0..3 { + let expect = [ax[i], ay[i], az[i]]; + let got = [ax_mlfmm[i], ay_mlfmm[i], az_mlfmm[i]]; + for j in 0..3 { + let denom = expect[j].abs().max(1.0); + let err = (got[j] - expect[j]).abs() / denom; + assert!( + err <= tol, + "component mismatch at target {i} axis {j}: got {}, expected {}, rel err {}", + got[j], + expect[j], + err + ); + } + } + } + + #[test] + fn compare_mlfmm_with_vector_potential_van_lanen_fmm() { + let xfil = [0.0, 0.5]; + let yfil = [0.0, 0.0]; + let zfil = [0.0, 0.0]; + let dlx = [0.5, 0.5]; + let dly = [0.0, 0.0]; + let dlz = [0.0, 0.0]; + let ifil = [10.0, 10.0]; + + let eps = [1e-3, 1e-3]; + + let targets_x = [0.25, 0.75, 0.5]; + let targets_y = [2.0, -2.5, 3.0]; + let targets_z = [0.5, 1.0, -1.5]; + + let xp = targets_x; + let yp = targets_y; + let zp = targets_z; + + let mut ax = vec![0.0; xp.len()]; + let mut ay = vec![0.0; xp.len()]; + let mut az = vec![0.0; xp.len()]; + vector_potential_linear_filament( + (&xp, &yp, &zp), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &vec![0.0; ifil.len()], + (&mut ax, &mut ay, &mut az), + ) + .expect("vector potential calc failed"); + + let opts = MlfmmOptions { + use_linear_filament: true, + direct_threshold: Some(1), + order: None, + }; + let mut bx_mlfmm = vec![0.0; xp.len()]; + let mut by_mlfmm = vec![0.0; xp.len()]; + let mut bz_mlfmm = vec![0.0; xp.len()]; + let mut ax_mlfmm = vec![0.0; xp.len()]; + let mut ay_mlfmm = vec![0.0; xp.len()]; + let mut az_mlfmm = vec![0.0; xp.len()]; + fields_linear_filament_mlfmm( + (&targets_x, &targets_y, &targets_z), + (&xfil, &yfil, &zfil), + (&dlx, &dly, &dlz), + &ifil, + &eps, + Some(&opts), + (&mut bx_mlfmm, &mut by_mlfmm, &mut bz_mlfmm), + (&mut ax_mlfmm, &mut ay_mlfmm, &mut az_mlfmm), + ) + .expect("mlfmm compute failed"); + + let tol = 5e-3_f64; + for i in 0..3 { + let expect = [ax[i], ay[i], az[i]]; + let got = [ax_mlfmm[i], ay_mlfmm[i], az_mlfmm[i]]; + for j in 0..3 { + let denom = expect[j].abs().max(1.0); + let err = (got[j] - expect[j]).abs() / denom; + assert!( + err <= tol, + "component mismatch at target {i} axis {j}: got {}, expected {}, rel err {}", + got[j], + expect[j], + err + ); + } + } + } +} diff --git a/src/physics/biotsavart.rs b/src/physics/biotsavart.rs deleted file mode 100644 index 7aa39d35..00000000 --- a/src/physics/biotsavart.rs +++ /dev/null @@ -1,7 +0,0 @@ -//! Biot-Savart calculations for B-field from current filaments. - -// Preserve references for backwards-compatibility -pub use crate::physics::linear_filament::{ - flux_density_linear_filament as flux_density_biot_savart, - flux_density_linear_filament_par as flux_density_biot_savart_par, -}; diff --git a/src/physics/circular_filament.rs b/src/physics/circular_filament.rs index 438e1231..34d16166 100644 --- a/src/physics/circular_filament.rs +++ b/src/physics/circular_filament.rs @@ -1019,7 +1019,8 @@ mod test { &xyzfil1.2[..n - 1], ), dl1, - &vec![1.0; x.len()][..], + &vec![1.0; n - 1][..], + &vec![0.0; n - 1][..], (&midpoints(&xi), &midpoints(&yi), &midpoints(&zi)), j2, (outxi, outyi, outzi), @@ -1099,6 +1100,7 @@ mod test { ), dlxyzfil, &ifil[..], + &vec![0.0; ifil.len()], (xcontrib, ycontrib, zcontrib), ) .unwrap(); diff --git a/src/physics/linear_filament.rs b/src/physics/linear_filament.rs index 7248a3bf..ed110a4a 100644 --- a/src/physics/linear_filament.rs +++ b/src/physics/linear_filament.rs @@ -7,11 +7,14 @@ use rayon::{ use crate::{ chunksize, - math::{cross3, cross3f, decompose_filament, dot3, dot3f, rss3}, + math::{cross3, dot3, rss3}, }; use crate::{MU0_OVER_4PI, macros::*}; +/// (m) minimum representable nonzero wire thickness. +const MIN_WIRE_THICKNESS: f64 = 1e-10; + /// Estimate the mutual inductance between two piecewise-linear current filaments. /// /// Uses filament midpoints as field source and target. @@ -149,12 +152,14 @@ pub fn inductance_piecewise_linear_filaments( /// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `out`: (T) bx, by, bz at observation points, each length `n` pub fn flux_density_linear_filament_par( xyzp: (&[f64], &[f64], &[f64]), xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result<(), &'static str> { // Chunk inputs @@ -166,7 +171,14 @@ pub fn flux_density_linear_filament_par( (bxc, byc, bzc, xpc, ypc, zpc) .into_par_iter() .try_for_each(|(bx, by, bz, xp, yp, zp)| { - flux_density_linear_filament((xp, yp, zp), xyzfil, dlxyzfil, ifil, (bx, by, bz)) + flux_density_linear_filament( + (xp, yp, zp), + xyzfil, + dlxyzfil, + ifil, + wire_radius, + (bx, by, bz), + ) })?; Ok(()) @@ -183,12 +195,14 @@ pub fn flux_density_linear_filament_par( /// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `out`: (T) bx, by, bz at observation points, each length `n` pub fn flux_density_linear_filament( xyzp: (&[f64], &[f64], &[f64]), xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result<(), &'static str> { // Unpack @@ -203,7 +217,17 @@ pub fn flux_density_linear_filament( let n = xfil.len(); let m = xp.len(); check_length!(m, xp, yp, zp, bx, by, bz); - check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil, ifil); + check_length!( + n, + xfil, + yfil, + zfil, + dlxfil, + dlyfil, + dlzfil, + ifil, + wire_radius + ); // Zero output bx.fill(0.0); @@ -225,7 +249,8 @@ pub fn flux_density_linear_filament( let obs = (xp[j], yp[j], zp[j]); // [m] // Field contributions - let (bxc, byc, bzc) = flux_density_linear_filament_scalar((fil0, fil1, current), obs); + let (bxc, byc, bzc) = + flux_density_linear_filament_scalar((fil0, fil1, current), wire_radius[i], obs); bx[j] += bxc; by[j] += byc; bz[j] += bzc; @@ -235,14 +260,54 @@ pub fn flux_density_linear_filament( Ok(()) } -/// Biot-Savart calculation for B-field contribution from many current filament -/// segments to many observation points. +/// Biot-Savart calculation for B-field contribution one filament +/// to one observation point. /// -/// Uses filament midpoint as field source. +/// Uses the formula for continuous current distribution and finite wire thickness. +/// Inside the wire radius, the field blends linearly to zero at the center. +/// +/// Draws from Griffiths eq'n 5.37 and Zahn eq'n 5.4.17 with inspiration from +/// rat-mlfmm's van Lanen kernel to replace the expensive sine functions with geometric +/// equivalents and to provide handling of the singularity near the filament axis. +/// +/// ```text +/// p (target) +/// * +/// /|\ +/// / | \ +/// ap/ | \bp +/// / ∠a|∠b \ +/// / | \ +/// a-----m-----b -> I +/// | +/// | d_perp (from line to p) +/// | +/// q (closest point on line) +///``` +/// +/// The base formula is +/// +/// $ |B| = \frac{\mu_0 I}{4 \pi r_\perp} (sin(\theta_b) - sin(\theta_a)) $ +/// +/// with the direction determined by $ \hat{dL} \times \hat r_\perp $ with $r_\perp$ defined perpendicular +/// from the axis of the filament to the target point. +/// +/// The otherwise-expensive sine functions are evaluated directly using distance magnitudes. +/// +/// Inside the wire radius, the formula is modified to blend linearly to zero at the wire center. +/// +/// ## References +/// +/// * \[1\] D. J. Griffiths, Introduction to electrodynamics, Fourth edition. Boston: Pearson, 2014. +/// * \[2\] J. van Nugteren and N. Deelen, “rat-mlfmm,” GitLab repository. Accessed: Jan. 16, 2026. [Online]. +/// Available: https://gitlab.com/Project-Rat/rat-mlfmm/-/tree/1e1d387522fafac50c0540af1ebb15d1d506d33d +/// * \[3\] M. Zahn, “5.4: The Vector Potential,” Engineering LibreTexts. Accessed: Jan. 20, 2026. [Online]. +/// Available: https://eng.libretexts.org/Bookshelves/Electrical_Engineering/Electro-Optics/Electromagnetic_Field_Theory%3A_A_Problem_Solving_Approach_(Zahn)/05%3A_The_Magnetic_Field/5.04%3A_The_Vector_Potential /// /// # Arguments /// /// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `wire_radius`: (m) (Half-) thickness of conductor. /// * `xyzp`: (m) Observation point coords /// /// # Returns @@ -250,51 +315,72 @@ pub fn flux_density_linear_filament( /// * `b`: (T) Magnetic flux density (B-field) pub fn flux_density_linear_filament_scalar( xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + wire_radius: f64, xyzobs: (f64, f64, f64), ) -> (f64, f64, f64) { + use crate::math::{PointLineDistance, point_line_distance_with_endpoints}; + // Unpack - let (xyz0, xyz1, ifil) = xyzifil; + let (start, end, ifil) = xyzifil; let (xp, yp, zp) = xyzobs; - // Get filament midpoint and length vector - let ((xmid, ymid, zmid), dl) = decompose_filament(xyz0, xyz1); - - // Get distance from middle of the filament segment to the observation point - let rx: f64 = xp - xmid; // [m] - let ry = yp - ymid; // [m] - let rz = zp - zmid; // [m] - - // Now that we've resolved the part of the calculation that involves a wide dynamic range, - // which drives the need for 64-bit floats to control roundoff error, - // we can switch to 32-bit floats for the majority of the calculation without incurring - // excessive error, before converting back to 64-bit float so that we maintain - // acceptable error during summation downstream. - let (rx, ry, rz) = (rx as f32, ry as f32, rz as f32); - let dl = (dl.0 as f32, dl.1 as f32, dl.2 as f32); - let ifil = ifil as f32; - - // Do 1/r^3 operation with an ordering that improves float error by eliminating - // the actual cube operation and using fused multiply-add to reduce roundoff events, - // then rolling the result into the factor that is constant between all contributions. - let sumsq = dot3f(rx, ry, rz, rx, ry, rz); - let rnorm3_inv = sumsq.powf(-1.5); // [m^-3] - - // This factor is constant across all x, y, and z components - let c = (MU0_OVER_4PI as f32) * ifil * rnorm3_inv; - - // Evaluate the cross products for each axis component - // separately using mul_add which would not be assumed usable - // in a more general implementation. - let (cx, cy, cz) = cross3f(dl.0, dl.1, dl.2, rx, ry, rz); - - // Assemble final B-field components - // and upcast back to 64-bit float so that summation operations - // downstream do not incur excessive roundoff error. - let bx = (c * cx) as f64; // [T] - let by = (c * cy) as f64; - let bz = (c * cz) as f64; - - (bx, by, bz) + // Get perpendicular distance and distance from each endpoint to the target, + // and a fraction between 0 and 1 representing how far the point is from the center of the wire + // to the edge of the wire. + // All 3 distances are clamped to at least the wire radius. + let PointLineDistance { + perp, + dist_a, + dist_b, + frac, + para_a, + para_b, + ab_norm: dlhat, + } = point_line_distance_with_endpoints(start, end, xyzobs, wire_radius); + + // Perpendicular unit vector from the line to the target for field direction. + let ap = (xp - start.0, yp - start.1, zp - start.2); + let ap_para = dot3(ap.0, ap.1, ap.2, dlhat.0, dlhat.1, dlhat.2); + let perp_vec = ( + ap.0 - ap_para * dlhat.0, + ap.1 - ap_para * dlhat.1, + ap.2 - ap_para * dlhat.2, + ); + let perp_inv = if perp > 0.0 { 1.0 / perp } else { 0.0 }; + let perp_hat = ( + perp_vec.0 * perp_inv, + perp_vec.1 * perp_inv, + perp_vec.2 * perp_inv, + ); + + // Sine of the angle formed by the lines from the target to each endpoint + // and the line of the filament. + let sin_theta_a = para_a / dist_a; // (dimensionless) + let sin_theta_b = para_b / dist_b; // (dimensionless) + + // Geometric component of B-field magnitude, + // including linear falloff inside finite-thickness wire. + let geometric_factor = -frac * (sin_theta_b - sin_theta_a); // (dimensionless) + + // This factor is constant across all x, y, and z components. + let c = geometric_factor * MU0_OVER_4PI * ifil / perp; // (A/m) + + // Direction of cross(dL, r), the direction of the field. + let (cx, cy, cz) = cross3( + dlhat.0, dlhat.1, dlhat.2, perp_hat.0, perp_hat.1, perp_hat.2, + ); // (dimensionless) + + // Assemble final B-field components. + let bx = c * cx; // [T] + let by = c * cy; // [T] + let bz = c * cz; // [T] + + // Finally, determine whether we are clipping to zero. + if frac > 1e6 * f64::EPSILON && perp > MIN_WIRE_THICKNESS { + return (bx, by, bz); + } else { + return (0.0, 0.0, 0.0); + } } /// Vector potential calculation for A-field contribution from many current filament @@ -310,12 +396,14 @@ pub fn flux_density_linear_filament_scalar( /// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `out`: (V-s/m) ax, ay, az at observation points, each length `n` pub fn vector_potential_linear_filament_par( xyzp: (&[f64], &[f64], &[f64]), xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result<(), &'static str> { // Chunk inputs @@ -327,7 +415,14 @@ pub fn vector_potential_linear_filament_par( (bxc, byc, bzc, xpc, ypc, zpc) .into_par_iter() .try_for_each(|(bx, by, bz, xp, yp, zp)| { - vector_potential_linear_filament((xp, yp, zp), xyzfil, dlxyzfil, ifil, (bx, by, bz)) + vector_potential_linear_filament( + (xp, yp, zp), + xyzfil, + dlxyzfil, + ifil, + wire_radius, + (bx, by, bz), + ) })?; Ok(()) @@ -344,12 +439,14 @@ pub fn vector_potential_linear_filament_par( /// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `out`: (V-s/m) ax, ay, az at observation points, each length `n` pub fn vector_potential_linear_filament( xyzp: (&[f64], &[f64], &[f64]), xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], out: (&mut [f64], &mut [f64], &mut [f64]), ) -> Result<(), &'static str> { // Unpack @@ -364,7 +461,17 @@ pub fn vector_potential_linear_filament( let n = xfil.len(); let m = xp.len(); check_length!(m, xp, yp, zp, ax, ay, az); - check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil, ifil); + check_length!( + n, + xfil, + yfil, + zfil, + dlxfil, + dlyfil, + dlzfil, + ifil, + wire_radius + ); // Zero output ax.fill(0.0); @@ -387,7 +494,7 @@ pub fn vector_potential_linear_filament( // Field contributions let (axc, ayc, azc) = - vector_potential_linear_filament_scalar((fil0, fil1, current), obs); + vector_potential_linear_filament_scalar((fil0, fil1, current), wire_radius[i], obs); ax[j] += axc; ay[j] += ayc; az[j] += azc; @@ -397,10 +504,31 @@ pub fn vector_potential_linear_filament( Ok(()) } -/// Vector potential (A-field) from a linear current -/// filament segment to an observation point. +/// Vector potential (A-field) from a linear current filament segment to an observation point. /// -/// Uses filament midpoint as field source. +/// Uses the formula for finite segment length and finite wire thickness. +/// +/// The base formula implemented here is: +/// +/// $$ +/// A_z +/// = \frac{\mu_0 I}{4\pi}\int_{-L/2}^{L/2} +/// \frac{dz'}{\sqrt{(z-z')^2+r^2}} +/// = \frac{\mu_0 I}{4\pi}\ln( +/// \frac{ +/// -z + \frac{L}{2} + \sqrt{(z-\frac{L}{2})^2 + r^2} +/// }{ +/// -(z+\frac{L}{2}) + \sqrt{(z+\frac{L}{2})^2 + r^2} +/// } +/// ) +/// $$ +/// +/// This has been manipulated to formulate in terms of components of the distance from the +/// filament endpoints and filament axis to the target point: +/// +/// $$ k1 = -||bp_\parallel|| + ||bp|| $$ +/// $$ k2 = -||ap_\parallel|| + ||ap|| $$ +/// $$ A_\parallel = \frac{\mu_0 I}{4 \pi} \ln (\frac{k1}{k2}) $$ /// /// # Arguments /// @@ -413,28 +541,43 @@ pub fn vector_potential_linear_filament( #[inline] pub fn vector_potential_linear_filament_scalar( xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + wire_radius: f64, xyzobs: (f64, f64, f64), ) -> (f64, f64, f64) { - // Unpack - let (xyz0, xyz1, ifil) = xyzifil; - - // Get filament midpoint and length vector - let ((xmid, ymid, zmid), dl) = decompose_filament(xyz0, xyz1); - - // [m] vector from filament midpoint to obs point - let (rx, ry, rz) = (xyzobs.0 - xmid, xyzobs.1 - ymid, xyzobs.2 - zmid); - let rnorm = rss3(rx, ry, rz); - - // Scale factor shared between all components of A - let c = MU0_OVER_4PI * (ifil / rnorm); - - // Vector potential is linear in the current and segment length - // and goes like 1/R from the segment to the observation point. - let ax = c * dl.0; - let ay = c * dl.1; - let az = c * dl.2; + use crate::math::{PointLineDistance, point_line_distance_with_endpoints}; - (ax, ay, az) + // Unpack + let (start, end, ifil) = xyzifil; + + // Get perpendicular distance and distance from each endpoint to the target, + // and a fraction between 0 and 1 representing how far the point is from the center of the wire + // to the edge of the wire. + // All 3 distances are clamped to at least the wire radius. + let PointLineDistance { + perp, + dist_a, + dist_b, + frac, + para_a, + para_b, + ab_norm: dlhat, + } = point_line_distance_with_endpoints(start, end, xyzobs, wire_radius); + + // Finite segment length log-form with quadratic blend to zero at axis. + let k1 = -para_b + dist_b; + let k2 = -para_a + dist_a; + let frac2 = frac * frac; // Quadratic fall-off (as opposed to linear for B-field) + let a_mag = frac2 * MU0_OVER_4PI * ifil * libm::log((k1 / k2).max(0.0)); + + // Direction is always aligned with the segment. + let (ax, ay, az) = (a_mag * dlhat.0, a_mag * dlhat.1, a_mag * dlhat.2); + + // Finally, determine whether we are clipping to zero. + if frac > 1e6 * f64::EPSILON && perp > MIN_WIRE_THICKNESS { + return (ax, ay, az); + } else { + return (0.0, 0.0, 0.0); + } } /// JxB (Lorentz) body force density (per volume) due to a linear current @@ -445,6 +588,7 @@ pub fn vector_potential_linear_filament_scalar( /// # Arguments /// /// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `wire_radius`: (m) (Half-) thickness of conductor. /// * `xyzobs`: (m) Observation point coords /// * `jobs`: (A/m^2) Current density vector at observation point /// @@ -453,11 +597,12 @@ pub fn vector_potential_linear_filament_scalar( /// * `jxb`: (N/m^3) Body force density pub fn body_force_density_linear_filament_scalar( xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + wire_radius: f64, xyzobs: (f64, f64, f64), jobs: (f64, f64, f64), ) -> (f64, f64, f64) { // Get magnetic flux density at target point - let (bx, by, bz) = flux_density_linear_filament_scalar(xyzifil, xyzobs); // [T] + let (bx, by, bz) = flux_density_linear_filament_scalar(xyzifil, wire_radius, xyzobs); // [T] // Take JxB Lorentz force cross3(jobs.0, jobs.1, jobs.2, bx, by, bz) // [N/m^3] @@ -473,6 +618,7 @@ pub fn body_force_density_linear_filament_scalar( /// * `xyzifil`: (m, m, A) Filament start and end coords and current /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `xyzobs`: (m) Observation point coords /// * `jobs`: (A/m^2) Current density vector at observation point /// * `out`: (N/m^3) Body force density x, y, z components @@ -480,6 +626,7 @@ pub fn body_force_density_linear_filament( xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], xyzobs: (&[f64], &[f64], &[f64]), jobs: (&[f64], &[f64], &[f64]), out: (&mut [f64], &mut [f64], &mut [f64]), @@ -497,7 +644,17 @@ pub fn body_force_density_linear_filament( let n = xfil.len(); let m = xp.len(); - check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil); + check_length!( + n, + xfil, + yfil, + zfil, + dlxfil, + dlyfil, + dlzfil, + ifil, + wire_radius + ); check_length!(m, xp, yp, zp, jx, jy, jz, outx, outy, outz); // Zero output @@ -518,8 +675,12 @@ pub fn body_force_density_linear_filament( let jj = (jx[j], jy[j], jz[j]); // [A/m^2] current density vector at obs point // [V-s/m] vector potential contribution of this filament to this observation point - let (jxbx, jxby, jxbz) = - body_force_density_linear_filament_scalar((fil0, fil1, ifil[i]), obs, jj); + let (jxbx, jxby, jxbz) = body_force_density_linear_filament_scalar( + (fil0, fil1, ifil[i]), + wire_radius[i], + obs, + jj, + ); outx[j] += jxbx; outy[j] += jxby; outz[j] += jxbz; @@ -541,6 +702,7 @@ pub fn body_force_density_linear_filament( /// * `xyzifil`: (m, m, A) Filament start and end coords and current /// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` /// * `ifil`: (A) Filament current, length `m` +/// * `wire_radius`: (m) (Half-) thickness of conductor, length `m` /// * `xyzobs`: (m) Observation point coords /// * `jobs`: (A/m^2) Current density vector at observation point /// * `out`: (N/m^3) Body force density x, y, z components @@ -548,6 +710,7 @@ pub fn body_force_density_linear_filament_par( xyzfil: (&[f64], &[f64], &[f64]), dlxyzfil: (&[f64], &[f64], &[f64]), ifil: &[f64], + wire_radius: &[f64], xyzobs: (&[f64], &[f64], &[f64]), jobs: (&[f64], &[f64], &[f64]), out: (&mut [f64], &mut [f64], &mut [f64]), @@ -566,6 +729,7 @@ pub fn body_force_density_linear_filament_par( xyzfil, dlxyzfil, ifil, + wire_radius, (xp, yp, zp), (jx, jy, jz), (outx, outy, outz), @@ -580,6 +744,9 @@ mod test { use std::f64::consts::PI; use super::*; + use crate::physics::point_source::segment::{ + flux_density_point_segment, vector_potential_point_segment, + }; use crate::testing::*; /// Make sure the forces have the right sign @@ -607,7 +774,8 @@ mod test { body_force_density_linear_filament( (&x[..ndiscr - 1], &y[..ndiscr - 1], &z[..ndiscr - 1]), dl, - &vec![ni; x.len()][..], + &vec![ni; ndiscr - 1][..], + &vec![0.0; ndiscr - 1][..], (&x[..ndiscr - 1], &y[..ndiscr - 1], &z[..ndiscr - 1]), dl, (jxbx, jxby, jxbz), @@ -661,6 +829,7 @@ mod test { (&xi[..ndiscr - 1], &yi[..ndiscr - 1], &zi[..ndiscr - 1]), dli, &vec![ni * nj; xi.len() - 1][..], + &vec![0.0; xi.len() - 1][..], mid, dlj, (jxbx, jxby, jxbz), @@ -679,6 +848,212 @@ mod test { } } + /// Compare single-segment Biot-Savart against discretized point-source segments. + #[test] + fn test_flux_density_against_point_segment_discretization() { + let (rtol, atol) = (1e-6, 1e-12); + + let start = (0.0, 0.0, -0.5); + let end = (0.0, 0.0, 0.5); + let ifil = [1.0]; + + let xfil = [start.0]; + let yfil = [start.1]; + let zfil = [start.2]; + let dlx = [end.0 - start.0]; + let dly = [end.1 - start.1]; + let dlz = [end.2 - start.2]; + let xyzfil = (&xfil[..], &yfil[..], &zfil[..]); + let dlxyz = (&dlx[..], &dly[..], &dlz[..]); + + let ngrid = 100; + let span = 10.0; + let xvals: Vec = (0..ngrid) + .map(|i| -span + (2.0 * span) * (i as f64) / (ngrid as f64 - 1.0)) + .collect(); + let yvals = xvals.clone(); + let zvals = xvals.clone(); + + let total = ngrid * ngrid * ngrid; + let mut xp = Vec::with_capacity(total); + let mut yp = Vec::with_capacity(total); + let mut zp = Vec::with_capacity(total); + for &x in &xvals { + for &y in &yvals { + for &z in &zvals { + xp.push(x); + yp.push(y); + zp.push(z); + } + } + } + let xyzp = (&xp[..], &yp[..], &zp[..]); + + let mut bx = vec![0.0; total]; + let mut by = vec![0.0; total]; + let mut bz = vec![0.0; total]; + flux_density_linear_filament( + xyzp, + xyzfil, + dlxyz, + &ifil, + &[0.0], + (&mut bx, &mut by, &mut bz), + ) + .unwrap(); + + let nseg = 1000; + let dz = (end.2 - start.2) / nseg as f64; + let mut xfil_ps = Vec::with_capacity(nseg); + let mut yfil_ps = Vec::with_capacity(nseg); + let mut zfil_ps = Vec::with_capacity(nseg); + for i in 0..nseg { + xfil_ps.push(start.0); + yfil_ps.push(start.1); + zfil_ps.push(start.2 + dz * i as f64); + } + let dlx = vec![0.0; nseg]; + let dly = vec![0.0; nseg]; + let dlz = vec![dz; nseg]; + let ifil_ps = vec![1.0; nseg]; + + let mut bx_ps = vec![0.0; total]; + let mut by_ps = vec![0.0; total]; + let mut bz_ps = vec![0.0; total]; + flux_density_point_segment( + xyzp, + (&xfil_ps, &yfil_ps, &zfil_ps), + (&dlx, &dly, &dlz), + &ifil_ps, + (&mut bx_ps, &mut by_ps, &mut bz_ps), + ) + .unwrap(); + + for i in 0..xp.len() { + assert!( + approx(bx[i], bx_ps[i], rtol, atol), + "bx is {}, should be {}", + bx[i], + bx_ps[i] + ); + assert!( + approx(by[i], by_ps[i], rtol, atol), + "by is {}, should be {}", + by[i], + by_ps[i] + ); + assert!( + approx(bz[i], bz_ps[i], rtol, atol), + "bz is {}, should be {}", + bz[i], + bz_ps[i] + ); + } + } + + /// Compare single-segment vector potential against discretized point-source segments. + #[test] + fn test_vector_potential_against_point_segment_discretization() { + let (rtol, atol) = (1e-6, 1e-15); + + let start = (0.0, 0.0, -0.5); + let end = (0.0, 0.0, 0.5); + let ifil = [1.0]; + + let xfil = [start.0]; + let yfil = [start.1]; + let zfil = [start.2]; + let dlx = [end.0 - start.0]; + let dly = [end.1 - start.1]; + let dlz = [end.2 - start.2]; + let xyzfil = (&xfil[..], &yfil[..], &zfil[..]); + let dlxyz = (&dlx[..], &dly[..], &dlz[..]); + + let ngrid = 100; + let span = 10.0; + let xvals: Vec = (0..ngrid) + .map(|i| -span + (2.0 * span) * (i as f64) / (ngrid as f64 - 1.0)) + .collect(); + let yvals = xvals.clone(); + let zvals = xvals.clone(); + + let total = ngrid * ngrid * ngrid; + let mut xp = Vec::with_capacity(total); + let mut yp = Vec::with_capacity(total); + let mut zp = Vec::with_capacity(total); + for &x in &xvals { + for &y in &yvals { + for &z in &zvals { + xp.push(x); + yp.push(y); + zp.push(z); + } + } + } + let xyzp = (&xp[..], &yp[..], &zp[..]); + + let mut ax = vec![0.0; total]; + let mut ay = vec![0.0; total]; + let mut az = vec![0.0; total]; + vector_potential_linear_filament( + xyzp, + xyzfil, + dlxyz, + &ifil, + &[0.0], + (&mut ax, &mut ay, &mut az), + ) + .unwrap(); + + let nseg = 1000; + let dz = (end.2 - start.2) / nseg as f64; + let mut xfil_ps = Vec::with_capacity(nseg); + let mut yfil_ps = Vec::with_capacity(nseg); + let mut zfil_ps = Vec::with_capacity(nseg); + for i in 0..nseg { + xfil_ps.push(start.0); + yfil_ps.push(start.1); + zfil_ps.push(start.2 + dz * i as f64); + } + let dlx = vec![0.0; nseg]; + let dly = vec![0.0; nseg]; + let dlz = vec![dz; nseg]; + let ifil_ps = vec![1.0; nseg]; + + let mut ax_ps = vec![0.0; total]; + let mut ay_ps = vec![0.0; total]; + let mut az_ps = vec![0.0; total]; + vector_potential_point_segment( + xyzp, + (&xfil_ps, &yfil_ps, &zfil_ps), + (&dlx, &dly, &dlz), + &ifil_ps, + (&mut ax_ps, &mut ay_ps, &mut az_ps), + ) + .unwrap(); + + for i in 0..xp.len() { + assert!( + approx(ax[i], ax_ps[i], rtol, atol), + "ax is {}, should be {}", + ax[i], + ax_ps[i] + ); + assert!( + approx(ay[i], ay_ps[i], rtol, atol), + "ay is {}, should be {}", + ay[i], + ay_ps[i] + ); + assert!( + approx(az[i], az_ps[i], rtol, atol), + "az is {}, should be {}", + az[i], + az_ps[i] + ); + } + } + /// Check that B = curl(A) #[test] fn test_vector_potential() { @@ -729,6 +1104,7 @@ mod test { (&xyz, &xyz, &xyz), (&dlxyz, &dlxyz, &dlxyz), &[1.0], + &[0.0], (outx, outy, outz), ) .unwrap(); @@ -738,7 +1114,11 @@ mod test { // magnetic flux through a surface bounded by filament 2. The flux through // filament 2 due to 1 ampere of current in filament 1 is the mutual inductance. // (We are stretching the applicability of Stokes' therorem because the filaments - // are not closed loops) + // are not closed loops). + // + // Because inductance_piecewise_linear_filaments uses Neumann's formula, which is + // exactly equivalent to the point-source formulation of the vector potential, + // we expect a small amount of error to the finite-length segment formula here. let a_dot_dl: Vec = (0..NFIL - 1) .map(|i| outx[i] * dlxfil2[i] + outy[i] * dlyfil2[i] + outz[i] * dlzfil2[i]) .collect(); @@ -751,7 +1131,12 @@ mod test { false, ) .unwrap(); - assert!(approx(m, m_from_a, 1e-10, 1e-15)); + assert!( + approx(m, m_from_a, 1e-2, 1e-15), + "m = {:.3e}, m_from_a = {:.3e}", + m, + m_from_a + ); let vp = |x: f64, y: f64, z: f64| { let mut outx = [0.0]; @@ -763,6 +1148,7 @@ mod test { (&xyz, &xyz, &xyz), (&dlxyz, &dlxyz, &dlxyz), &[1.0], + &[0.0], (&mut outx, &mut outy, &mut outz), ) .unwrap(); @@ -771,19 +1157,29 @@ mod test { }; let vals = [ - 0.25, 0.5, 2.5, 10.0, 100.0, 1000.0, -1000.0, -100.0, -10.0, -2.5, -0.5, -0.25, + 0.25, 0.5, 2.1, 10.0, 100.0, 1000.0, -1000.0, -100.0, -10.0, -2.0, -0.5, -0.25, ]; // finite diff delta needs to be small enough to be accurate // but large enough that we can tell the difference between adjacent points // that are very far from the origin - let eps = 1e-7; for x in vals.iter() { for y in vals.iter() { for z in vals.iter() { + // Skip the diagonal, which will land exactly on a filament + // several times, and the field is non-smooth on-axis. + if x == y && x == z { + continue; + } + let x = &(x + 1e-2); // Slightly adjust to avoid nans let y = &(y + 1e-2); let z = &(z - 1e-2); + // Scale tolerance and step size based on distance + let r = rss3(*x, *y, *z); + let atol = 1e-12 / r.max(1.0); // Smaller absolute tolerance as field falls off + let eps = 1e-8 * r; // Larger finite difference delta in far-field for resolution + // Brute-force jac because we're only using it once let mut da = [[0.0; 3]; 3]; // da/dx @@ -824,17 +1220,34 @@ mod test { let mut by = [0.0]; let mut bz = [0.0]; flux_density_linear_filament( - (&[*x], &[*y], &[*z]), - (&xyz, &xyz, &xyz), - (&dlxyz, &dlxyz, &dlxyz), + (&[*x][..], &[*y][..], &[*z][..]), + (&xyz[..], &xyz[..], &xyz[..]), + (&dlxyz[..], &dlxyz[..], &dlxyz[..]), &[1.0], + &[0.0], (&mut bx, &mut by, &mut bz), ) .unwrap(); - assert!(approx(bx[0], ca[0], 1e-6, 1e-15)); - assert!(approx(by[0], ca[1], 1e-6, 1e-15)); - assert!(approx(bz[0], ca[2], 1e-6, 1e-15)); + println!("x,y,z = {:.2},{:.2},{:.2}", x, y, z); + assert!( + approx(bx[0], ca[0], 1e-6, atol), + "bx = {:.6e}, ca[0] = {:.6e}", + bx[0], + ca[0] + ); + assert!( + approx(by[0], ca[1], 1e-6, atol), + "by = {:.6e}, ca[1] = {:.6e}", + by[0], + ca[1] + ); + assert!( + approx(bz[0], ca[2], 1e-6, atol), + "bz = {:.6e}, ca[2] = {:.6e}", + bz[0], + ca[2], + ); } } } @@ -878,8 +1291,25 @@ mod test { let out5 = &mut [5.0; NOBS]; // Flux density - flux_density_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, (out0, out1, out2)).unwrap(); - flux_density_linear_filament_par(xyzp, xyzfil, dlxyzfil, ifil, (out3, out4, out5)).unwrap(); + let wire_radius = vec![0.0; ifil.len()]; + flux_density_linear_filament( + xyzp, + xyzfil, + dlxyzfil, + ifil, + &wire_radius, + (out0, out1, out2), + ) + .unwrap(); + flux_density_linear_filament_par( + xyzp, + xyzfil, + dlxyzfil, + ifil, + &wire_radius, + (out3, out4, out5), + ) + .unwrap(); for i in 0..NOBS { assert_eq!(out0[i], out3[i]); assert_eq!(out1[i], out4[i]); @@ -895,9 +1325,25 @@ mod test { let out5 = &mut [5.0; NOBS]; // Vector potential - vector_potential_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, (out0, out1, out2)).unwrap(); - vector_potential_linear_filament_par(xyzp, xyzfil, dlxyzfil, ifil, (out3, out4, out5)) - .unwrap(); + let wire_radius = vec![0.0; ifil.len()]; + vector_potential_linear_filament( + xyzp, + xyzfil, + dlxyzfil, + ifil, + &wire_radius, + (out0, out1, out2), + ) + .unwrap(); + vector_potential_linear_filament_par( + xyzp, + xyzfil, + dlxyzfil, + ifil, + &wire_radius, + (out3, out4, out5), + ) + .unwrap(); for i in 0..NOBS { assert_eq!(out0[i], out3[i]); assert_eq!(out1[i], out4[i]); diff --git a/src/physics/mesh_filament.rs b/src/physics/mesh_filament.rs index d8b0d5f5..c751a817 100644 --- a/src/physics/mesh_filament.rs +++ b/src/physics/mesh_filament.rs @@ -14,7 +14,7 @@ use crate::{ }; use crate::{ math::{decompose_filament, dot3}, - physics::linear_filament::vector_potential_linear_filament_scalar, + physics::point_source::segment::vector_potential_point_segment_scalar, }; /// Mutual inductance from each edge in mesh 1 to each edge in mesh 2. @@ -89,7 +89,7 @@ where // with unit current in edge 1 in order to extract A/I. let (midpoint2, dl2) = decompose_filament(start2, end2); let (ax_per_amp, ay_per_amp, az_per_amp) = - vector_potential_linear_filament_scalar(xyzifil1, midpoint2); + vector_potential_point_segment_scalar(xyzifil1, midpoint2); // Take M = dot((A/I), dL) let m = dot3(ax_per_amp, ay_per_amp, az_per_amp, dl2.0, dl2.1, dl2.2); diff --git a/src/physics/mod.rs b/src/physics/mod.rs index 68e0f954..a813fe1c 100644 --- a/src/physics/mod.rs +++ b/src/physics/mod.rs @@ -1,5 +1,4 @@ //! Electromagnetics calculations. -pub mod biotsavart; pub mod circular_filament; pub mod gradshafranov; pub mod linear_filament; diff --git a/src/physics/point_source.rs b/src/physics/point_source/dipole.rs similarity index 100% rename from src/physics/point_source.rs rename to src/physics/point_source/dipole.rs diff --git a/src/physics/point_source/mod.rs b/src/physics/point_source/mod.rs new file mode 100644 index 00000000..584fd26f --- /dev/null +++ b/src/physics/point_source/mod.rs @@ -0,0 +1,4 @@ +pub mod dipole; +pub use dipole::*; + +pub mod segment; diff --git a/src/physics/point_source/segment.rs b/src/physics/point_source/segment.rs new file mode 100644 index 00000000..f4e4ffad --- /dev/null +++ b/src/physics/point_source/segment.rs @@ -0,0 +1,784 @@ +//! Magnetics calculations for piecewise-linear current filaments +//! in point-source form. + +use rayon::{ + iter::{IntoParallelIterator, ParallelIterator}, + slice::{ParallelSlice, ParallelSliceMut}, +}; + +use crate::{ + chunksize, + math::{cross3, decompose_filament, dot3, rss3}, +}; + +use crate::{MU0_OVER_4PI, macros::*}; + +/// Biot-Savart calculation for B-field contribution from many current filament +/// segments to many observation points. +/// +/// Uses filament midpoint as field source. +/// +/// This variant of the function is parallelized over chunks of observation points. +/// +/// # Arguments +/// +/// * `xyzp`: (m) Observation point coords, each length `n` +/// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `out`: (T) bx, by, bz at observation points, each length `n` +pub fn flux_density_point_segment_par( + xyzp: (&[f64], &[f64], &[f64]), + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Chunk inputs + let n = chunksize(xyzp.0.len()); + let (xpc, ypc, zpc) = par_chunks_3tup!(xyzp, n); + let (bxc, byc, bzc) = mut_par_chunks_3tup!(out, n); + + // Run calcs + (bxc, byc, bzc, xpc, ypc, zpc) + .into_par_iter() + .try_for_each(|(bx, by, bz, xp, yp, zp)| { + flux_density_point_segment((xp, yp, zp), xyzfil, dlxyzfil, ifil, (bx, by, bz)) + })?; + + Ok(()) +} + +/// Biot-Savart calculation for B-field contribution from many current filament +/// segments to many observation points. +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzp`: (m) Observation point coords, each length `n` +/// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `out`: (T) bx, by, bz at observation points, each length `n` +pub fn flux_density_point_segment( + xyzp: (&[f64], &[f64], &[f64]), + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Unpack + let (xp, yp, zp) = xyzp; + let (xfil, yfil, zfil) = xyzfil; + let (dlxfil, dlyfil, dlzfil) = dlxyzfil; + + let (bx, by, bz) = out; + + // Check lengths; if there is any possibility of mismatch, + // the compiler will bypass vectorization + let n = xfil.len(); + let m = xp.len(); + check_length!(m, xp, yp, zp, bx, by, bz); + check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil, ifil); + + // Zero output + bx.fill(0.0); + by.fill(0.0); + bz.fill(0.0); + + // For each filament, evaluate the contribution to each observation point. + // + // The inner function is inlined, so values that are reused between iterations + // can be pulled to the outer scope by the compiler and do not affect performance. + for i in 0..n { + for j in 0..m { + // Filament + let fil0 = (xfil[i], yfil[i], zfil[i]); // [m] start point + let fil1 = (fil0.0 + dlxfil[i], fil0.1 + dlyfil[i], fil0.2 + dlzfil[i]); // [m] end point + let current = ifil[i]; + + // Observation point + let obs = (xp[j], yp[j], zp[j]); // [m] + + // Field contributions + let (bxc, byc, bzc) = flux_density_point_segment_scalar((fil0, fil1, current), obs); + bx[j] += bxc; + by[j] += byc; + bz[j] += bzc; + } + } + + Ok(()) +} + +/// Biot-Savart calculation for B-field contribution one filament +/// to one observation point. +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `xyzp`: (m) Observation point coords +/// +/// # Returns +/// +/// * `b`: (T) Magnetic flux density (B-field) +pub fn flux_density_point_segment_scalar( + xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + xyzobs: (f64, f64, f64), +) -> (f64, f64, f64) { + // Unpack + let (xyz0, xyz1, ifil) = xyzifil; + let (xp, yp, zp) = xyzobs; + + // Get filament midpoint and length vector + let ((xmid, ymid, zmid), dl) = decompose_filament(xyz0, xyz1); + + // Get distance from middle of the filament segment to the observation point + let rx: f64 = xp - xmid; // [m] + let ry = yp - ymid; // [m] + let rz = zp - zmid; // [m] + + // Now that we've resolved the part of the calculation that involves a wide dynamic range, + // which drives the need for 64-bit floats to control roundoff error, + // we can switch to 32-bit floats for the majority of the calculation without incurring + // excessive error, before converting back to 64-bit float so that we maintain + // acceptable error during summation downstream. + let (rx, ry, rz) = (rx, ry, rz); + let dl = (dl.0, dl.1, dl.2); + let ifil = ifil; + + // Do 1/r^3 operation with an ordering that improves float error by eliminating + // the actual cube operation and using fused multiply-add to reduce roundoff events, + // then rolling the result into the factor that is constant between all contributions. + let sumsq = dot3(rx, ry, rz, rx, ry, rz); + let rnorm3_inv = sumsq.powf(-1.5); // [m^-3] + + // This factor is constant across all x, y, and z components + let c = (MU0_OVER_4PI) * ifil * rnorm3_inv; + + // Evaluate the cross products for each axis component + // separately using mul_add which would not be assumed usable + // in a more general implementation. + let (cx, cy, cz) = cross3(dl.0, dl.1, dl.2, rx, ry, rz); + + // Assemble final B-field components + // and upcast back to 64-bit float so that summation operations + // downstream do not incur excessive roundoff error. + let bx = c * cx; // [T] + let by = c * cy; + let bz = c * cz; + + (bx, by, bz) +} + +/// Vector potential calculation for A-field contribution from many current filament +/// segments to many observation points. +/// +/// Uses filament midpoint as field source. +/// +/// This variant of the function is parallelized over chunks of observation points. +/// +/// # Arguments +/// +/// * `xyzp`: (m) Observation point coords, each length `n` +/// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `out`: (V-s/m) ax, ay, az at observation points, each length `n` +pub fn vector_potential_point_segment_par( + xyzp: (&[f64], &[f64], &[f64]), + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Chunk inputs + let n = chunksize(xyzp.0.len()); + let (xpc, ypc, zpc) = par_chunks_3tup!(xyzp, n); + let (bxc, byc, bzc) = mut_par_chunks_3tup!(out, n); + + // Run calcs + (bxc, byc, bzc, xpc, ypc, zpc) + .into_par_iter() + .try_for_each(|(bx, by, bz, xp, yp, zp)| { + vector_potential_point_segment((xp, yp, zp), xyzfil, dlxyzfil, ifil, (bx, by, bz)) + })?; + + Ok(()) +} + +/// Vector potential calculation for A-field contribution from many current filament +/// segments to many observation points. +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzp`: (m) Observation point coords, each length `n` +/// * `xyzfil`: (m) Filament origin coords (start of segment), each length `m` +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `out`: (V-s/m) ax, ay, az at observation points, each length `n` +pub fn vector_potential_point_segment( + xyzp: (&[f64], &[f64], &[f64]), + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Unpack + let (xp, yp, zp) = xyzp; + let (xfil, yfil, zfil) = xyzfil; + let (dlxfil, dlyfil, dlzfil) = dlxyzfil; + + let (ax, ay, az) = out; + + // Check lengths; if there is any possibility of mismatch, + // the compiler will bypass vectorization + let n = xfil.len(); + let m = xp.len(); + check_length!(m, xp, yp, zp, ax, ay, az); + check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil, ifil); + + // Zero output + ax.fill(0.0); + ay.fill(0.0); + az.fill(0.0); + + // For each filament, evaluate the contribution to each observation point. + // + // The inner function is inlined, so values that are reused between iterations + // can be pulled to the outer scope by the compiler and do not affect performance. + for i in 0..n { + for j in 0..m { + // Filament + let fil0 = (xfil[i], yfil[i], zfil[i]); // [m] start point + let fil1 = (fil0.0 + dlxfil[i], fil0.1 + dlyfil[i], fil0.2 + dlzfil[i]); // [m] end point + let current = ifil[i]; + + // Observation point + let obs = (xp[j], yp[j], zp[j]); // [m] + + // Field contributions + let (axc, ayc, azc) = vector_potential_point_segment_scalar((fil0, fil1, current), obs); + ax[j] += axc; + ay[j] += ayc; + az[j] += azc; + } + } + + Ok(()) +} + +/// Vector potential (A-field) from a linear current +/// filament segment to an observation point. +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `xyzobs`: (m) Observation point coords +/// +/// # Returns +/// +/// * `a`: (V-s/m) Vector potential x, y, z components +#[inline] +pub fn vector_potential_point_segment_scalar( + xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + xyzobs: (f64, f64, f64), +) -> (f64, f64, f64) { + // Unpack + let (xyz0, xyz1, ifil) = xyzifil; + + // Get filament midpoint and length vector + let ((xmid, ymid, zmid), dl) = decompose_filament(xyz0, xyz1); + + // [m] vector from filament midpoint to obs point + let (rx, ry, rz) = (xyzobs.0 - xmid, xyzobs.1 - ymid, xyzobs.2 - zmid); + let rnorm = rss3(rx, ry, rz); + + // Scale factor shared between all components of A + let c = MU0_OVER_4PI * (ifil / rnorm); + + // Vector potential is linear in the current and segment length + // and goes like 1/R from the segment to the observation point. + let ax = c * dl.0; + let ay = c * dl.1; + let az = c * dl.2; + + (ax, ay, az) +} + +/// JxB (Lorentz) body force density (per volume) due to a linear current +/// filament segment at an observation point with some current density (per area). +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `xyzobs`: (m) Observation point coords +/// * `jobs`: (A/m^2) Current density vector at observation point +/// +/// # Returns +/// +/// * `jxb`: (N/m^3) Body force density +pub fn body_force_density_point_segment_scalar( + xyzifil: ((f64, f64, f64), (f64, f64, f64), f64), + xyzobs: (f64, f64, f64), + jobs: (f64, f64, f64), +) -> (f64, f64, f64) { + // Get magnetic flux density at target point + let (bx, by, bz) = flux_density_point_segment_scalar(xyzifil, xyzobs); // [T] + + // Take JxB Lorentz force + cross3(jobs.0, jobs.1, jobs.2, bx, by, bz) // [N/m^3] +} + +/// JxB (Lorentz) body force density (per volume) due to a linear current +/// filament segment at an observation point with some current density (per area). +/// +/// Uses filament midpoint as field source. +/// +/// # Arguments +/// +/// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `xyzobs`: (m) Observation point coords +/// * `jobs`: (A/m^2) Current density vector at observation point +/// * `out`: (N/m^3) Body force density x, y, z components +pub fn body_force_density_point_segment( + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + xyzobs: (&[f64], &[f64], &[f64]), + jobs: (&[f64], &[f64], &[f64]), + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Unpack + let (xp, yp, zp) = xyzobs; + let (jx, jy, jz) = jobs; + let (xfil, yfil, zfil) = xyzfil; + let (dlxfil, dlyfil, dlzfil) = dlxyzfil; + + let (outx, outy, outz) = out; + + // Check lengths; if there is any possibility of mismatch, + // the compiler will bypass vectorization + let n = xfil.len(); + let m = xp.len(); + + check_length!(n, xfil, yfil, zfil, dlxfil, dlyfil, dlzfil); + check_length!(m, xp, yp, zp, jx, jy, jz, outx, outy, outz); + + // Zero output + outx.fill(0.0); + outy.fill(0.0); + outz.fill(0.0); + + // For each filament, evaluate the contribution to each observation point. + // + // The inner function is inlined, so values that are reused between iterations + // can be pulled to the outer scope by the compiler and do not affect performance. + for i in 0..n { + for j in 0..m { + // Geometry + let fil0 = (xfil[i], yfil[i], zfil[i]); // [m] this filament start + let fil1 = (fil0.0 + dlxfil[i], fil0.1 + dlyfil[i], fil0.2 + dlzfil[i]); // [m] this filament end + let obs = (xp[j], yp[j], zp[j]); // [m] this observation point + let jj = (jx[j], jy[j], jz[j]); // [A/m^2] current density vector at obs point + + // [V-s/m] vector potential contribution of this filament to this observation point + let (jxbx, jxby, jxbz) = + body_force_density_point_segment_scalar((fil0, fil1, ifil[i]), obs, jj); + outx[j] += jxbx; + outy[j] += jxby; + outz[j] += jxbz; + } + } + + Ok(()) +} + +/// JxB (Lorentz) body force density (per volume) due to a linear current +/// filament segment at an observation point with some current density (per area). +/// +/// Uses filament midpoint as field source. +/// +/// This variant is parallelized over chunks of observation points. +/// +/// # Arguments +/// +/// * `xyzifil`: (m, m, A) Filament start and end coords and current +/// * `dlxyzfil`: (m) Filament segment length deltas, each length `m` +/// * `ifil`: (A) Filament current, length `m` +/// * `xyzobs`: (m) Observation point coords +/// * `jobs`: (A/m^2) Current density vector at observation point +/// * `out`: (N/m^3) Body force density x, y, z components +pub fn body_force_density_point_segment_par( + xyzfil: (&[f64], &[f64], &[f64]), + dlxyzfil: (&[f64], &[f64], &[f64]), + ifil: &[f64], + xyzobs: (&[f64], &[f64], &[f64]), + jobs: (&[f64], &[f64], &[f64]), + out: (&mut [f64], &mut [f64], &mut [f64]), +) -> Result<(), &'static str> { + // Chunk inputs + let n = chunksize(xyzobs.0.len()); + let (xpc, ypc, zpc) = par_chunks_3tup!(xyzobs, n); + let (jxc, jyc, jzc) = par_chunks_3tup!(jobs, n); + let (outxc, outyc, outzc) = mut_par_chunks_3tup!(out, n); + + // Run calcs + (outxc, outyc, outzc, xpc, ypc, zpc, jxc, jyc, jzc) + .into_par_iter() + .try_for_each(|(outx, outy, outz, xp, yp, zp, jx, jy, jz)| { + body_force_density_point_segment( + xyzfil, + dlxyzfil, + ifil, + (xp, yp, zp), + (jx, jy, jz), + (outx, outy, outz), + ) + })?; + + Ok(()) +} + +#[cfg(test)] +mod test { + use std::f64::consts::PI; + + use super::*; + use crate::physics::linear_filament::inductance_piecewise_linear_filaments; + use crate::testing::*; + + /// Make sure the forces have the right sign + /// and self-forces sum to zero within discretization error + #[test] + fn test_body_force_density() { + let (rtol, atol) = (1e-9, 1e-10); + + let ndiscr = 100; // Discretizations of circular filament into linear filaments + + // Make some circular filaments + let (rfil, zfil, nfil) = example_circular_filaments(); + + // For filament self-field, use the filament roots as the observation points + for i in 0..rfil.len() { + let (ri, zi, ni) = (rfil[i], zfil[i], nfil[i]); + + let (x, y, z) = discretize_circular_filament(ri, zi, ndiscr); + let dl = (&diff(&x)[..], &diff(&y)[..], &diff(&z)[..]); + let (jxbx, jxby, jxbz) = ( + &mut vec![0.0; ndiscr - 1], + &mut vec![0.0; ndiscr - 1], + &mut vec![0.0; ndiscr - 1], + ); + body_force_density_point_segment( + (&x[..ndiscr - 1], &y[..ndiscr - 1], &z[..ndiscr - 1]), + dl, + &vec![ni; x.len()][..], + (&x[..ndiscr - 1], &y[..ndiscr - 1], &z[..ndiscr - 1]), + dl, + (jxbx, jxby, jxbz), + ) + .unwrap(); + + // Make sure the totals sum to zero - a magnet can't put a net force on itself + let jxbx_sum: f64 = jxbx.iter().sum(); + let jxby_sum: f64 = jxby.iter().sum(); + let jxbz_sum: f64 = jxbz.iter().sum(); + assert!(approx(0.0, jxbx_sum, rtol, atol)); + assert!(approx(0.0, jxby_sum, rtol, atol)); + assert!(approx(0.0, jxbz_sum, rtol, atol)); + + // Make sure jxb points outward everywhere + for j in 0..ndiscr - 1 { + let r: (f64, f64, f64) = (x[j], y[j], 0.0); + let rxjxb = cross3(r.0, r.1, r.2, jxbx[j], jxby[j], jxbz[j]); + // Linear filaments aren't perfectly aligned, so we need a slighter wider tolerance here + assert!(approx(0.0, rss3(rxjxb.0, rxjxb.1, rxjxb.2), rtol, 1e-8)); + } + } + + // For filament pairs, make sure axial force is pulling them together + for i in 0..rfil.len() { + let (ri, zif, ni) = (rfil[i], zfil[i], nfil[i]); + let (xi, yi, zi) = discretize_circular_filament(ri, zif, ndiscr); + let dli = (&diff(&xi)[..], &diff(&yi)[..], &diff(&zi)[..]); + + for j in 0..rfil.len() { + // Self-field examined separately + if j == i { + continue; + } + + let (rj, zjf, nj) = (rfil[j], zfil[j], nfil[j]); + let (xj, yj, zj) = discretize_circular_filament(rj, zjf, ndiscr); + let dlj = (&diff(&xj)[..], &diff(&yj)[..], &diff(&zj)[..]); + let mid = ( + &midpoints(&xj)[..], + &midpoints(&yj)[..], + &midpoints(&zj)[..], + ); + + let (jxbx, jxby, jxbz) = ( + &mut vec![0.0; ndiscr - 1], + &mut vec![0.0; ndiscr - 1], + &mut vec![0.0; ndiscr - 1], + ); + body_force_density_point_segment( + (&xi[..ndiscr - 1], &yi[..ndiscr - 1], &zi[..ndiscr - 1]), + dli, + &vec![ni * nj; xi.len() - 1][..], + mid, + dlj, + (jxbx, jxby, jxbz), + ) + .unwrap(); + + // Expect attracting force from j toward i, + // and no centering force because the loops are coaxial + let jxbx_sum: f64 = jxbx.iter().sum(); + let jxby_sum: f64 = jxby.iter().sum(); + let jxbz_sum: f64 = jxbz.iter().sum(); + assert!(approx(0.0, jxbx_sum, rtol, atol)); + assert!(approx(0.0, jxby_sum, rtol, atol)); + assert!(jxbz_sum.signum() == (zif - zjf).signum()); + } + } + } + + /// Check that B = curl(A) + #[test] + fn test_vector_potential() { + // One super basic filament as the source + let xyz = [0.0]; + let dlxyz = [1.0]; + + // Build a second scattering of filament locations as the target + const NFIL: usize = 10; + let xfil2: Vec = (0..NFIL).map(|i| (i as f64).sin() + PI).collect(); + let yfil2: Vec = (0..NFIL).map(|i| (i as f64).cos() - PI).collect(); + let zfil2: Vec = (0..NFIL) + .map(|i| (i as f64) - (NFIL as f64) / 2.0 + PI) + .collect(); + let xyzfil2 = ( + &xfil2[..=NFIL - 2], + &yfil2[..=NFIL - 2], + &zfil2[..=NFIL - 2], + ); + + let dlxfil2: Vec = (0..=NFIL - 2).map(|i| xfil2[i + 1] - xfil2[i]).collect(); + let dlyfil2: Vec = (0..=NFIL - 2).map(|i| yfil2[i + 1] - yfil2[i]).collect(); + let dlzfil2: Vec = (0..=NFIL - 2).map(|i| zfil2[i + 1] - zfil2[i]).collect(); + let dlxyzfil2 = (&dlxfil2[..], &dlyfil2[..], &dlzfil2[..]); + + let xmid2: Vec = xfil2 + .iter() + .zip(dlxfil2.iter()) + .map(|(x, dx)| x + dx / 2.0) + .collect(); + let ymid2: Vec = yfil2 + .iter() + .zip(dlyfil2.iter()) + .map(|(x, dx)| x + dx / 2.0) + .collect(); + let zmid2: Vec = zfil2 + .iter() + .zip(dlzfil2.iter()) + .map(|(x, dx)| x + dx / 2.0) + .collect(); + + // Check against Neumann's formula for mutual inductance + let outx = &mut [0.0; NFIL - 1]; + let outy = &mut [0.0; NFIL - 1]; + let outz = &mut [0.0; NFIL - 1]; + vector_potential_point_segment( + (&xmid2, &ymid2, &zmid2), + (&xyz, &xyz, &xyz), + (&dlxyz, &dlxyz, &dlxyz), + &[1.0], + (outx, outy, outz), + ) + .unwrap(); + // Here the mutual inductance of the two filaments is calculated from the + // vector potential at filament 2 due to 1 ampere of current flowing in filament 1. + // By Stokes' theorem, the line integral of A over filament 2 is equal to the + // magnetic flux through a surface bounded by filament 2. The flux through + // filament 2 due to 1 ampere of current in filament 1 is the mutual inductance. + // (We are stretching the applicability of Stokes' therorem because the filaments + // are not closed loops) + let a_dot_dl: Vec = (0..NFIL - 1) + .map(|i| outx[i] * dlxfil2[i] + outy[i] * dlyfil2[i] + outz[i] * dlzfil2[i]) + .collect(); + let m_from_a = a_dot_dl.iter().sum(); + let m = inductance_piecewise_linear_filaments( + (&xyz, &xyz, &xyz), + (&dlxyz, &dlxyz, &dlxyz), + xyzfil2, + dlxyzfil2, + false, + ) + .unwrap(); + assert!(approx(m, m_from_a, 1e-10, 1e-15)); + + let vp = |x: f64, y: f64, z: f64| { + let mut outx = [0.0]; + let mut outy = [0.0]; + let mut outz = [0.0]; + + vector_potential_point_segment( + (&[x], &[y], &[z]), + (&xyz, &xyz, &xyz), + (&dlxyz, &dlxyz, &dlxyz), + &[1.0], + (&mut outx, &mut outy, &mut outz), + ) + .unwrap(); + + (outx[0], outy[0], outz[0]) + }; + + let vals = [ + 0.25, 0.5, 2.5, 10.0, 100.0, 1000.0, -1000.0, -100.0, -10.0, -2.5, -0.5, -0.25, + ]; + // finite diff delta needs to be small enough to be accurate + // but large enough that we can tell the difference between adjacent points + // that are very far from the origin + let eps = 1e-7; + for x in vals.iter() { + for y in vals.iter() { + for z in vals.iter() { + let x = &(x + 1e-2); // Slightly adjust to avoid nans + let y = &(y + 1e-2); + let z = &(z - 1e-2); + + // Brute-force jac because we're only using it once + let mut da = [[0.0; 3]; 3]; + // da/dx + let (ax0, ay0, az0) = vp(*x - eps, *y, *z); + let (ax1, ay1, az1) = vp(*x + eps, *y, *z); + da[0][0] = (ax1 - ax0) / (2.0 * eps); + da[0][1] = (ay1 - ay0) / (2.0 * eps); + da[0][2] = (az1 - az0) / (2.0 * eps); + + // da/dy + let (ax0, ay0, az0) = vp(*x, *y - eps, *z); + let (ax1, ay1, az1) = vp(*x, *y + eps, *z); + da[1][0] = (ax1 - ax0) / (2.0 * eps); + da[1][1] = (ay1 - ay0) / (2.0 * eps); + da[1][2] = (az1 - az0) / (2.0 * eps); + + // da/dz + let (ax0, ay0, az0) = vp(*x, *y, *z - eps); + let (ax1, ay1, az1) = vp(*x, *y, *z + eps); + da[2][0] = (ax1 - ax0) / (2.0 * eps); + da[2][1] = (ay1 - ay0) / (2.0 * eps); + da[2][2] = (az1 - az0) / (2.0 * eps); + + // B = curl(A) + let daz_dy = da[1][2]; + let day_dz = da[2][1]; + + let daz_dx = da[0][2]; + let dax_dz = da[2][0]; + + let day_dx = da[0][1]; + let dax_dy = da[1][0]; + + let ca = [daz_dy - day_dz, dax_dz - daz_dx, day_dx - dax_dy]; + + // B via biot-savart + let mut bx = [0.0]; + let mut by = [0.0]; + let mut bz = [0.0]; + flux_density_point_segment( + (&[*x], &[*y], &[*z]), + (&xyz, &xyz, &xyz), + (&dlxyz, &dlxyz, &dlxyz), + &[1.0], + (&mut bx, &mut by, &mut bz), + ) + .unwrap(); + + assert!(approx(bx[0], ca[0], 1e-6, 1e-15)); + assert!(approx(by[0], ca[1], 1e-6, 1e-15)); + assert!(approx(bz[0], ca[2], 1e-6, 1e-15)); + } + } + } + } + + /// Check that parallel variants of functions produce the same result as serial. + /// This also incidentally tests defensive zeroing of input slices. + #[test] + fn test_serial_vs_parallel() { + const NFIL: usize = 10; + const NOBS: usize = 100; + + // Build a scattering of filament locations + let xfil: Vec = (0..NFIL).map(|i| (i as f64).sin()).collect(); + let yfil: Vec = (0..NFIL).map(|i| (i as f64).cos()).collect(); + let zfil: Vec = (0..NFIL) + .map(|i| (i as f64) - (NFIL as f64) / 2.0) + .collect(); + let xyzfil = (&xfil[..=NFIL - 2], &yfil[..=NFIL - 2], &zfil[..=NFIL - 2]); + + let dlxfil: Vec = (0..=NFIL - 2).map(|i| xfil[i + 1] - xfil[i]).collect(); + let dlyfil: Vec = (0..=NFIL - 2).map(|i| yfil[i + 1] - yfil[i]).collect(); + let dlzfil: Vec = (0..=NFIL - 2).map(|i| zfil[i + 1] - zfil[i]).collect(); + let dlxyzfil = (&dlxfil[..], &dlyfil[..], &dlzfil[..]); + + let ifil: &[f64] = &(0..NFIL - 1).map(|i| i as f64).collect::>()[..]; + + // Build a scattering of observation locations + let xp: Vec = (0..NOBS).map(|i| 2.0 * (i as f64).sin() + 2.1).collect(); + let yp: Vec = (0..NOBS).map(|i| 4.0 * (2.0 * i as f64).cos()).collect(); + let zp: Vec = (0..NOBS).map(|i| (0.1 * i as f64).exp()).collect(); + let xyzp = (&xp[..], &yp[..], &zp[..]); + + // Some output storage + // Initialize with different values for each buffer to test zeroing + let out0 = &mut [0.0; NOBS]; + let out1 = &mut [1.0; NOBS]; + let out2 = &mut [2.0; NOBS]; + let out3 = &mut [3.0; NOBS]; + let out4 = &mut [4.0; NOBS]; + let out5 = &mut [5.0; NOBS]; + + // Flux density + flux_density_point_segment(xyzp, xyzfil, dlxyzfil, ifil, (out0, out1, out2)).unwrap(); + flux_density_point_segment_par(xyzp, xyzfil, dlxyzfil, ifil, (out3, out4, out5)).unwrap(); + for i in 0..NOBS { + assert_eq!(out0[i], out3[i]); + assert_eq!(out1[i], out4[i]); + assert_eq!(out2[i], out5[i]); + } + + // Reinit to test zeroing + let out0 = &mut [0.0; NOBS]; + let out1 = &mut [1.0; NOBS]; + let out2 = &mut [2.0; NOBS]; + let out3 = &mut [3.0; NOBS]; + let out4 = &mut [4.0; NOBS]; + let out5 = &mut [5.0; NOBS]; + + // Vector potential + vector_potential_point_segment(xyzp, xyzfil, dlxyzfil, ifil, (out0, out1, out2)).unwrap(); + vector_potential_point_segment_par(xyzp, xyzfil, dlxyzfil, ifil, (out3, out4, out5)) + .unwrap(); + for i in 0..NOBS { + assert_eq!(out0[i], out3[i]); + assert_eq!(out1[i], out4[i]); + assert_eq!(out2[i], out5[i]); + } + } +} diff --git a/src/python.rs b/src/python.rs index 1855c768..5d4a06da 100644 --- a/src/python.rs +++ b/src/python.rs @@ -4,6 +4,8 @@ use pyo3::exceptions; use pyo3::prelude::*; use std::fmt::Debug; +#[cfg(feature = "rat-mlfmm")] +use crate::mlfmm::MlfmmOptions; use crate::{math, mesh, physics}; /// Errors from mismatch between python and rust @@ -209,7 +211,7 @@ fn flux_density_circular_filament( } /// Python bindings for cfsemrs::physics::linear_filament::flux_density_linear_filament -#[pyfunction] +#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par=true))] fn flux_density_linear_filament( xyzp: ( PyReadonlyArray1, @@ -226,7 +228,8 @@ fn flux_density_linear_filament( PyReadonlyArray1, PyReadonlyArray1, ), // [m] Filament length delta - ifil: PyReadonlyArray1, // [A] filament current + ifil: PyReadonlyArray1, // [A] filament current + wire_radius: PyReadonlyArray1, // [m] filament radius par: bool, ) -> PyResult<(Py>, Py>, Py>)> { // Get references to contiguous data as slice @@ -235,6 +238,7 @@ fn flux_density_linear_filament( _3tup_slice_ro!(xyzfil); _3tup_slice_ro!(dlxyzfil); let ifil = ifil.as_slice()?; + let wire_radius = wire_radius.as_slice()?; // Do calculations let n = xyzp.0.len(); @@ -244,6 +248,57 @@ fn flux_density_linear_filament( true => physics::linear_filament::flux_density_linear_filament_par, false => physics::linear_filament::flux_density_linear_filament, }; + match func( + xyzp, + xyzfil, + dlxyzfil, + ifil, + wire_radius, + (&mut bx, &mut by, &mut bz), + ) { + Ok(x) => x, + Err(x) => { + let err: PyErr = PyInteropError::DimensionalityError { msg: x.to_string() }.into(); + return Err(err); + } + }; + + _3tup_ret!((bx, f64), (by, f64), (bz, f64)) +} + +/// Python bindings for cfsemrs::physics::point_source::segment::flux_density_point_segment +#[pyfunction] +fn flux_density_point_segment( + xyzp: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Test point coords + xyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament origin coords (start of segment) + dlxyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament length delta + ifil: PyReadonlyArray1, // [A] filament current + par: bool, +) -> PyResult<(Py>, Py>, Py>)> { + _3tup_slice_ro!(xyzp); + _3tup_slice_ro!(xyzfil); + _3tup_slice_ro!(dlxyzfil); + let ifil = ifil.as_slice()?; + + let n = xyzp.0.len(); + let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + + let func = match par { + true => physics::point_source::segment::flux_density_point_segment_par, + false => physics::point_source::segment::flux_density_point_segment, + }; match func(xyzp, xyzfil, dlxyzfil, ifil, (&mut bx, &mut by, &mut bz)) { Ok(x) => x, Err(x) => { @@ -255,8 +310,85 @@ fn flux_density_linear_filament( _3tup_ret!((bx, f64), (by, f64), (bz, f64)) } +/// Python bindings for cfsemrs::mlfmm::fields_linear_filament_mlfmm +#[cfg(feature = "rat-mlfmm")] +#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, eps, use_linear_filament = true, direct_threshold = 10000000, order = None))] +fn fields_linear_filament_mlfmm( + xyzp: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Target coords + xyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament segment start coords + dlxyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament delta from start to end + ifil: PyReadonlyArray1, // [A] filament current + eps: PyReadonlyArray1, // [m] van Lanen softening parameter + use_linear_filament: bool, + direct_threshold: u64, + order: Option, +) -> PyResult<( + (Py>, Py>, Py>), + (Py>, Py>, Py>), +)> { + _3tup_slice_ro!(xyzp); + _3tup_slice_ro!(xyzfil); + _3tup_slice_ro!(dlxyzfil); + let ifil = ifil.as_slice()?; + let eps = eps.as_slice()?; + + let mut opts = MlfmmOptions::default(); + opts.use_linear_filament = use_linear_filament; + if direct_threshold > 0 { + opts.direct_threshold = Some(direct_threshold); + } + if let Some(order) = order { + if order > 0 { + opts.order = Some(order); + } + } + + let n = xyzp.0.len(); + let (mut bx, mut by, mut bz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + let (mut ax, mut ay, mut az) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + + match crate::mlfmm::fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + eps, + Some(&opts), + (&mut bx, &mut by, &mut bz), + (&mut ax, &mut ay, &mut az), + ) { + Ok(x) => x, + Err(x) => { + let err: PyErr = PyInteropError::DimensionalityError { msg: x.to_string() }.into(); + return Err(err); + } + }; + + Python::attach(|py| { + let bx: Py> = PyArray1::from_vec(py, bx).unbind(); + let by: Py> = PyArray1::from_vec(py, by).unbind(); + let bz: Py> = PyArray1::from_vec(py, bz).unbind(); + let ax: Py> = PyArray1::from_vec(py, ax).unbind(); + let ay: Py> = PyArray1::from_vec(py, ay).unbind(); + let az: Py> = PyArray1::from_vec(py, az).unbind(); + Ok(((bx, by, bz), (ax, ay, az))) + }) +} + /// Python bindings for cfsemrs::physics::linear_filament::vector_potential_linear_filament -#[pyfunction] +#[pyfunction(signature = (xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par=true))] fn vector_potential_linear_filament( xyzp: ( PyReadonlyArray1, @@ -273,7 +405,8 @@ fn vector_potential_linear_filament( PyReadonlyArray1, PyReadonlyArray1, ), // [m] Filament length delta - ifil: PyReadonlyArray1, // [A] filament current + ifil: PyReadonlyArray1, // [A] filament current + wire_radius: PyReadonlyArray1, // [m] filament radius par: bool, ) -> PyResult<(Py>, Py>, Py>)> { // Get references to contiguous data as slice @@ -282,6 +415,7 @@ fn vector_potential_linear_filament( _3tup_slice_ro!(xyzfil); _3tup_slice_ro!(dlxyzfil); let ifil = ifil.as_slice()?; + let wire_radius = wire_radius.as_slice()?; // Do calculations let n = xyzp.0.len(); @@ -291,6 +425,57 @@ fn vector_potential_linear_filament( true => physics::linear_filament::vector_potential_linear_filament_par, false => physics::linear_filament::vector_potential_linear_filament, }; + match func( + xyzp, + xyzfil, + dlxyzfil, + ifil, + wire_radius, + (&mut outx, &mut outy, &mut outz), + ) { + Ok(x) => x, + Err(x) => { + let err: PyErr = PyInteropError::DimensionalityError { msg: x.to_string() }.into(); + return Err(err); + } + }; + + _3tup_ret!((outx, f64), (outy, f64), (outz, f64)) +} + +/// Python bindings for cfsemrs::physics::point_source::segment::vector_potential_point_segment +#[pyfunction] +fn vector_potential_point_segment( + xyzp: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Test point coords + xyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament origin coords (start of segment) + dlxyzfil: ( + PyReadonlyArray1, + PyReadonlyArray1, + PyReadonlyArray1, + ), // [m] Filament length delta + ifil: PyReadonlyArray1, // [A] filament current + par: bool, +) -> PyResult<(Py>, Py>, Py>)> { + _3tup_slice_ro!(xyzp); + _3tup_slice_ro!(xyzfil); + _3tup_slice_ro!(dlxyzfil); + let ifil = ifil.as_slice()?; + + let n = xyzp.0.len(); + let (mut outx, mut outy, mut outz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); + + let func = match par { + true => physics::point_source::segment::vector_potential_point_segment_par, + false => physics::point_source::segment::vector_potential_point_segment, + }; match func( xyzp, xyzfil, @@ -644,6 +829,7 @@ fn body_force_density_circular_filament_cartesian( /// Python bindings for cfsemrs::physics::body_force_density_linear_filament #[pyfunction] +#[pyo3(signature = (xyzfil, dlxyzfil, ifil, obs, j, wire_radius, par=true))] fn body_force_density_linear_filament( xyzfil: ( PyReadonlyArray1, @@ -666,6 +852,7 @@ fn body_force_density_linear_filament( PyReadonlyArray1, PyReadonlyArray1, ), // [A/m^2] current density at observation points + wire_radius: PyReadonlyArray1, // [m] filament radius par: bool, ) -> PyResult<(Py>, Py>, Py>)> { // Get references to contiguous data as slice @@ -675,6 +862,7 @@ fn body_force_density_linear_filament( let ifil = ifil.as_slice()?; _3tup_slice_ro!(obs); _3tup_slice_ro!(j); + let wire_radius = wire_radius.as_slice()?; // Select variant let func = match par { @@ -687,7 +875,7 @@ fn body_force_density_linear_filament( let (mut outx, mut outy, mut outz) = (vec![0.0; n], vec![0.0; n], vec![0.0; n]); let out = (&mut outx[..], &mut outy[..], &mut outz[..]); - match func(xyzfil, dlxyzfil, ifil, obs, j, out) { + match func(xyzfil, dlxyzfil, ifil, wire_radius, obs, j, out) { Ok(_) => (), Err(x) => { let err: PyErr = PyInteropError::DimensionalityError { msg: x.to_string() }.into(); @@ -726,10 +914,14 @@ fn _cfsem<'py>(_py: Python, m: Bound<'py, PyModule>) -> PyResult<()> { // Linear filaments m.add_function(wrap_pyfunction!(flux_density_linear_filament, m.clone())?)?; + m.add_function(wrap_pyfunction!(flux_density_point_segment, m.clone())?)?; m.add_function(wrap_pyfunction!( vector_potential_linear_filament, m.clone() )?)?; + m.add_function(wrap_pyfunction!(vector_potential_point_segment, m.clone())?)?; + #[cfg(feature = "rat-mlfmm")] + m.add_function(wrap_pyfunction!(fields_linear_filament_mlfmm, m.clone())?)?; m.add_function(wrap_pyfunction!( inductance_piecewise_linear_filaments, m.clone() diff --git a/test/test_electromagnetics.py b/test/test_electromagnetics.py index e1459bff..76564d0b 100644 --- a/test/test_electromagnetics.py +++ b/test/test_electromagnetics.py @@ -29,11 +29,20 @@ def test_body_force_density(r, z, par): ifil = np.ones_like(xyzfil[0]) jxbx, jxby, jxbz = cfsem.body_force_density_circular_filament_cartesian([1.0], [r], [z], obs, j, par) - jxbx1, jxby1, jxbz1 = cfsem.body_force_density_linear_filament(xyzfil, dlxyzfil, ifil, obs, j, par) + wire_radius = np.zeros_like(ifil) + jxbx1, jxby1, jxbz1 = cfsem.body_force_density_linear_filament( + xyzfil, dlxyzfil, ifil, obs, j, wire_radius, par=par + ) + jxbx2, jxby2, jxbz2 = cfsem.body_force_density_linear_filament( + xyzfil, dlxyzfil, ifil, obs, j, 0.0, par=par + ) assert np.allclose(jxbx, jxbx1, rtol=1e-2, atol=1e-9) assert np.allclose(jxby, jxby1, rtol=1e-2, atol=1e-9) assert np.allclose(jxbz, jxbz1, rtol=1e-2, atol=1e-9) + assert np.allclose(jxbx1, jxbx2, rtol=1e-12, atol=1e-12) + assert np.allclose(jxby1, jxby2, rtol=1e-12, atol=1e-12) + assert np.allclose(jxbz1, jxbz2, rtol=1e-12, atol=1e-12) @mark.parametrize("r", [0.775 * 2, np.pi]) @@ -259,8 +268,14 @@ def test_biot_savart_against_flux_density_ideal_solenoid(r, par): xyzfil = (x1[:-1], y1[:-1], z1[:-1]) # Get B-field at the origin zero = np.array([0.0]) - bx, _by, _bz = cfsem.flux_density_biot_savart( - xyzp=(zero, zero, zero), xyzfil=xyzfil, dlxyzfil=dlxyzfil, ifil=ifil, par=par + wire_radius = np.zeros_like(ifil) + bx, _by, _bz = cfsem.flux_density_linear_filament( + xyzp=(zero, zero, zero), + xyzfil=xyzfil, + dlxyzfil=dlxyzfil, + ifil=ifil, + wire_radius=wire_radius, + par=par, ) b_bs = bx[0] # [T] First and only element on the axis of the solenoid @@ -299,7 +314,10 @@ def test_biot_savart_against_flux_density_circular_filament(r, z, par): xyzfil = (xfils[1:], yfils[1:], zfils[1:]) dlxyzfil = (xfils[1:] - xfils[:-1], yfils[1:] - yfils[:-1], zfils[1:] - zfils[:-1]) ifil = np.ones_like(xfils[1:]) - Br_bs, By_bs, Bz_bs = cfsem.flux_density_biot_savart(xyzp, xyzfil, dlxyzfil, ifil, par) # [T] + wire_radius = np.zeros_like(ifil) + Br_bs, By_bs, Bz_bs = cfsem.flux_density_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par + ) # [T] assert np.allclose(Br_circular, Br_bs, rtol=1e-6, atol=1e-7) # Should match circular calc assert np.allclose(Bz_circular, Bz_bs, rtol=1e-6, atol=1e-7) # ... @@ -756,7 +774,10 @@ def test_vector_potential_linear_against_circular_filament(r, z, par): xyzfil = (xfils[1:], yfils[1:], zfils[1:]) dlxyzfil = (xfils[1:] - xfils[:-1], yfils[1:] - yfils[:-1], zfils[1:] - zfils[:-1]) ifil = np.ones_like(xfils[1:]) - ax, ay, az = cfsem.vector_potential_linear_filament(xyzp, xyzfil, dlxyzfil, ifil, par) # [V-s/m] + wire_radius = np.zeros_like(ifil) + ax, ay, az = cfsem.vector_potential_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par + ) # [V-s/m] assert np.allclose(a_phi, ay, rtol=1e-12, atol=1e-12) # Should match circular calc assert np.allclose(az, np.zeros_like(az), atol=1e-9) # Should sum to zero everywhere @@ -874,4 +895,4 @@ def test_inductance_matrix_axisymmetric_coaxial_rectangular_coils(): td = [1.0, 1.0], nr = [10, 10], nz = [10, 10], - ) \ No newline at end of file + ) diff --git a/test/test_mlfmm.py b/test/test_mlfmm.py new file mode 100644 index 00000000..5e7ed60e --- /dev/null +++ b/test/test_mlfmm.py @@ -0,0 +1,119 @@ +import numpy as np +import pytest +from pytest import mark + +import cfsem + + +def _build_helix_sources(): + n_path = int(1e4) + z = np.linspace(0.0, 1.0, n_path) + path = (np.zeros_like(z), np.zeros_like(z), z) + helix = cfsem.filament_helix_path( + path=path, + helix_start_offset=(5.0, 0.0, 0.0), + twist_pitch=1.0, + angle_offset=0.0, + ) + helix = np.array(helix) + xyzfil = helix[:, :-1] + dlxyzfil = helix[:, 1:] - helix[:, :-1] + ifil = np.ones(dlxyzfil.shape[1]) + eps = np.full_like(ifil, 1e-6) + return xyzfil, dlxyzfil, ifil, eps + + +def _cube_targets(center, half_len, points_per_axis=11): + grid = np.linspace(-half_len, half_len, points_per_axis) + xg, yg, zg = np.meshgrid(grid, grid, grid, indexing="ij") + x = xg.ravel() + center[0] + y = yg.ravel() + center[1] + z = zg.ravel() + center[2] + return x, y, z + + +@mark.parametrize("half_len", [5.025, 100.0]) +@mark.parametrize("direct_threshold", [0, int(1e3), int(1e12)]) +@mark.parametrize("use_linear_filament", [False, True]) +def test_mlfmm_fields_against_direct(half_len, direct_threshold, use_linear_filament): + try: + fields_linear_filament_mlfmm = cfsem.fields_linear_filament_mlfmm + except AttributeError: + pytest.skip("rat-mlfmm feature is not enabled in this build") + + xyzfil, dlxyzfil, ifil, eps = _build_helix_sources() + center = (0.0, 0.0, 0.5) + xyzp = _cube_targets(center=center, half_len=half_len, points_per_axis=11) + + wire_radius = np.zeros_like(ifil) + b_direct = cfsem.flux_density_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par=False + ) + a_direct = cfsem.vector_potential_linear_filament( + xyzp, xyzfil, dlxyzfil, ifil, wire_radius, par=False + ) + + b_mlfmm, a_mlfmm = fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + eps, + use_linear_filament=use_linear_filament, + direct_threshold=direct_threshold, + order=12, + ) + + for got, exp in zip(b_mlfmm, b_direct, strict=True): + assert np.allclose(got, exp, rtol=1e-6, atol=1e-10) + for got, exp in zip(a_mlfmm, a_direct, strict=True): + assert np.allclose(got, exp, rtol=1e-6, atol=1e-10) + + +def test_mlfmm_linear_filament_single_segment_matches_subdivided_direct(): + try: + fields_linear_filament_mlfmm = cfsem.fields_linear_filament_mlfmm + except AttributeError: + pytest.skip("rat-mlfmm feature is not enabled in this build") + + xyzfil = (np.array([0.0]), np.array([0.0]), np.array([0.0])) + dlxyzfil = (np.array([1.0]), np.array([0.0]), np.array([0.0])) + ifil = np.array([1.0]) + eps = np.array([1e-3]) + + t = np.linspace(0.0, 1.0, 101) + x0 = t[:-1] + dx = np.diff(t) + xyzfil_sub = (x0, np.zeros_like(x0), np.zeros_like(x0)) + dlxyzfil_sub = (dx, np.zeros_like(dx), np.zeros_like(dx)) + ifil_sub = np.full_like(dx, ifil[0]) + + xyzp = ( + np.array([0.5, 0.5, 0.5, 0.25, 0.75]), + np.array([0.2, -0.2, 0.3, -0.4, 0.1]), + np.array([0.1, 0.1, -0.2, 0.2, -0.3]), + ) + + wire_radius_sub = np.zeros_like(ifil_sub) + b_direct = cfsem.flux_density_linear_filament( + xyzp, xyzfil_sub, dlxyzfil_sub, ifil_sub, wire_radius_sub, par=False + ) + a_direct = cfsem.vector_potential_linear_filament( + xyzp, xyzfil_sub, dlxyzfil_sub, ifil_sub, wire_radius_sub, par=False + ) + + b_mlfmm, a_mlfmm = fields_linear_filament_mlfmm( + xyzp, + xyzfil, + dlxyzfil, + ifil, + eps, + use_linear_filament=True, + direct_threshold=1, + order=None, + ) + + for got, exp in zip(b_mlfmm, b_direct, strict=True): + assert np.allclose(got, exp, rtol=1e-6, atol=1e-10) + for got, exp in zip(a_mlfmm, a_direct, strict=True): + assert np.allclose(got, exp, rtol=1e-6, atol=1e-10) diff --git a/uv.lock b/uv.lock index abaac76f..78d5d7f4 100644 --- a/uv.lock +++ b/uv.lock @@ -1,5 +1,5 @@ version = 1 -revision = 2 +revision = 3 requires-python = ">=3.10, <3.14" resolution-markers = [ "python_full_version >= '3.11'", @@ -48,7 +48,6 @@ wheels = [ [[package]] name = "cfsem" -version = "3.1.0" source = { editable = "." } dependencies = [ { name = "findiff" }, @@ -75,7 +74,7 @@ dev = [ [package.metadata] requires-dist = [ { name = "findiff", specifier = ">=0.12.1" }, - { name = "interpn", extras = ["pydantic"], specifier = ">=0.8.2,<0.9" }, + { name = "interpn", extras = ["pydantic"], specifier = ">=0.8.2,<0.12" }, { name = "numpy", specifier = ">=2" }, { name = "pydantic", specifier = ">=2" }, { name = "pydantic-numpy", specifier = ">=6" }, diff --git a/vendor/armadillo-15.2.3.zip b/vendor/armadillo-15.2.3.zip new file mode 100644 index 00000000..99205f8d Binary files /dev/null and b/vendor/armadillo-15.2.3.zip differ diff --git a/vendor/boost-boost-1.90.0 b/vendor/boost-boost-1.90.0 new file mode 160000 index 00000000..94c16ace --- /dev/null +++ b/vendor/boost-boost-1.90.0 @@ -0,0 +1 @@ +Subproject commit 94c16ace0d426ffe169364197856be6bbc0af33b diff --git a/vendor/jsoncpp b/vendor/jsoncpp new file mode 160000 index 00000000..b511d9e6 --- /dev/null +++ b/vendor/jsoncpp @@ -0,0 +1 @@ +Subproject commit b511d9e64956db998b74909df112ac8c8f41d6ff diff --git a/vendor/rat-common b/vendor/rat-common new file mode 160000 index 00000000..a3211219 --- /dev/null +++ b/vendor/rat-common @@ -0,0 +1 @@ +Subproject commit a3211219951bffbcf0607c71055facb62855a4c7 diff --git a/vendor/rat-mlfmm b/vendor/rat-mlfmm new file mode 160000 index 00000000..d85cbee1 --- /dev/null +++ b/vendor/rat-mlfmm @@ -0,0 +1 @@ +Subproject commit d85cbee14dc090cd69b0859c7daeb8acd50a263d diff --git a/vendor/tclap b/vendor/tclap new file mode 160000 index 00000000..8b35dd1c --- /dev/null +++ b/vendor/tclap @@ -0,0 +1 @@ +Subproject commit 8b35dd1c23922231ec2dcb95edea0285189de039