Skip to content

add GVEC as MHD optimizable#481

Draft
Rykath wants to merge 39 commits intohiddenSymmetries:masterfrom
gvec-group:gvec
Draft

add GVEC as MHD optimizable#481
Rykath wants to merge 39 commits intohiddenSymmetries:masterfrom
gvec-group:gvec

Conversation

@Rykath
Copy link

@Rykath Rykath commented Mar 18, 2025

This PR adds an Optimizable-class for the Galerkin Variational Equilibrium Code (GVEC). GVEC is a flexible 3D ideal MHD equilibrium solver inspired by VMEC with a few notable differences. The two most prominent are:

  • The radial direction is discretized using B-Splines rather than finite differences on a fixed grid. This allows representing the region of the magnetic axis accurately and evaluations anywhere in the domain.
  • The coordinate mapping from flux aligned coordinates to real space is flexible and not fixed to $(R,Z,\phi)$. This also allows using a flexible curve-following coordinate frame, e.g. to compute equilibria for a figure-8 stellarator

GVEC is developed at the division of Numerical Methods of the Max Planck Institute for Plasma Physics in Garching, Germany. It has recently been released publicly and open source (MIT License).
Details on the implementation of GVEC are partially described in the documentation.


This PR adds a simsopt.mhd.gvec.Gvec optimizable, similar to the Vmec and Spec optimizables.

  • DoFs: phiedge, boundary, pressure and iota or current profile
  • GVEC is parallelized with OpenMP and MPI. For now, parallel least_squares_mpi_solve is only supported with one MPI process for each GVEC run.
  • fixed boundary
  • prescribing a iota or (toroidal) current profile
  • basic return functions for iota, iota_axis, volume, etc.
  • evaluating other equilibrium quantities (e.g. boozer transform) can use GVEC's python API, the GVECQuantity optimizable can be used to access any of the computable quantities of GVEC.

This PR is still work in progress, but we wanted to show our progress as soon as possible. Open ToDo's:

  • prescribe toroidal current
  • documentation
  • examples / stellarator benchmarks
    • 1DOF_circularCrossSection_varyAxis_targetIota_gvec.py
    • 1DOF_circularCrossSection_varyR0_targetVolume_gvec.py
    • 2DOF_gvecOnly_targetIotaAndVolume.py
  • unit tests
  • add SurfaceGVECFourier to represent surface DoFs when not using cylindrical coordinates
    • this is not a Surface object, e.g. does not support the gamma() method

Let us know what else you need.

Rykath added 3 commits March 6, 2025 16:34
adds a GVEC optimizable class, similar to the VMEC optimizable
* refactor Gvec optimizable to use the run_stages script: allows running GVEC with refinement and current constraint (optimization using picard iterations)
* profile objects are now available as iota_profile, pressure_profiel and current_profile, while iota and I_tor are return functions
* Gvec init using dictionary of parameters, from_parameter_file classmethod additionally
* adapted example parameter file & run script for 1DoF circularCrossSection target Volume
* from_parameters & from_parameter_file classmethods are also available
* allow setting iota & current profile afterwards
* add default parameters (minimization and radial resolution)
* the varyR0_targetVolume example now no longer needs an input file
Copy link
Contributor

@smiet smiet left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looking very nice! Only major issue is with the use of SurfaceRZFourier as the boundary object, which I believe should be replaced with a SurfaceGVEC or SurfaceRZFourierOnFrame lest users try to access, plot, evaluate, and optimize properties of the boundary directly.

Aditionally, unit tests need to be written to test all the new code. See tests/mhd/vmec.py for a good example to adapt.
This will also involve updating the workflow files to install GVEC.

