Integrate veqpy as an optional tokamak equilibrium solver

[chencx04]

Motivation

PROCESS is a tokamak system code designed for fast system-level modelling and optimisation. Many equilibrium-related quantities used by a system code, however, are currently obtained from analytical models, parameterisations, or prescribed profiles rather than from a self-consistent Grad–Shafranov equilibrium calculation.

This can introduce inconsistencies between the assumed plasma shape, pressure/current profiles, safety-factor profile, and other equilibrium quantities. A fast equilibrium solver could provide a more self-consistent description of the plasma equilibrium while retaining the computational efficiency required for repeated evaluations in a system-code optimisation loop.

I propose integrating veqpy as an optional equilibrium solver in PROCESS.

[veqpy](https://github.com/FusionAlpha/veqpy) is a Python implementation of VEQ (Veloce EQuilibrium), providing a fast parametric Grad–Shafranov solver for fixed-boundary, axisymmetric tokamak equilibria. It represents the equilibrium using flux-surface harmonics and radial profile/source coefficients and provides a continuous equilibrium representation that can be resampled and diagnosed. The package is designed for repeated modelling calls, parameter scans, transport coupling, and control-oriented applications.

I am willing to implement this feature and submit a pull request.

Proposed functionality

The initial goal would be to make veqpy available as an optional equilibrium solver within PROCESS, without changing the existing default workflow.

Conceptually, the workflow would be:

PROCESS plasma/system parameters
              |
              v
       Equilibrium interface
              |
       +------+------+
       |             |
       v             v
Existing approach   veqpy
       |             |
       +------+------+
              |
              v
   Equilibrium quantities
              |
              v
      PROCESS physics

When veqpy is not selected, the existing PROCESS calculations should remain unchanged.

When veqpy is selected, PROCESS would construct the required equilibrium inputs, call the veqpy solver, and make the resulting equilibrium quantities available to downstream physics models.

The exact interface and location within the PROCESS architecture should be discussed with the PROCESS developers before implementation.

Potential equilibrium quantities

The initial integration could expose a limited set of quantities that are useful to PROCESS, such as:

• poloidal flux and normalized poloidal flux;
• pressure profiles;
• toroidal and/or poloidal current profiles;
• safety-factor profile q(rho);
• q95;
• plasma current distribution;
• magnetic-field quantities;
• flux-surface averaged quantities;
• equilibrium geometry and derived quantities.

The initial implementation would not necessarily replace all existing equilibrium-related calculations. Instead, it would establish a clean interface through which equilibrium quantities can be provided to PROCESS.

Why veqpy

veqpy is particularly interesting for PROCESS because it is designed for repeated equilibrium calculations rather than only for single high-fidelity equilibrium reconstruction.

The package provides:

• a compact parametric representation of fixed-boundary equilibria;
• a public Kernel interface for equilibrium solves;
• pressure/current/safety-factor based source routes;
• continuous equilibrium snapshots;
• GEQDSK input/output support;
• Numba and optional C++ backends;
• support for repeated modelling and parameter-scan workflows.

The public API separates equilibrium topology, runtime boundary/source data, solver configuration, and the resulting equilibrium representation. This may provide a suitable basis for integrating an equilibrium solver without tightly coupling PROCESS to the internal implementation of veqpy.

Proposed integration strategy

I suggest initially treating veqpy as an external Python dependency rather than copying its implementation into PROCESS.

A possible high-level interface could be:

equilibrium = equilibrium_solver.solve(
    boundary=boundary,
    profiles=profiles,
    plasma_current=plasma_current,
    ...
)

The exact interface would be determined after discussion with the PROCESS maintainers.

The integration should avoid depending unnecessarily on veqpy's internal implementation details. PROCESS should interact with a small, well-defined equilibrium interface and obtain only the quantities required by its physics models.

Initial scope

The initial implementation would focus on:

1. Adding an optional veqpy equilibrium-solver dependency/interface.
2. Translating PROCESS equilibrium inputs into the corresponding veqpy inputs.
3. Calling veqpy through its public API.
4. Extracting a defined set of equilibrium quantities.
5. Making these quantities available to PROCESS physics modules.
6. Providing an option to enable/disable veqpy.
7. Adding unit and integration tests.
8. Adding documentation and an example demonstrating the workflow.

The initial implementation would not attempt to replace every existing equilibrium-related model in PROCESS.

Further integration with individual physics models could be addressed in follow-up issues once the basic equilibrium interface has been established.

Validation

The integrated equilibrium results should be validated against an independent Grad–Shafranov equilibrium calculation.

Possible validation quantities include:

• plasma boundary;
• magnetic flux surfaces;
• normalized poloidal flux;
• pressure profile;
• current profile;
• safety-factor profile;
• q95;
• plasma current;
• relevant flux-surface averaged quantities.

Where appropriate, GEQDSK-based cases could also be used for validation because veqpy already supports GEQDSK workflows.

The initial validation should include representative tokamak configurations covering different plasma shaping and profile cases.

PROCESS performance

Because equilibrium calculations may be called repeatedly during PROCESS optimisation, computational cost is an important consideration.

The integration should therefore evaluate:

• equilibrium solve time;
• overhead associated with PROCESS-to-veqpy data conversion;
• repeated-call performance;
• behaviour when the equilibrium solver is called many times during optimisation.

If appropriate, warm-start or continuation capabilities provided by veqpy could be considered in a later stage to reduce the cost of repeated equilibrium calculations.

Backward compatibility

The existing PROCESS workflow should remain available when veqpy is disabled.

In particular:

• existing input files should continue to work;
• existing equilibrium-related calculations should remain available;
• veqpy should not become a mandatory dependency unless agreed by the PROCESS maintainers;
• enabling veqpy should be explicit;
• unrelated PROCESS calculations should not be changed by the initial integration.

Testing

The implementation would include appropriate tests for the new functionality, including:

Unit tests

• PROCESS-to-veqpy input conversion;
• equilibrium solver interface;
• extraction and conversion of equilibrium quantities;
• handling of invalid or unsupported equilibrium inputs.

Integration tests

A representative PROCESS case would be used to verify the complete workflow:

PROCESS
   |
   v
veqpy equilibrium
   |
   v
equilibrium quantities
   |
   v
PROCESS physics calculations

Regression tests

Existing PROCESS regression cases would be evaluated to ensure that the new optional functionality does not unintentionally alter the existing default behaviour.

Any intentional regression changes resulting from using veqpy would be documented and justified in the corresponding pull request.

Documentation

The implementation should include documentation covering:

• the purpose of the veqpy equilibrium option;
• installation/dependency requirements;
• relevant PROCESS input parameters;
• the equilibrium quantities provided by veqpy;
• an example PROCESS case;
• known limitations and assumptions;
• validation results.

Questions for discussion

Before implementation, I would appreciate feedback from the PROCESS developers on the following points:

1. Would an optional veqpy-based equilibrium solver be compatible with the current PROCESS architecture and development direction?
2. Where would an equilibrium-solver interface be most appropriately located in the current PROCESS codebase?
3. Which equilibrium quantities should be included in the initial implementation?
4. Should veqpy initially be implemented as an alternative equilibrium calculation, or should it first provide equilibrium information to selected existing physics models?
5. Should veqpy be an optional dependency, or would including it as a standard PROCESS dependency be preferable?
6. Are there existing or planned equilibrium-solver integrations that should be considered before proceeding?

Expected outcome

The intended outcome is an optional veqpy-based equilibrium capability in PROCESS that can provide self-consistent equilibrium quantities while preserving the existing PROCESS workflow and computational efficiency.

This would establish a foundation for future equilibrium-dependent physics models in PROCESS, including more self-consistent calculations of current profiles, safety-factor profiles, bootstrap current, MHD stability quantities, transport-related quantities, and other equilibrium-dependent physics.

I propose we discuss this feature with the PROCESS developers, and I am happy to implement it and create a pull request if the approach is considered suitable.

[grmtrkngtn]

Thank you for opening this discussion. I have read some of the validation in your paper. I was wondering how well the reduced VEQ representation performs when moving away from reconstruction of a known reference equilibrium and into the wider parameter studies you would encounter in a systems-code optimisation.

In PROCESS, the optimiser may vary quantities such as aspect ratio, shaping, beta quite substantially between evaluations, and some trial points may be far from a well-behaved reference equilibrium. Have you tested the robustness and accuracy of VEQ in that kind of situation?

[chencx04]

Thank you for raising this important point. We agree that the robustness of the reduced representation over a broad parameter space is particularly important for a systems-code application, where the optimiser may explore configurations substantially different from a single reference equilibrium.

For the earlier VEQ model, in which three parameters were used as the variational parameters, we performed systematic one-at-a-time parameter scans to examine its behaviour over a relatively broad range of geometric and physical inputs. In these scans, all other parameters were fixed at the EHL-2 reference configuration while one parameter was varied. The scanned parameters included aspect ratio, elongation, triangularity, plasma current, vacuum magnetic field, density/temperature/current peaking factors, and the on-axis ion density and temperature.

For example, the parameter ranges included:

• aspect ratio: 1.2–1.85;
• elongation: 1.0–2.5;
• triangularity: 0.1–0.8;
• plasma current: 3–6 MA;
• vacuum magnetic field: 3–5 T;
• density peaking factor: 1.0–2.5;
• temperature peaking factor: 1.0–3.0;
• current peaking factor: 0.1–3.0;
• on-axis ion density: (1.0–1.8\times10^{20},\mathrm{m^{-3}});
• on-axis ion temperature: 20–60 keV.

Within these scans, the VEQ model generally produced good equilibrium solutions and showed stable behaviour over the tested parameter ranges. We also generated nearly 1500 equilibrium cases for the q-profile analysis, with more than 1000 cases achieving a q-profile fitting error below 10%.

There is, however, an important distinction between solver robustness and the existence of a physical equilibrium. For example, when the elongation was fixed at 2.2, no physically acceptable Grad–Shafranov equilibrium could be obtained once the aspect ratio decreased to approximately 1.55 or below. We therefore reduced the elongation to 1.2 for the corresponding scan to ensure that the parameter space contained physically admissible equilibria. This also illustrates that some apparently problematic optimisation points may correspond to configurations for which no physical equilibrium exists, rather than a failure of the VEQ representation itself.

The current veqpy implementation is a more general and more complete implementation of the VEQ approach, with a flexible source representation and a unified API. We therefore expect it to provide at least comparable, and potentially better, robustness and applicability for broad parameter scans. However, we agree that this should be demonstrated quantitatively rather than assumed.

For the PROCESS integration, we think this robustness test can be carried out after veqpy has been integrated into PROCESS, together with the integration tests. This would allow us to evaluate veqpy directly under the conditions encountered by the systems code, including variations of the relevant optimisation parameters, and to assess both the convergence of the equilibrium solver and the accuracy of key equilibrium quantities against the corresponding reference results. We would be happy to include such tests as part of the integration and validation work.