Cbs/fieldine integrator #555

smiet · 2025-09-19T15:01:29Z

Revamp of fieldline integration. Some tests still need to be written, but this should simplify the user experience in generating Poincare plots.

What’s new

Integrator base class
with subclasses:

SimsoptFieldlineIntegrator (interface to the sopp integration routines
ScipyFieldlineIntegrator (integration of the magnetic field using scipy.solve_ivp methods).

These subclasses provide methods to return Poincare sections, three dimensional field line trajectories, as well as the ability to integrate from any starting point over a given toroidal distance.

The ScipyFieldlineIntegrator provides methods to integrate $R$ and $Z$ using $\phi$ as the 'time' variable in the ODE, providing crisper Poincare sections and integrates about as fast as the sopp methods (which evolves the three coordinates independently).

The Integrator subclasses also provide an interface to the PoincarePlotter class, which provides an easy interface to create Poincare plots of the field.

The PoincarePlotter uses caching and the simsopt dependency graph (parent Integrator with parent MagneticField) to trigger recomputes when necessary.
Manipulation of the res_phi_hits and res_tys arrays is a thing of the past!

give it a spin!:

from simsopt.configs import get_data
from simsopt.field import SimsoptFieldlineIntegrator, PoincarePlotter
from simsopt.geo import plot
import numpy as np

base_coils, base_currents, ma, nfp, bs = get_data('ncsx')

axis_rz = ma.gamma()[0][::2]
# startpoints linspace from axis out: 
start_points_RZ = np.linspace(axis_rz, axis_rz+np.array([0.1, 0]), 10)


integrator = SimsoptFieldlineIntegrator(bs, stellsym=True, nfp=nfp, R0=axis_rz[0], tol=1e-9)
pplotter = PoincarePlotter(integrator, start_points_RZ, n_transits=200, phis=4)

fig, axs = pplotter.plot_poincare_all()


pplotter.plot_fieldline_trajectories_3d(engine='mayavi', show=False, opacity=0.2, tube_radius=0.001)
plot(bs.coils, engine='mayavi', show=False)
pplotter.plot_poincare_in_3d(engine='mayavi', show=True, scale_factor=0.01)

…ator

leading zero copied over from VMEC convention, causing all coefficients off-by-one. See issue #554

codecov · 2025-09-29T15:11:44Z

Codecov Report

❌ Patch coverage is 89.14186% with 62 lines in your changes missing coverage. Please review.
✅ Project coverage is 92.45%. Comparing base (0da1f07) to head (bdc0ab5).
⚠️ Report is 13 commits behind head on master.

Files with missing lines	Patch %	Lines
src/simsopt/field/tracing.py	89.14%	62 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #555      +/-   ##
==========================================
- Coverage   92.56%   92.45%   -0.12%     
==========================================
  Files          83       83              
  Lines       16229    16798     +569     
==========================================
+ Hits        15023    15530     +507     
- Misses       1206     1268      +62

Flag	Coverage Δ
unittests	`92.45% <89.14%> (-0.12%)`	⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

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

andrewgiuliani

Thanks Chris!

There are a few lines in code coverage that are missing, we should have unit tests for these.

By the way, we still cover plotting code in unit tests, but the canvas is just not plotted by skipping the mlab.show() command, like we do here

simsopt/src/simsopt/geo/surface.py

Line 266 in d053f90

mlab.show()

BoozerMagneticField is not in cartesian coordinates, and the Integrator subclasses are not set up to handle it.

smiet · 2025-10-22T13:01:54Z

@andrewgiuliani Thank you for your review, I have added tests to cover all the missed lines that I could. The ones that are not hit are:

the mayavi plotting routines (because we do not have a runner with mayavi installed, but they will run when that is the case)
The lines that gather data from several processes (though these should be hit in the examples, which are run with multiprocessing, but maybe coverage is unaware)
a few others for which I could not find an elegant way to hit, like the exception to not make proc0print fail on non-mpi but mpi-capable systems.

I hope that we can accept the few misses, the diff is only a few percent below the current coverage.

There are some other failing tests because SPEC install on pytnon 3.9 is currently broken.

@landreman @mishapadidar, could I ask you for a review if you have a bit of time to spare? A few people already using it have told me this is very pleasant to work with, and I would like to make these features available in the not-too-distant-future.

I know that there are many more features to add, but I will deal with these in later PRs, otherwise this one will bloat even larger than it currently is.

Future plans include:

make InterpolatedField aware of its underlying MagneticField DoF's and trigger re-computation on first access after DoF change.
include a SimsoptParticleIntegrator (requires me to first familiarize myself with how particle tracing works... or I could delegate...)
deal with BoozerMagneticField in integration and in particle integration.
methods for Integrator to generate Surface and Curve objects from the integration results.

src/simsopt/util/mpi.py

mishapadidar · 2025-10-22T13:28:01Z

src/simsopt/field/tracing.py

+    and being an Optimizable, the results of integrations can be
+    optimized. 
+    """
+    def __init__(self, field: MagneticField, comm=None, nfp=None, stellsym=False, R0=None, test_symmetries=True):


can we remove test_symmetries?

As in completely remove the functionality?
We can, but I thought it prudent to have a few basic checks before we create an object that a user will blindly trust. It is too easy to pass the wrong NFP or assume something is stellarator-symmetric.
These properties are relied upon by PoincarePlotter, and the test is also super fast compared to the time it takes to integrate.

I would like to keep this functionality until in a next PR we can move it to the MagneticField class (which should carry information on NFP, symmetries, and a method to test them).

src/simsopt/field/tracing.py

mishapadidar

The Integrator class is proving two general benefits: it makes integrators Optimizable and it provides convenience methods for unpacking the res_phi_hits data. Not all of the instances of the class adhere to this framework, so I suggest we make an attempt to align the usage. In particular, currently Simsopt traces fieldlines through the compute_fieldlines function which calls the sopp integrators. The development of the ScipyFieldlineIntegrator class suggests that there may be desire for using scipy integrators in place of the sopp integrators. To align with the current simsopt style, we should implement a function similar to sopp.fieldline_tracing but for the scipy integrators. Then we modify the compute_fieldlines function to accept a method=... argument. Finally, the Integrator class can take a method=... argument as well. This would allow us to eliminate the instances of the FieldlineIntegrator that are specific to integration methods, i.e. the ScipyFieldLineInegrator, sliming the code down to have a single FieldLineIntegrator class. In this framework, other instances of the Integrator class would solve other ODES. For example, a GuidingCenterIntegrator class.

A second point is that, since we are setting up field line integration as an Optimizable, we really should be implementing derivatives. Differentiating the field line integration is tractable and greatly enhances the usability and applicability of this class.

I have added a lot of comments, and I have not finished reviewing the PoincarePlotter.
Here are some overarching themes to the comments.

Add clean doc strings to all functions with a succinct description of the method, arguments and returns.
The integrator classes have some convenience functions which lead to duplication. Consider removing the duplication.

src/simsopt/field/tracing.py

mishapadidar · 2025-10-22T13:42:32Z

src/simsopt/field/tracing.py

+        if not periodicity_test:
+            raise ValueError(f"The field does not adhere to the {self.nfp} field periods it is set to be.")
+
+    def Interpolate_field(self, rrange, phirange, zrange, degree=2, skip=None, extrapolate=True): 


The main purpose of this function is to build an interpolatedfield from a field. Since the method only depends on the field and not the inegrator, this convenience function would be most suited to the MagneticField class rather than the Integrator class.
Removing this method simplifies the Integrator class and should not effect user experience much (they could initialize the Integrator from the interpolatedField easily).
Consider removing this method.

I remove this classmethod and will in the future make InterpolatedField self-updating.

src/simsopt/field/tracing.py

mishapadidar · 2025-10-22T14:24:12Z

src/simsopt/field/tracing.py

+    def generate_symmetry_planes(phis, nfp=1):
+        """
+        Given a list of phis in [0, 2pi/nfp], generate the full list of phis
+        in [0, 2pi] by adding the symmetry planes. 


mishapadidar · 2025-10-22T14:24:20Z

src/simsopt/field/tracing.py

+            return self._randomcolors
+
+    @staticmethod
+    def generate_phis(nplanes, nfp=1):


mishapadidar · 2025-10-22T14:25:30Z

src/simsopt/field/tracing.py

+        return self._res_tys
+
+    @property
+    def res_phi_hits(self):


mishapadidar · 2025-10-22T14:25:56Z

src/simsopt/field/tracing.py

+        """
+        Generate a hash from dofs, self.phis, self.start_points_RZ, and self.n_transits
+        """
+        hash_list = self.integrator.field.full_x.tolist() + self.phis.tolist() + self.start_points_RZ.flatten().tolist() + [self.n_transits]


use the uuid class or similar

That would cause recompute if a second plotter with an identical field and startpoints, wanted to plot the same. Or in the next session or different kernel. This way, if a file has been generated, and has the same field with dofs, it will not recompute. Except if pythons hash function changes (it had recently), but that is too much of an edge case.

mishapadidar · 2025-10-22T14:27:49Z

src/simsopt/field/tracing.py

+        poincare_hash = hash(tuple(hash_list))
+        return poincare_hash
+
+    def save_to_disk(self, filename=None, name=None):


It seems that the functionality here is to initialize a poincare plotter from data. Rather than save the class, which potentially conflicts with Simsopt save/load structure, we could write a from_phi_hits method than initializes the class from the necessary data.

I can call it differently, this is indeed confusing. It is only to store/retrieve poincare data to/from a file, not the entire object.

smiet · 2025-10-23T11:49:04Z

Thanks @mishapadidar for your points so far. You raise a few points for discussion:

Implement `compute_fieldlines` with the scipy integration?

I would like to move away completely from the functional paradigm for Poincare plotting, and suggest users to only use class-based Integrators. I left the compute_fieldlines method for reverse compatibility, but want to deprecate that method of performing integration. @landreman @andrewgiuliani, what do you think, update the compute_fieldlines function or move towards deprecating it (and leave it be for now)?

Limiting the scope of each integrator

Scipy and Simsopt fieldline integrators perform different computations, and use different methods, but should we limit each integrator to only one way of solving the equations? Should we merge the two integrators and let it switch between whether sopp or scipy is used based on the requested data from the PoincarePlotter?

I would hope to leave this work to a later PR, and I am not equipped to modify the sopp integrator to implement the features I want

Derivatives of integration

This is in the pipeline. I am trying to get this PR out to cement the class interfaces, and am developing a topological optimisation toolset in parallel that will do this. This will involve giving the ScipyFieldlineIntegrator methods to differentiate integrations, evaluate action integrals, etc. and defining new Optimizables with proper J and dJ functions.

smiet · 2026-01-07T14:59:00Z

@andrewgiuliani @mishapadidar @landreman Think we're good to go! lmk if more changes are needed.

smiet added 17 commits May 20, 2025 17:41

start on integrator and poincareplotting classes

c57c28e

further work on SimsoptIntegrator

2b30d48

Merge remote-tracking branch 'origin/master' into cbs/fieldine_integr…

12c3564

…ator

Merge remote-tracking branch 'origin/master' into cbs/fieldine_integr…

6cbd7c9

…ator

tracer updates

635af31

Merge remote-tracking branch 'origin/master' into cbs/fieldine_integr…

55d2d3a

…ator

update tracing routines

5df54fd

update tracing fieldlines example to use integrator

d958c68

update QA tracing example

9363e99

few more fixes to tracing

69c68ee

add docstring to proc0_print

9e2236d

updates to tracing

58235b7

test the integrators

34becef

use the scipy integrator for the QA example

c89a188

harmonize angle integration call signature

0415134

improve test clarity

ae06e3e

streamline call signature for integrate_in_phi

6847fe0

smiet mentioned this pull request Sep 29, 2025

w7x-axis does not conform to field line tracing #554

Closed

smiet added 5 commits September 29, 2025 15:47

Correct w7x axis coefficients

bc6b75c

leading zero copied over from VMEC convention, causing all coefficients off-by-one. See issue #554

fix color kwarg in plotting and from_field classmethod signature

dd59fc7

fix tests, do not skip w7x

35b99ef

tests for PoincarePlotter class

f150e21

fix lint

c297ba7

smiet requested a review from mishapadidar September 29, 2025 14:00

fix missed lint complaint

bca7ab1

smiet marked this pull request as ready for review September 29, 2025 14:02

smiet added 3 commits September 29, 2025 17:30

fix example fail when comm_world is None

fd58a46

leave recomputing interpolated field to future me

d57f950

add pplotter saving to and loading from file with store_results flag

d170c90

andrewgiuliani requested changes Oct 6, 2025

View reviewed changes

smiet added 9 commits October 9, 2025 16:42

use mpi in integrator examples

e718bac

skip input validation for integrator interpolation fn

5abad1d

fix [0,1] simsopt angle for surface in pplotter

65fbcdc

remove BoozerMagneticField integration and add tests

c30861e

BoozerMagneticField is not in cartesian coordinates, and the Integrator subclasses are not set up to handle it.

moar coverage

b10f954

fix examples to run on systems without MPIPartition

ac44a18

test most all untested lines

e584014

fix lints

9550cea

test mark lost plotting

521ad3d

smiet requested a review from andrewgiuliani October 22, 2025 13:07

mishapadidar reviewed Oct 22, 2025

View reviewed changes

src/simsopt/util/mpi.py Show resolved Hide resolved

mishapadidar reviewed Oct 22, 2025

View reviewed changes

src/simsopt/field/tracing.py Outdated Show resolved Hide resolved

mishapadidar reviewed Oct 22, 2025

View reviewed changes

smiet added 5 commits October 22, 2025 16:52

expand proc0_print docstring

bc5d0eb

clarify Integrator docstring

8464c98

mishas' suggestions implemented

bdc0ab5

better docstrings

bc0aac5

fix all docstrings

7af9d44

fix poincare load/save name usage

fadeed0

smiet requested a review from mishapadidar October 27, 2025 16:15

smiet added 5 commits January 7, 2026 11:51

comment and fix phi0

88551d3

Merge remote-tracking branch 'origin' into cbs/fieldine_integrator

fe770e5

reduce code duplication in scipy integrator pp gen

960765f

fix last tests

055a376

ruff ruff (lintfix)

e6493e2

Cbs/fieldine integrator #555

Are you sure you want to change the base?

Cbs/fieldine integrator #555

Uh oh!

Conversation

smiet commented Sep 19, 2025

Uh oh!

codecov bot commented Sep 29, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

andrewgiuliani left a comment • 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!

smiet commented Oct 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

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!

Uh oh!

mishapadidar left a comment • 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!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

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!

smiet Oct 22, 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!

smiet commented Oct 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Implement compute_fieldlines with the scipy integration?

Limiting the scope of each integrator

Derivatives of integration

Uh oh!

smiet commented Jan 7, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Sep 29, 2025 •

edited

Loading

andrewgiuliani left a comment •

edited

Loading

smiet commented Oct 22, 2025 •

edited

Loading

mishapadidar left a comment •

edited

Loading

smiet Oct 22, 2025 •

edited

Loading

smiet commented Oct 23, 2025 •

edited

Loading

Implement `compute_fieldlines` with the scipy integration?