ToDo: doc
"""

def __init__(
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In order for the object to be GSONable (to be saved and then loaded again using simsopt methods), all parameters to the __init__ function must be stored either as self.<param> or self._<param>. If you want this feature, parameters needs to be an attribute of self.

if f"{key}_mn_max" not in params:
params[f"{key}_mn_max"] = (boundary.mpol, boundary.ntor)

if not isinstance(boundary, SurfaceRZFourier):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Transformation to SurfaceRZFourier already done above

self._state.compute(ds, "V")
return ds.V.item()

def aspect(self) -> float:
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Aspect() should also be a property of the boundary Surface object


assert np.abs(surf.get_rc(0, 0) - 0.7599088773175) < 1.0e-5
assert np.abs(equil.volume() - 0.15) < 1.0e-6
assert np.abs(surf.volume() - 0.15) < 1.0e-6
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are evaluating the surfaces' volume, but this is only equal to equil.volume in this case, and is misleading. See other comments on creating a SurfaceGVEC class for representing the surface.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

again only following the VMEC example. But I agree that the current lack of a SurfaceGVEC can be misleading.

Rykath added 16 commits May 8, 2025 17:27
* improved handling of parameters:
  * surface can be specified via the parameters or as an object
  * profiles can be specified via the parameters, an object or a constant
* iota profile is always used, but if current is specified it is only used as the initial guess
* removed from_parameters constructor - all features now available in init
* add logger property to access the GVEC logger
* make paths in the parameterfile absolute when writing parameters for each run
* new argument `keep_gvec_intermediates`: passed to gvec.run with values "all", "stages" or None
* new decorator `return_function` to mark all return functions which are then registered with `return_fn_map`
* cleaned return functions
  * vacuum_well is now the vacuum magnetic well depth with the same  formual as vmec
  * removed current & iota return functions that returned arrays as the resolution is not customizable
  * return functions now use the .state().evaluate(...)
* new function `boundary_to_params` & overhaul of `boundary_from_params`: now supports non-stellarator-symmetry
* fix current optimization by setting `picard_current` correctly
* Gvec.state is now a property and not a return function
* new Gvec.beta_avg return function
* new MirrorRatio optimizable based on gvec
* new Elongation optimizable based on gvec
* new Gvec.from_rundir constructor
* Gvec.run_count now starts at 0
* changed Gvec option restart_file -> restart
* prototype for restart functionality (wip)
@Rykath
Copy link
Author

Rykath commented Dec 2, 2025

I apologize for not updating this in a very long while.
In the meantime we've done quite some work behind the scenes, in particular extending the utilities and improving the stability of GVEC, which was limiting this integration. With GVEC v1.3 I think we have everything ready on the GVEC side, and only the open ToDos here remain.

@Rykath
Copy link
Author

Rykath commented Jan 30, 2026

Only major issue is with the use of SurfaceRZFourier as the boundary object, which I believe should be replaced with a SurfaceGVEC or SurfaceRZFourierOnFrame lest users try to access, plot, evaluate, and optimize properties of the boundary directly.

I have now added a SurfaceGVECFourier object to hold the surface DoFs when not using cylindrical coordinates with GVEC. I have also added from_gvec_parameters and to_gvec_parameters methods to both SurfaceRZFourier and SurfaceGVECFourier (rather than doing that conversion within the gvec optimizable object).

@smiet
Copy link
Contributor

smiet commented Jan 30, 2026

Nice work! This cleans up the dealing with degrees of freedom.

For full integration with simsopt, the SurfaceGvecFourier could be made to inherit from the base Surface class. Then it would need to expose a few more methods such as gamma() (map from the unit square to R3) so it can be plotted, fields evaluated on it, etc, integrated with Coil objects.

For fancy optimization like Rogerio's single-stage, you would want dgamma_by_dcoeff (change of points if dofs are changed) etc, but that is not required for core functionality.

If we do not want to make it a full surface, it should still be its own Optimizable object, but maybe called GVECSurfaceDofs, and live in mhd, as it does not implement any geometry.

I urge to implement as a Surface, because that would allow for so many nice things, but I know this is much more work and do not know how easy it is to expose these functions.

@codecov
Copy link

codecov bot commented Jan 30, 2026

Codecov Report

❌ Patch coverage is 70.91837% with 171 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.66%. Comparing base (4883116) to head (b12d774).
⚠️ Report is 2 commits behind head on master.

Files with missing lines Patch % Lines
src/simsopt/mhd/gvec.py 66.58% 137 Missing ⚠️
src/simsopt/mhd/gvec_dofs.py 81.89% 21 Missing ⚠️
src/simsopt/geo/surfacerzfourier.py 78.68% 13 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #481      +/-   ##
==========================================
- Coverage   90.23%   89.66%   -0.57%     
==========================================
  Files          84       86       +2     
  Lines       17896    18484     +588     
==========================================
+ Hits        16148    16574     +426     
- Misses       1748     1910     +162     
Flag Coverage Δ
unittests 89.66% <70.91%> (-0.57%) ⬇️

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.

@Rykath
Copy link
Author

Rykath commented Feb 2, 2026

For full integration with simsopt, the SurfaceGvecFourier could be made to inherit from the base Surface class. Then it would need to expose a few more methods such as gamma() (map from the unit square to R3) so it can be plotted, fields evaluated on it, etc, integrated with Coil objects.

For fancy optimization like Rogerio's single-stage, you would want dgamma_by_dcoeff (change of points if dofs are changed) etc, but that is not required for core functionality.

I urge to implement as a Surface, because that would allow for so many nice things, but I know this is much more work and do not know how easy it is to expose these functions.

I agree having a proper Surface object would be very nice to have, but exposing these GVEC functions cleanly is not that easy, maybe a re-implementation of the G-Frame (e.g. in JAX) would be faster. It's also not strictly necessary for "normal" stage-1 optimization, so I would really like to postpone this for some later date.

If we do not want to make it a full surface, it should still be its own Optimizable object, but maybe called GVECSurfaceDofs, and live in mhd, as it does not implement any geometry.

Yes, I'll do that.

Rykath added 16 commits February 2, 2026 14:56
* new gvec parameter file for testing: LandremanPaul2021_QA (converted from vmec input)
* added import guard clause for qsc for test_surface_rzfourier
* allows quickly evaluating the root-mean-square
* allows targeting straight-fieldline quantities
was previously moved to mhd.gvec_dofs but accidently and should have been deleted
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants