Biot-Savart volume integral #996

ddudt · 2024-04-16T14:22:35Z

Computes B_plasma in a volume as:

𝐁ᵥ(𝐫) = μ₀/4π ∫ 𝐉(𝐫') × (𝐫 − 𝐫')/|𝐫 − 𝐫'|³ d³𝐫'

This is useful for when you need B_total = B_plasma + B_external in the full computational domain, not just on the boundary surface.

codecov · 2024-04-16T15:39:04Z

Codecov Report

❌ Patch coverage is 96.40411% with 21 lines in your changes missing coverage. Please review.
✅ Project coverage is 94.56%. Comparing base (7aa703f) to head (8f18cf4).
⚠️ Report is 1 commits behind head on master.

Files with missing lines	Patch %	Lines
desc/basis.py	95.23%	7 Missing ⚠️
desc/grid.py	95.31%	6 Missing ⚠️
desc/magnetic_fields/_core.py	75.00%	5 Missing ⚠️
desc/magnetic_fields/_plasma_field.py	97.33%	2 Missing ⚠️
desc/equilibrium/equilibrium.py	98.24%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #996      +/-   ##
==========================================
+ Coverage   94.52%   94.56%   +0.03%     
==========================================
  Files         102      105       +3     
  Lines       28712    29272     +560     
==========================================
+ Hits        27141    27680     +539     
- Misses       1571     1592      +21

Files with missing lines	Coverage Δ
desc/coils.py	`97.96% <100.00%> (-0.01%)`	⬇️
desc/compute/_fast_ion.py	`100.00% <100.00%> (ø)`
desc/compute/nabla.py	`100.00% <100.00%> (ø)`
desc/integrals/virtual_casing.py	`100.00% <100.00%> (ø)`
desc/magnetic_fields/__init__.py	`100.00% <100.00%> (ø)`
desc/transform.py	`94.21% <100.00%> (+1.56%)`	⬆️
desc/equilibrium/equilibrium.py	`96.21% <98.24%> (+0.13%)`	⬆️
desc/magnetic_fields/_plasma_field.py	`97.33% <97.33%> (ø)`
desc/magnetic_fields/_core.py	`96.10% <75.00%> (-0.52%)`	⬇️
desc/grid.py	`94.76% <95.31%> (+0.03%)`	⬆️
... and 1 more

... and 3 files with indirect coverage changes

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

github-actions · 2024-04-16T18:39:34Z

|             benchmark_name             |         dt(%)          |         dt(s)          |        t_new(s)        |        t_old(s)        | 
| -------------------------------------- | ---------------------- | ---------------------- | ---------------------- | ---------------------- |
 test_build_transform_fft_lowres         |     +0.65 +/- 3.17     | +3.28e-03 +/- 1.60e-02 |  5.08e-01 +/- 1.0e-02  |  5.05e-01 +/- 1.3e-02  |
 test_build_transform_fft_midres         |     +0.48 +/- 1.36     | +2.80e-03 +/- 8.02e-03 |  5.91e-01 +/- 4.7e-03  |  5.89e-01 +/- 6.5e-03  |
 test_build_transform_fft_highres        |     +0.19 +/- 2.26     | +1.84e-03 +/- 2.21e-02 |  9.79e-01 +/- 1.0e-02  |  9.77e-01 +/- 2.0e-02  |
 test_equilibrium_init_lowres            |     +0.51 +/- 0.61     | +1.86e-02 +/- 2.21e-02 |  3.63e+00 +/- 1.5e-02  |  3.61e+00 +/- 1.6e-02  |
 test_equilibrium_init_medres            |     +0.10 +/- 1.63     | +4.19e-03 +/- 6.64e-02 |  4.09e+00 +/- 5.4e-02  |  4.09e+00 +/- 3.9e-02  |
 test_equilibrium_init_highres           |     +0.32 +/- 0.52     | +1.76e-02 +/- 2.87e-02 |  5.50e+00 +/- 2.1e-02  |  5.48e+00 +/- 1.9e-02  |
 test_objective_compile_dshape_current   |     +0.18 +/- 0.92     | +6.74e-03 +/- 3.45e-02 |  3.76e+00 +/- 2.5e-02  |  3.75e+00 +/- 2.4e-02  |
 test_objective_compile_atf              |     +0.08 +/- 1.25     | +6.46e-03 +/- 1.02e-01 |  8.12e+00 +/- 7.9e-02  |  8.11e+00 +/- 6.4e-02  |
 test_objective_compute_dshape_current   |     +1.45 +/- 5.20     | +1.82e-05 +/- 6.52e-05 |  1.27e-03 +/- 5.9e-05  |  1.26e-03 +/- 2.7e-05  |
 test_objective_compute_atf              |     -0.32 +/- 4.53     | -1.38e-05 +/- 1.92e-04 |  4.23e-03 +/- 1.5e-04  |  4.25e-03 +/- 1.2e-04  |
 test_objective_jac_dshape_current       |     -5.67 +/- 11.28    | -2.18e-03 +/- 4.33e-03 |  3.62e-02 +/- 3.1e-03  |  3.84e-02 +/- 3.0e-03  |
 test_objective_jac_atf                  |     +0.14 +/- 2.41     | +2.62e-03 +/- 4.53e-02 |  1.89e+00 +/- 3.8e-02  |  1.88e+00 +/- 2.4e-02  |
 test_perturb_1                          |     +0.50 +/- 0.52     | +6.47e-02 +/- 6.73e-02 |  1.30e+01 +/- 4.8e-02  |  1.30e+01 +/- 4.7e-02  |
 test_perturb_2                          |     +0.51 +/- 0.55     | +9.10e-02 +/- 9.81e-02 |  1.79e+01 +/- 8.3e-02  |  1.78e+01 +/- 5.3e-02  |
 test_proximal_jac_atf                   |     +0.16 +/- 0.96     | +1.18e-02 +/- 7.01e-02 |  7.31e+00 +/- 5.7e-02  |  7.30e+00 +/- 4.1e-02  |
 test_proximal_freeb_compute             |     +1.08 +/- 0.70     | +1.90e-03 +/- 1.23e-03 |  1.78e-01 +/- 7.3e-04  |  1.76e-01 +/- 9.8e-04  |
 test_proximal_freeb_jac                 |     +0.39 +/- 0.99     | +2.83e-02 +/- 7.19e-02 |  7.31e+00 +/- 5.2e-02  |  7.28e+00 +/- 4.9e-02  |
 test_solve_fixed_iter                   |     +0.80 +/- 8.77     | +1.17e-01 +/- 1.29e+00 |  1.48e+01 +/- 9.3e-01  |  1.47e+01 +/- 8.9e-01  |

dpanici · 2024-04-24T19:48:40Z

use this fxn

DESC/desc/magnetic_fields/_core.py

Line 31 in b10f8c8

def biot_savart_general(re, rs, J, dV):

dpanici · 2024-07-03T20:01:15Z

scipy.integrate.tpl_quad to check this

dpanici · 2024-07-03T20:01:40Z

IAEA2020_JGeiger_Manuscript_8p_finalversion (1).pdf
Contributions to Plasma Physics - 2010 - Geiger - Effects of Net Currents on the Magnetic Configuration of W7‐X (2).pdf

f0uriest · 2024-07-10T20:45:06Z

Several people at the scidac meeting today are using the BMW code: https://github.com/ORNL-Fusion/BMW/tree/master

refs: https://scipub.euro-fusion.org/wp-content/uploads/eurofusion/WPS1CPR17_17542_submitted-4.pdf
https://bpb-us-e2.wpmucdn.com/wordpress.auburn.edu/dist/5/118/files/2022/05/ISHW2017.pdf

basically it does a biot-savart volume integral to compute A (not B), then takes the finite difference curl of A to get B. As far as I can tell from the code the source grid is hardcoded to ntheta=101, ns, nzeta taken from the vmec/mgrid file. Allegedly it works well and gives accurate results?

Could try something similar (and maybe use spectral derivatives for going A->B)

ddudt · 2024-07-18T04:00:08Z

Try implementing Eq. 17-18 of Peraza-Rodriguez et al. Phys. Plasmas 2017
It should be a simple linear system solve to get the Fourier-Zernike coefficients of the vector potential as:
$\nabla \times (\nabla \times \mathbf{A}) = \mu_0 \mathbf{J}$
Eq. 16 might be needed for boundary conditions on $\mathbf{A}$ at $\rho=0$ and $\rho=1$. Then the magnetic field is:
$\mathbf{B} = \nabla \times \mathbf{A}$
And all of this can be done in the native $(\rho,\theta,\zeta)$ coordinate system.

github-actions · 2025-05-21T21:05:28Z

Memory benchmark result

