Spec refactoring/improvement (PR 1 of 3 for single-stage) #418

smiet · 2024-05-30T13:53:13Z

Updates to Spec and NormalField and a few helper functions .

Spec class has been updated with the following features:

detection when free-boundary solution is not converged, re-set boundary field guess when detected.
default free-boundary configuration and classmethod to generate.
re-written loops pythonically where possible (for vol in range(self.inputlist.mvol+1) -> for surface in self.initial_guess)
implemented array_translator to reduce looping over index counters.
activate_profile helper function.

The NormalField class has been augmented to:

Always include the boundary on which it is defined (reference to associated SPEC boundary object).
Store the Fourier components that describe it in array-form.
include getters and setters for the FCs in array form.
set the dofs in their entirety even when setting single (fixes bug where changing a DOF did not trigger recompute bell).

The Surface base class has two new methods: fourier_transform_field and inverse_fourier_transform_field to efficiently calculate fourier transforms in the VMEC convention (-ntor to ntor toroidal components, 0:mpol poloidal). Allows for different Fourier normalization conventions through the 'normalization' parameter.

This is PR 1 of 3, the following two:

CoilSet and ReducedCoilSet: helper classes to enable interchanging of coils, and create new optimization variables through SVD on coil DOFS
CoilNormalField and misc changes: class inheriting form NormalField that interfaces to the CoilSets, all other misc changes needed to make the single-stage optimizations work.

codecov · 2024-05-30T15:04:58Z

Codecov Report

Attention: Patch coverage is 92.56018% with 34 lines in your changes missing coverage. Please review.

Project coverage is 92.52%. Comparing base (2c0f057) to head (a2a05ae).

Files	Patch %	Lines
src/simsopt/mhd/spec.py	90.29%	23 Missing ⚠️
src/simsopt/geo/surfacerzfourier.py	88.70%	7 Missing ⚠️
src/simsopt/field/normal_field.py	97.41%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #418      +/-   ##
==========================================
+ Coverage   92.00%   92.52%   +0.52%     
==========================================
  Files          75       75              
  Lines       13504    13833     +329     
==========================================
+ Hits        12424    12799     +375     
+ Misses       1080     1034      -46

Flag	Coverage Δ
unittests	`92.52% <92.56%> (+0.52%)`	⬆️

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.

landreman

Having the new Fourier transform and inverse transform functions go in the base Surface class feels awkward to me. These functions need mpol and ntor, which not all Surface classes have (e.g. SurfaceHenneberg, SurfaceGarabedian). Also the transform and inverse-transform functions don't use any of the surface geometry, just the grid points. So I would suggest moving those two functions to simsopt.util (probably not in a class, just standalone functions).

It also feels awkward for NormalField to carry around computational_boundary. See the comments attached to specific line numbers for some specific issues.

Codecov is flagging a lot of lines not covered by tests, 71% coverage for the changes compared to 91% for simsopt overall.

landreman · 2024-06-02T19:23:32Z

src/simsopt/field/normal_field.py

+
+    The Fourier harmonics are the degees of freedom, the computational boundary is kept fixed. 
+    The normal field should not be normalized to unit area, i.e. it is the 
+    fourier components of B.(\grad\theta \times \nabla\zeta) on the surface.


$(\partial \vec{r}/\partial\theta) \times (\partial \vec{r}/\partial\phi)$ rather than $\grad\theta \times \nabla\zeta$

But stepping back, why not choose the opposite convention, B dot n where n is the unit vector? Then computational_boundary.normal() could be removed on line 487, making it mostly unnecessary for NormalField to keep or refer to computational_boundary. The only other dependence is that get_real_space_field() uses quadrature points but those could be arguments.

I use this convention because this is the way that SPEC requires the Fourier components of the normal field. This is currently the only use-case for NormalField, and it is simplest if the degrees-of-freedom correspond directly with the variables passed to SPEC.
It is possible to remove the convention and add have SPEC do the translation. or to pass to the class a convention='SPEC', which would raise a NotImplementedError if another convention is requested. would that be a solution?

