Integrate veqpy as an optional tokamak equilibrium solver #4593

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.