|               Test Name                |      %Δ      |    Master (MB)     |      PR (MB)       |    Δ (MB)    |    Time PR (s)     |  Time Master (s)   |
| -------------------------------------- | ------------ | ------------------ | ------------------ | ------------ | ------------------ | ------------------ |
  test_objective_jac_w7x                 |    4.53 %    |     3.796e+03      |     3.968e+03      |    171.93    |       39.19        |       36.44        |
  test_proximal_jac_w7x_with_eq_update   |    0.82 %    |     6.440e+03      |     6.492e+03      |    52.57     |       161.59       |       161.18       |
  test_proximal_freeb_jac                |   -0.29 %    |     1.324e+04      |     1.320e+04      |    -38.53    |       82.12        |       83.13        |
  test_proximal_freeb_jac_blocked        |   -0.81 %    |     7.549e+03      |     7.487e+03      |    -61.39    |       72.69        |       72.61        |
  test_proximal_freeb_jac_batched        |    0.37 %    |     7.465e+03      |     7.493e+03      |    27.90     |       72.59        |       71.79        |
  test_proximal_jac_ripple               |    1.04 %    |     3.458e+03      |     3.494e+03      |    35.98     |       64.79        |       65.09        |
  test_proximal_jac_ripple_bounce1d      |   -0.96 %    |     3.517e+03      |     3.483e+03      |    -33.85    |       76.15        |       76.24        |
  test_eq_solve                          |    5.97 %    |     2.016e+03      |     2.136e+03      |    120.25    |       93.58        |       93.37        |

For the memory plots, go to the summary of Memory Benchmarks workflow and download the artifact.

…nto dd/B_plasma

ddudt · 2025-06-03T00:16:09Z

Update: I fixed a silly bug and now this is working. Direct Biot-Savart volume integral to compute B and the similar volume integral to compute A both agree well with BMW results at moderate resolution. But more benchmarking is needed.

rahulgaur104 · 2025-06-03T00:31:05Z

Thank you! Could be really useful for my umbilic stuff!

dpanici · 2025-06-05T19:25:40Z

desc/integrals/singularities.py