I do prefer to have the NormalField carry a reference to the surface on which it is calculated; without this the DoF's of NormalField are meaningless. If the user needs to keep track of which surface object the field is evaluated on, it makes sense to me to have NormalField to reference it. Modifying its DoFs I agree is not warranted.
The surface info is also used in helper functions such as get_real_space_field, and I would consider adding a .plot funcionality where it plots the field on the surface.

landreman · 2024-06-02T19:39:47Z

src/simsopt/geo/surface.py

@@ -839,6 +839,130 @@ def interpolate_on_arclength_grid(self, function, theta_evaluate):

        return function_interpolated

+    def fourier_transform_field(self, field, mpol=None, ntor=None, normalization=None, **kwargs):


Throughout this function and the inverse-transform function, the use of the word "field" is potentially confusing. It is natural to interpret it to mean the magnetic field, which is a vector. But I guess you mean any scalar field? It would be good to clarify this in the function names, docstrings, etc.

I have changed the name to scalar as it is a general method of transforming scalars, and adapted it in the docstrings. It is now also moved to surfaceRZFourier and not the Surface base class.

The calculation of such a transform always requires information from the surface on which the scalar quantity is evaluated. The function could be in util, but I find it conceptually cleaner to have it as a method of the surfaceRZFourier.

tests/geo/test_surface_rzfourier.py

landreman · 2024-06-02T20:10:37Z

tests/geo/test_surface_rzfourier.py

+
+        # Test that all other elements are close to zero
+        sines_mask = np.ones_like(ft_sines, dtype=bool)
+        cosines_mask = sines_mask


You want a copy() here I think. Presently cosines_mask and sines_mask are two names for the same single object, so setting an entry in one modifies both.

tests/mhd/test_spec.py

src/simsopt/mhd/spec.py

src/simsopt/field/normal_field.py

landreman · 2024-06-02T20:45:17Z

src/simsopt/geo/surface.py

+                if not stellsym:
+                    cosangle = np.cos(angle)
+                field = field + A_mns[m, n + ntor] * sinangle
+                if not stellsym:
+                    field = field + A_mnc[m, n + ntor] * cosangle


The two if statements can be combined into one.

src/simsopt/geo/surface.py

…rsists in tests)

smiet · 2024-06-17T14:01:10Z

I have adapted most changes and would like to move the review forward.

Almost all new untested lines are now covered.

Especially awaiting comments including a SurfaceRZFourier object with the NormalField. I think it is useful to have here, but especially in the next PR, where I create a CoilNormalField that inherits from NormalField it will involve many more functions that have to access the Surface (to calculate the coefficients etc).

Adapted NormalField.from_spec method to generate a SurfaceRZFourier class from the SPEC data.

P.S. I might have accidentally 'started a review' myself, my apologies.

…py 2.0

landreman

Hi @smiet, I see the testing was expanded but some of the issues above are still open. The main one is the location of the Fourier transform routines - I still think it doesn't make sense for them to be in the base Surface class, for the reasons explained above. There are also the issues about explaining that "field" is not a magnetic field vector, cosines_mask = sines_mask, and combining if statements

landreman · 2024-06-24T14:31:57Z

tests/geo/test_surface_rzfourier.py

+
+        # create a test field where only Fourier elements [m=2, n=3] 
+        # and [m=4,n=5] are nonzero:
+        field = field = 0.8 * np.sin(2*theta2d - 3*s.nfp*phi2d) + 0.2*np.sin(4*theta2d - 5*s.nfp*phi2d)+ 0.7*np.cos(3*theta2d - 3*s.nfp*phi2d)


Duplicate field =

landreman · 2024-06-24T14:44:20Z

src/simsopt/field/normal_field.py

            vnc = np.zeros((self.mpol + 1, 2 * self.ntor + 1))
+
+        if surface is None:
+            from simsopt.geo import SurfaceRZFourier


Why not put this import with the other imports at the top of the file?

landreman · 2024-06-24T14:44:45Z

src/simsopt/field/normal_field.py

@@ -91,6 +106,7 @@ def from_spec(cls, filename):
        if py_spec is None:
            raise RuntimeError(
                "Initialization from Spec requires py_spec to be installed.")
+        from simsopt.geo import SurfaceRZFourier


Why not put this import with the others at the top of the file?