+
+    # shape = (source_grid.num_nodes, eval_grid.num_nodes, 3)  # noqa: E800
+    dx = rpz2xyz(eval_grid)[None, :, :] - rpz2xyz(data["x"])[:, None, :]
+    B = np.sum(  # sum over source_grid nodes


use biot_savart_general instead

dpanici · 2025-06-05T19:26:12Z

check error for numpy vs jax.numpy
show plots of error
errors may be due to large matrix dot product errors?

ddudt · 2025-06-05T19:52:28Z

Possibly redundant with #1436, should merge them once we decide on a public API.

dpanici · 2025-09-08T18:47:11Z

desc/basis.py

            self._set_up()


+class DoubleChebyshevFourierBasis(_Basis):


Do we want to make things clearer that this is a real space basis (as opposed to flux coordinates) in the name/docstring?

dpanici · 2025-09-08T18:51:00Z

desc/basis.py

+        Maximum vertical resolution.
+    NFP : int
+        Number of field periods.
+    sym : {``'cos'``, ``'sin'``, ``False``}


there is not 't' in real space, what do these symmetries refer to now?

dpanici · 2025-09-08T18:53:50Z

desc/basis.py

+    ----------
+    L : int
+        Maximum radial resolution.
+    M : int


I don't love this, almost everywhere in the code N is toroidal, yet here we change it so M is toroidal. I would prefer if here we don't use M, and keep N as toroidal and use something else for the vertical resolution (Or be explicit and write N_R N_phi N_Z as the resolutions

dpanici · 2025-09-08T18:54:19Z

desc/basis.py

+        Parameters
+        ----------
+        grid : Grid or ndarray of float, size(num_nodes,3)
+            Node coordinates, in (rho,phi,Z).


dpanici · 2025-09-08T18:54:25Z

desc/basis.py

+        grid : Grid or ndarray of float, size(num_nodes,3)
+            Node coordinates, in (rho,phi,Z).
+        derivatives : ndarray of int, shape(num_derivatives,3)
+            Order of derivatives to compute in (rho,phi,Z).


dpanici · 2025-09-08T18:56:19Z

desc/basis.py

+            return np.array([]).reshape((grid.num_nodes, 0))
+
+        # Get the nodes of the grid
+        r, p, z = grid.nodes.T


nitpick but if r is really cylindrical R, use R, same for Z

dpanici · 2025-09-08T18:56:46Z

desc/basis.py

+        l, m, n = modes.T
+
+        # The coordinates in the grid are hardcoded to be called rho, theta, zeta,
+        # but these vars are just unique values for 1st, 2nd, and 3rd coordinates


This will change, presumedly, when we figure out our new grid stuff @ddudt #1840

dpanici · 2025-09-08T19:18:24Z

desc/basis.py

        )


+def chebfit(y, axis):


we may aready have some utility functions that do this, see

DESC/desc/integrals/basis.py

Line 260 in 702ae49

def compute_cheb(self, x):

dpanici · 2025-09-08T19:20:21Z

desc/basis.py

+    return idct(y_c / f.reshape(f_shape), axis=axis, type=1, norm=None)
+
+
+def fftfit(y, axis, n=None):


if these are to be public, some more safeguards should be added. If these are not meant to be public, add _ before the name

dpanici · 2025-09-08T19:21:21Z

desc/basis.py

+    """
+    Real fast fourier transform along axis.
+
+    Designed to convert coefficients into the form expected by the


what is meant here? takes in a function evaluated at linearly spaced nodes, then fourier transforms it, then adjusts the coefficients to fit what DESC's fourier function requires?

dpanici · 2025-09-08T19:21:49Z

desc/basis.py

+    Parameters
+    ----------
+    y : ndarray, shape(...,2N+1,...)
+        Function to decompose into basis of Chebyshev functions.


is this docstring accurate?

dpanici · 2025-09-08T19:24:52Z

desc/compute/nabla.py

+        transform from the spectral basis on which A is calculated to the
+        real grid on which the curl is to be evaluated.
+        the transform must have been built, with derivs>=1.
+    scales: jnp.ndarray, shape (3,)


shouldn't this be done in the transform itself, not here?

dpanici · 2025-09-08T19:25:19Z

desc/compute/nabla.py

+    Returns
+    -------
+    curl_A : ndarray, shape(n,3)
+        The curl of the vector field, in cylindrical coordinates.


Suggested change

The curl of the vector field, in cylindrical coordinates.

The curl of the vector field, in cylindrical coordinates, at the points specified by `out_transform`.

dpanici · 2025-09-08T19:25:38Z

desc/compute/nabla.py

+    out_R : ndarray, shape(n,)
+        radial location for each point at which to calculate
+        the curl. out_R should never be 0.
+    in_transform: Transform


specify these should be bases and grids of R, phi, Z

dpanici · 2025-09-08T19:36:00Z

desc/compute/nabla.py

+        for d in range(3):
+            if c == 1 and d == 0:
+                # partial (R*A_phi)/partial R instead of partial A_phi/partial R
+                term_cd = out_transform.transform(RA_phi_coeff, dr=1)


cant you just use chain rule and do d/dR (R Aphi) = Aphi + R d(Aphi)/dR? then you only need to do one fit above, and all of the partial derivatives are "consistent" in that they are all derivaitves of A, instead of derivatives of A and some derivatves of R*A (which may have a different spectral width than A i.e. same resolution basis may be a worse for for R*A then for A)

dpanici · 2025-09-08T19:40:58Z

desc/compute/nabla.py

+    A_coeff = in_transform.fit(
+        A * jnp.stack([in_R, jnp.ones_like(in_R), jnp.ones_like(in_R)], axis=-1)
+    )
+    # Calculate matrix of terms for the curl


Suggested change

# Calculate matrix of terms for the curl

# Calculate matrix of terms for the div

dpanici · 2025-09-08T19:42:45Z

desc/compute/nabla.py

@@ -0,0 +1,163 @@
+"""
+Functions for taking the curl and divergence in cylindrical coordinates.


usually desc.compute has things for computing stuff in the compute dicts. This I could see as useful in the future when we have things in the compute dict using RpZ grids (#1805 ), for now let's call this like _nabla_utils.py maybe or _real_space_deriv_utils.py

dpanici · 2025-09-08T19:45:26Z

desc/equilibrium/equilibrium.py

+        coords = jnp.atleast_2d(coords)
+
+        # Biot-Savart method
+        if method == methods[0]:


Suggested change

if method == methods[0]:

if method == "biot-savart":

I prefer to be verbose in the code here

dpanici · 2025-09-08T19:45:38Z

desc/equilibrium/equilibrium.py

+                B = xyz2rpz_vec(B, phi=coords[:, 1])
+
+        # Virtual casing method
+        elif method == methods[1]:


Suggested change

elif method == methods[1]:

elif method == "virtual casing":

dpanici · 2025-09-08T19:47:26Z

desc/integrals/virtual_casing.py

+from desc.utils import errorif, xyz2rpz_vec
+
+
+def integrate_surface(coords, source_data, source_grid, kernel, chunk_size=None):


maybe different name like integrate_surface_at_external_point

it does not need to be strictly outside per-say right? just needs to be off-surface, the kernels don't really care about outside or inside as far as the singularity goes

dpanici

partially reviewed but already some changes noted

github-actions · 2026-01-27T23:39:38Z

desc/magnetic_fields/_core.py

+        field = args[0]
+        r = rpz[:, 0]
+        br, bp, bz = (
+            scale


[flake8] <821> _{reported by reviewdog 🐶}
undefined name 'scale'

daniel-dudt and others added 2 commits April 16, 2024 09:19

initial commit

9ded7bf

Merge branch 'master' into dd/B_plasma

d41b81d

do calculations in Cartesian coordinates

6ef3279

dpanici and others added 2 commits June 20, 2024 10:49

Merge branch 'master' into dd/B_plasma

6afebef

Merge branch 'master' into dd/B_plasma

f6dab47

Merge branch 'master' into dd/B_plasma

814942d

ddudt self-assigned this May 21, 2025

daniel-dudt added 2 commits May 21, 2025 14:52

merge with master

c91b4d7

fix x_Equilibrium computation

0e911b0

daniel-dudt and others added 7 commits May 22, 2025 10:33

remove duplicate compute 'x' fun

d9d9926

Merge branch 'master' into dd/B_plasma

2e270e7

Merge branch 'master' into dd/B_plasma

9b1ac99

Merge branch 'master' into dd/B_plasma

c00a625

fix norm axis bug

5723f13

Merge branch 'dd/B_plasma' of https://github.com/PlasmaControl/DESC i…

6b5b8f3

…nto dd/B_plasma

replace for loop with broadcasting

facedf6

Merge branch 'master' into dd/B_plasma

36b2ea9

dpanici reviewed Jun 5, 2025

View reviewed changes

added DoubleChebyshevFourierBasis to all

d8214b6

maya-avida marked this pull request as ready for review August 29, 2025 15:24

Merge branch 'master' into dd/B_plasma

1c82b70

maya-avida requested a review from dpanici September 2, 2025 17:45

dpanici reviewed Sep 8, 2025

View reviewed changes

dpanici requested changes Sep 8, 2025

View reviewed changes

dpanici added the waiting for other PRs label Sep 24, 2025

dpanici marked this pull request as draft September 24, 2025 19:53

f0uriest mentioned this pull request Dec 3, 2025

Guiding Center Particle Tracing #1820

Merged

Merge branch 'master' into dd/B_plasma

8f18cf4

github-actions bot reviewed Jan 27, 2026

View reviewed changes

		return idct(y_c / f.reshape(f_shape), axis=axis, type=1, norm=None)


		def fftfit(y, axis, n=None):

	The curl of the vector field, in cylindrical coordinates.
	The curl of the vector field, in cylindrical coordinates, at the points specified by `out_transform`.

	# Calculate matrix of terms for the curl
	# Calculate matrix of terms for the div

		@@ -0,0 +1,163 @@
		"""
		Functions for taking the curl and divergence in cylindrical coordinates.

		from desc.utils import errorif, xyz2rpz_vec


		def integrate_surface(coords, source_data, source_grid, kernel, chunk_size=None):

Biot-Savart volume integral #996

Are you sure you want to change the base?

Biot-Savart volume integral #996

Conversation

ddudt commented Apr 16, 2024

Uh oh!

codecov bot commented Apr 16, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

github-actions bot commented Apr 16, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

dpanici commented Apr 24, 2024

Uh oh!

dpanici commented Jul 3, 2024

Uh oh!

dpanici commented Jul 3, 2024

Uh oh!

f0uriest commented Jul 10, 2024

Uh oh!

ddudt commented Jul 18, 2024 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

github-actions bot commented May 21, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Memory benchmark result

Uh oh!

ddudt commented Jun 3, 2025

Uh oh!

rahulgaur104 commented Jun 3, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dpanici commented Jun 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

ddudt commented Jun 5, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dpanici Sep 8, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

dpanici Sep 8, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

codecov bot commented Apr 16, 2024 •

edited

Loading

github-actions bot commented Apr 16, 2024 •

edited

Loading

ddudt commented Jul 18, 2024 •

edited

Loading

github-actions bot commented May 21, 2025 •

edited

Loading

dpanici commented Jun 5, 2025 •

edited

Loading

dpanici Sep 8, 2025 •

edited

Loading

dpanici Sep 8, 2025 •

edited

Loading