landreman · 2024-06-24T14:45:36Z

src/simsopt/mhd/spec.py

+            copy_to_pwd: boolean, if True, the default input file will be copied to the current working directory. Has to be set True as free-boundary SPEC can only handle files in the current working directory.
+            verbose: boolean, if True, print statements will be printed
+        """
+        import shutil, os


Imports usually go at the top of the file unless necessary

landreman · 2024-06-24T14:46:11Z

tests/mhd/test_spec.py

+        test the default freeboundary file, and re-set if Picard
+        iteration fails
+        """
+        from simsopt._core.util import ObjectiveFailure


Imports to top of file

smiet

Hi Matt,

I am sorry, had entered many comments and discussion that are pending (aparently not for a response from you, but not submitted and hence not visible to you!

My apologies, submitting these now and I will put the imports at the top of the file where requested.

smiet · 2024-06-12T14:35:28Z

src/simsopt/field/normal_field.py

+
+    The Fourier harmonics are the degees of freedom, the computational boundary is kept fixed. 
+    The normal field should not be normalized to unit area, i.e. it is the 
+    fourier components of B.(\grad\theta \times \nabla\zeta) on the surface.


I use this convention because this is the way that SPEC requires the Fourier components of the normal field. This is currently the only use-case for NormalField, and it is simplest if the degrees-of-freedom correspond directly with the variables passed to SPEC.
It is possible to remove the convention and add have SPEC do the translation. or to pass to the class a convention='SPEC', which would raise a NotImplementedError if another convention is requested. would that be a solution?

I do prefer to have the NormalField carry a reference to the surface on which it is calculated; without this the DoF's of NormalField are meaningless. If the user needs to keep track of which surface object the field is evaluated on, it makes sense to me to have NormalField to reference it. Modifying its DoFs I agree is not warranted.
The surface info is also used in helper functions such as get_real_space_field, and I would consider adding a .plot funcionality where it plots the field on the surface.

src/simsopt/field/normal_field.py

smiet · 2024-06-14T08:15:55Z

src/simsopt/geo/surface.py

@@ -839,6 +839,130 @@ def interpolate_on_arclength_grid(self, function, theta_evaluate):

        return function_interpolated

+    def fourier_transform_field(self, field, mpol=None, ntor=None, normalization=None, **kwargs):


I have changed the name to scalar as it is a general method of transforming scalars, and adapted it in the docstrings. It is now also moved to surfaceRZFourier and not the Surface base class.

The calculation of such a transform always requires information from the surface on which the scalar quantity is evaluated. The function could be in util, but I find it conceptually cleaner to have it as a method of the surfaceRZFourier.

src/simsopt/geo/surface.py

smiet · 2024-06-14T08:22:08Z

src/simsopt/geo/surface.py

+                if not stellsym:
+                    cosangle = np.cos(angle)
+                field = field + A_mns[m, n + ntor] * sinangle
+                if not stellsym:
+                    field = field + A_mnc[m, n + ntor] * cosangle


tests/mhd/test_spec.py

smiet · 2024-06-14T09:55:38Z

tests/geo/test_surface_rzfourier.py

+
+        # Test that all other elements are close to zero
+        sines_mask = np.ones_like(ft_sines, dtype=bool)
+        cosines_mask = sines_mask


src/simsopt/mhd/spec.py

src/simsopt/field/normal_field.py

…for_numpy_2

… cbs/spec_update_rebase

…_2' into cbs/spec_update_rebase" This reverts commit a193ade, reversing changes made to 7c2a1f9.

mbkumar · 2024-06-27T18:57:31Z

@smiet
Can you please take a look at the SPEC failure here?

smiet · 2024-06-27T19:12:47Z

I had forgotten to include the python 3.8 exception to the numpy 2.0 build-time dependency in the python wrapper. Just merged the commit into SPEC master right before you triggered the tests.

Your commit just now should trigger the tests which will pull the update and they should now pass

smiet · 2024-07-03T08:46:01Z

@landreman, @mbkumar, could this PR be moved forward?

Because the discussion is mixed in with a lot of changed code, I am summarizing the main open threads:

NormalField with a reference to the surface

@landreman noted:

It also feels awkward for NormalField to carry around computational_boundary.

I believe this is necessary for the user, as just the Fourier components are meaningless without the surface. I have also added methods as NormalField.get_real_space_field() that performs the Fourier transform (using the surface properties) and is useful for plotting and analyzing the real space field (and passing to a Focus-like algorithm if you want to generate coils for your Spec). My insistence on this is also because of future development in PR3, where I will introduce a CoilNormalField that inherits from it, and requires the Surface. This way the two classes will have more attributes in common.

Location of Fourier-transforming routine

@landreman noted:

Having the new Fourier transform and inverse transform functions go in the base Surface class feels awkward to me. These functions need mpol and ntor, which not all Surface classes have (e.g. SurfaceHenneberg, SurfaceGarabedian). Also the transform and inverse-transform functions don't use any of the surface geometry, just the grid points. So I would suggest moving those two functions to simsopt.util (probably not in a class, just standalone functions).

I agree that having the transform in the Surface base class is not optimal, but I do think it fits well in the SurfaceRZFourier, as it uses the same VMEC Fourier convention. As a user who has a SurfaceRZFourier, and evaluates some quantity on it, it is only natural for the surface to provide the method to transform it into the same representation as it uses itself. Having a separate function in util that you pass both a Surface and the quantity you want transformed seems unnatural to me.
I have moved the function to the SurfaceRZFourier class
If you still prefer util, let me know and I will make the change

implementation of Fourier-transforming routine

@andrewgiuliani noted:

one suggestion to simplify this function is to use a prewritten FFT and inverse FFT algorithm here

I tried this, but failed, the current code works and is clear (though not optimally efficient), hope it is OK.

landreman · 2024-07-03T15:39:16Z

Presently, the fourier_transform and inverse_fourier_transform methods exist twice each, in both Surface and SurfaceRZFourier. Had you intended to remove them from Surface? Having these functions only in SurfaceRZFourier would be ok.

andrewgiuliani

Hi Chris, there are lines missing from code coverage and a merge conflict that needs to be resolved too.

smiet · 2024-07-05T12:27:25Z

I have covered the new and functionally significant added lines of code. The patch has increased coverage over the old state of the SPEC code (though still under target, sorry). Uncovered are mainly raise statements and specific cases regarding the running of nonstellaratorsymmetric equilibria.

I hope to cover those in the future, but it will take time to generate a suitable nonstellaratorsymmetric spec file that runs in a reasonable time (for now we only test readin, the nonstellaratorsymmetric case takes cpu-hours to run and is unsitbable).

With the increased coverage, passing tests, and removed duplicate transform method, I hope we are good to merge!

Thank you all for your input and suggestions and time reviewing.

andrewgiuliani · 2024-07-05T13:52:46Z

Here's an example unit test for exceptions, it's not a massive change and it makes sure no lines of code in production are untested:

simsopt/tests/geo/test_curve.py

Line 166 in 2c0f057

def test_curvexyzsymmetries_raisesexception(self):

thanks!

smiet · 2024-07-05T15:28:20Z

I added the tests to the raises in NormalField. I will be doing some further additions during my next PR's is it OK if I catch the straggling uncovered exceptions in my next PR?
(the push of these extra tests undid the approval, sorry @landreman)

smiet added 5 commits May 23, 2024 17:13

re-factor normal_field and spec, add FT fn to surface

e076915

fix linting

4983382

fix assertion logic in surface FT

06cfbce

test array translator and default_freeboundary spec

9c2db06

Merge remote-tracking branch 'origin/master' into cbs/spec_update_rebase

2486624

smiet requested review from mbkumar, landreman and abaillod May 30, 2024 13:53

landreman reviewed Jun 2, 2024

View reviewed changes

andrewgiuliani reviewed Jun 4, 2024

View reviewed changes

src/simsopt/geo/surface.py Outdated Show resolved Hide resolved

smiet added 8 commits June 12, 2024 16:57

updated test coverage and getter-setters in NormalField

504af28

make Surface's get_quadpoints staticmethods

59877c1

test transformation to real space

9e4ab7d

move fourier transformation to surfaceRZFourier class

5b440ab

change name field to scalar in fourier transform

1396746

add tests and Matt's improvements

9b3f626

read intial guess based on Linitialize not array presence (f90wrap pe…

9f28499

…rsists in tests)

FIX SPORADIC SPEC FAIL

b177ce8

smiet mentioned this pull request Jun 14, 2024

SPEC test_set_profile_* sporadically failing #357

Closed

smiet added 2 commits June 14, 2024 14:50

fix lint and expand test

a45278b

test non-stellarator-symmetric-dof-setters

cac8b20

Matt Landreman added 2 commits June 23, 2024 20:09

In CIs and containers, use the branch of vmec2000 that works with num…

3718a33

…py 2.0

use the branch of vmec2000 that works with numpy 2.0

53bf9fc

landreman reviewed Jun 24, 2024

View reviewed changes

smiet commented Jun 24, 2024

View reviewed changes

smiet added 3 commits June 24, 2024 17:23

imports to top where they belong

acf037d

test_surface copy mistake fix

b6e1fbd

Merge remote-tracking branch 'origin/numpytypes' into ml/vmec_update_…

03b36f1

…for_numpy_2

smiet added 9 commits June 25, 2024 13:46

second time updating pyproject is the charm

2f139bc

test python3.8 support with oldest-supported-numpy

efebae1

replace deprecated np.NINF with -np.inf

ec8af3d

np.product -> np.prod

75a1ddf

Merge branch 'master' into cbs/spec_update_rebase

65b6148

remove deprecated alltrue->all

7c2a1f9

Merge remote-tracking branch 'origin/ml/vmec_update_for_numpy_2' into…

a193ade

… cbs/spec_update_rebase

FIX SPEC SPORADIC FAIL FINALLY

7fb1729

Revert "Merge remote-tracking branch 'origin/ml/vmec_update_for_numpy…

c68c104

…_2' into cbs/spec_update_rebase" This reverts commit a193ade, reversing changes made to 7c2a1f9.

andrewgiuliani previously requested changes Jul 5, 2024

View reviewed changes

smiet added 5 commits July 5, 2024 09:25

Merge branch 'master' into cbs/spec_update_rebase

78e0777

remove vestigial surface fourier transforms

03b64b5

change last occurences of transform_field->transform_scalar

9d807fa

cover normalization of ft surface_rzfourier

ac81224

test normal field setter too

800060f

landreman self-requested a review July 5, 2024 14:08

landreman previously approved these changes Jul 5, 2024

View reviewed changes

Test raises in NormalField

a2a05ae

smiet dismissed landreman’s stale review via a2a05ae July 5, 2024 15:24

landreman self-requested a review July 5, 2024 17:02

landreman approved these changes Jul 5, 2024

View reviewed changes

smiet merged commit f6a8412 into master Jul 8, 2024
47 checks passed

		@@ -839,6 +839,130 @@ def interpolate_on_arclength_grid(self, function, theta_evaluate):

		return function_interpolated

		def fourier_transform_field(self, field, mpol=None, ntor=None, normalization=None, **kwargs):

Spec refactoring/improvement (PR 1 of 3 for single-stage) #418

Spec refactoring/improvement (PR 1 of 3 for single-stage) #418

Conversation

smiet commented May 30, 2024

codecov bot commented May 30, 2024 • edited Loading

Codecov Report

landreman left a comment

Choose a reason for hiding this comment

landreman Jun 2, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

smiet commented Jun 17, 2024

landreman left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

smiet left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mbkumar commented Jun 27, 2024

smiet commented Jun 27, 2024

smiet commented Jul 3, 2024 • edited Loading

NormalField with a reference to the surface

Location of Fourier-transforming routine

implementation of Fourier-transforming routine

landreman commented Jul 3, 2024 • edited Loading

andrewgiuliani left a comment

Choose a reason for hiding this comment

smiet commented Jul 5, 2024

andrewgiuliani commented Jul 5, 2024

smiet commented Jul 5, 2024

codecov bot commented May 30, 2024 •

edited

Loading

landreman Jun 2, 2024 •

edited

Loading

smiet left a comment •

edited

Loading

smiet commented Jul 3, 2024 •

edited

Loading

landreman commented Jul 3, 2024 •

edited

Loading