Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions process/core/io/plot_proc.py
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,11 @@
from process.models.physics.impurity_radiation import read_impurity_file
from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel
from process.models.physics.physics import BetaComponentLimits, BetaNormMaxModel
from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel
from process.models.physics.plasma_current import (
PlasmaCurrent,
PlasmaCurrentModel,
PlasmaDiamagneticCurrentModel,
)
from process.models.superconductors import SuperconductorModel
from process.models.tfcoil.base import TFCoilShapeModel

Expand Down Expand Up @@ -3080,9 +3084,9 @@ def plot_main_plasma_information(
# Add plasma current information
textstr_currents = (
f"$\\mathbf{{Plasma\\ currents:}}$\n\n"
f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA\n"
f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA \n"
Copy link

Copilot AI Apr 7, 2026

Choose a reason for hiding this comment

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

There are extra trailing spaces before the newline in the plasma current line (... MA \n). This can produce odd spacing/wrapping in the rendered matplotlib text and is inconsistent with the other lines. Remove the trailing spaces unless they are required for a specific alignment behavior.

Suggested change
f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA \n"
f"Plasma current ({PlasmaCurrentModel(int(mfile.get('i_plasma_current', scan=scan))).full_name}): {mfile.get('plasma_current_ma', scan=scan):.4f} MA\n"

Copilot uses AI. Check for mistakes.
f" - Bootstrap fraction {mfile.get('f_c_plasma_bootstrap', scan=scan):.4f}\n"
f" - Diamagnetic fraction {mfile.get('f_c_plasma_diamagnetic', scan=scan):.4f}\n"
f" - Diamagnetic fraction ({PlasmaDiamagneticCurrentModel(int(mfile.get('i_diamagnetic_current', scan=scan))).full_name}): {mfile.get('f_c_plasma_diamagnetic', scan=scan):.4f}\n"
f" - Pfirsch-Schlüter fraction {mfile.get('f_c_plasma_pfirsch_schluter', scan=scan):.4f}\n"
f" - Auxiliary fraction {mfile.get('f_c_plasma_auxiliary', scan=scan):.4f}\n"
f" - Inductive fraction {mfile.get('f_c_plasma_inductive', scan=scan):.4f}"
Expand All @@ -3105,7 +3109,7 @@ def plot_main_plasma_information(

# Add plasma current label
axis.text(
0.925,
0.93,
0.9,
"$I_{\\text{p}} $",
fontsize=23,
Expand Down
4 changes: 3 additions & 1 deletion process/main.py
Original file line number Diff line number Diff line change
Expand Up @@ -105,7 +105,7 @@
PlasmaBeta,
PlasmaInductance,
)
from process.models.physics.plasma_current import PlasmaCurrent
from process.models.physics.plasma_current import PlasmaCurrent, PlasmaDiamagneticCurrent
from process.models.physics.plasma_fields import PlasmaFields
from process.models.physics.plasma_geometry import PlasmaGeom
from process.models.physics.plasma_profiles import PlasmaProfile
Expand Down Expand Up @@ -721,6 +721,7 @@ def __init__(self, data: DataStructure):
self.plasma_transition = PlasmaConfinementTransition()
self.plasma_current = PlasmaCurrent()
self.plasma_fields = PlasmaFields()
self.plasma_dia_current = PlasmaDiamagneticCurrent()
self.physics = Physics(
plasma_profile=self.plasma_profile,
current_drive=self.current_drive,
Expand All @@ -733,6 +734,7 @@ def __init__(self, data: DataStructure):
plasma_transition=self.plasma_transition,
plasma_current=self.plasma_current,
plasma_fields=self.plasma_fields,
plasma_dia_current=self.plasma_dia_current,
)
self.physics_detailed = DetailedPhysics(
plasma_profile=self.plasma_profile,
Expand Down
14 changes: 0 additions & 14 deletions process/models/physics/bootstrap_current.py
Original file line number Diff line number Diff line change
Expand Up @@ -1307,20 +1307,6 @@ def output(self):
"OP ",
)

po.ovarrf(
self.outfile,
"Diamagnetic fraction (Hender)",
"(f_c_plasma_diamagnetic_hender)",
current_drive_variables.f_c_plasma_diamagnetic_hender,
"OP ",
)
po.ovarrf(
self.outfile,
"Diamagnetic fraction (SCENE)",
"(f_c_plasma_diamagnetic_scene)",
current_drive_variables.f_c_plasma_diamagnetic_scene,
"OP ",
)
po.ovarrf(
self.outfile,
"Pfirsch-Schlueter fraction (SCENE)",
Expand Down
67 changes: 5 additions & 62 deletions process/models/physics/physics.py
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@
from process.models.physics.density_limit import PlasmaDensityLimit
from process.models.physics.exhaust import PlasmaExhaust
from process.models.physics.l_h_transition import PlasmaConfinementTransition
from process.models.physics.plasma_current import PlasmaCurrent
from process.models.physics.plasma_current import PlasmaCurrent, PlasmaDiamagneticCurrent
from process.models.physics.plasma_fields import PlasmaFields

logger = logging.getLogger(__name__)
Expand Down Expand Up @@ -145,45 +145,6 @@ def rether(
# -----------------------------------------------------


@nb.jit(nopython=True, cache=True)
def diamagnetic_fraction_hender(beta: float) -> float:
"""Calculate the diamagnetic fraction based on the Hender fit.

Parameters
----------
beta :
the plasma beta value

Returns
-------
:
the diamagnetic fraction

"""
return beta / 2.8


@nb.jit(nopython=True, cache=True)
def diamagnetic_fraction_scene(beta: float, q95: float, q0: float) -> float:
"""Calculate the diamagnetic fraction based on the SCENE fit by Tim Hender.

Parameters
----------
beta :
the plasma beta value
q95 :
the normalized safety factor at 95% of the plasma radius
q0 :
the normalized safety factor at the magnetic axis

Returns
-------
:
the diamagnetic fraction
"""
return beta * (0.1 * q95 / q0 + 0.44) * 0.414


@nb.jit(nopython=True, cache=True)
def ps_fraction_scene(beta: float) -> float:
"""Calculate the Pfirsch-Schlüter fraction based on the SCENE fit by Tim Hender 2019.
Expand Down Expand Up @@ -217,6 +178,7 @@ def __init__(
plasma_transition: PlasmaConfinementTransition,
plasma_current: PlasmaCurrent,
plasma_fields: PlasmaFields,
plasma_dia_current: PlasmaDiamagneticCurrent,
):
self.outfile = constants.NOUT
self.mfile = constants.MFILE
Expand All @@ -231,6 +193,7 @@ def __init__(
self.plasma_transition = plasma_transition
self.current = plasma_current
self.fields = plasma_fields
self.dia_current = plasma_dia_current

def output(self):
self.calculate_effective_charge_ionisation_profiles()
Expand Down Expand Up @@ -501,28 +464,7 @@ def run(self):
# DIAMAGNETIC CURRENT #
# ***************************** #

# Hender scaling for diamagnetic current at tight physics_variables.aspect ratio
current_drive_variables.f_c_plasma_diamagnetic_hender = (
diamagnetic_fraction_hender(physics_variables.beta_total_vol_avg)
)

# SCENE scaling for diamagnetic current
current_drive_variables.f_c_plasma_diamagnetic_scene = (
diamagnetic_fraction_scene(
physics_variables.beta_total_vol_avg,
physics_variables.q95,
physics_variables.q0,
)
)

if physics_variables.i_diamagnetic_current == 1:
current_drive_variables.f_c_plasma_diamagnetic = (
current_drive_variables.f_c_plasma_diamagnetic_hender
)
elif physics_variables.i_diamagnetic_current == 2:
current_drive_variables.f_c_plasma_diamagnetic = (
current_drive_variables.f_c_plasma_diamagnetic_scene
)
self.dia_current.run()

# ***************************** #
# PFIRSCH-SCHLÜTER CURRENT #
Expand Down Expand Up @@ -3212,6 +3154,7 @@ def outplas(self):
self.inductance.output_volt_second_information()
if stellarator_variables.istell == 0:
self.plasma_bootstrap_current.output()
self.dia_current.output()

po.osubhd(self.outfile, "Fuelling :")
po.ovarre(
Expand Down
123 changes: 122 additions & 1 deletion process/models/physics/plasma_current.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,18 @@
import logging
from enum import IntEnum
from types import DynamicClassAttribute

import matplotlib.pyplot as plt
import numba as nb
import numpy as np

import process.core.io.mfile as mf
from process.core import constants
from process.core import process_output as po
from process.core.exceptions import ProcessValueError
from process.core.model import Model
from process.data_structure import (
current_drive_variables,
physics_variables,
)

Expand All @@ -29,9 +33,13 @@ class PlasmaCurrentModel(IntEnum):
def __new__(cls, value, full_name):
obj = int.__new__(cls, value)
obj._value_ = value
obj.full_name = full_name
obj._full_name_ = full_name
return obj

@DynamicClassAttribute
def full_name(self):
return self._full_name_


class PlasmaCurrent:
"""Class to hold plasma current calculations for plasma processing."""
Expand Down Expand Up @@ -925,3 +933,116 @@ def calculate_current_coefficient_fiesta(
Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530.
"""
return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060


class PlasmaDiamagneticCurrentModel(IntEnum):
"""Enum for plasma diamagnetic current method types"""

NONE = (0, "None")
HENDER_ST_FIT = (1, "Hender ST fit")
SCENE_FIT = (2, "SCENE fit")

def __new__(cls, value, full_name):
obj = int.__new__(cls, value)
obj._value_ = value
obj._full_name_ = full_name
return obj

@DynamicClassAttribute
def full_name(self):
return self._full_name_


class PlasmaDiamagneticCurrent(Model):
"""Class to hold plasma diamagnetic current calculations for plasma processing."""

def __init__(self):
self.outfile = constants.NOUT
self.mfile = constants.MFILE

def run(self):
# Hender scaling for diamagnetic current at tight physics_variables.aspect ratio
current_drive_variables.f_c_plasma_diamagnetic_hender = (
self.diamagnetic_fraction_hender(physics_variables.beta_total_vol_avg)
)

# SCENE scaling for diamagnetic current
current_drive_variables.f_c_plasma_diamagnetic_scene = (
self.diamagnetic_fraction_scene(
physics_variables.beta_total_vol_avg,
physics_variables.q95,
physics_variables.q0,
)
)

if (
physics_variables.i_diamagnetic_current
== PlasmaDiamagneticCurrentModel.HENDER_ST_FIT
):
current_drive_variables.f_c_plasma_diamagnetic = (
current_drive_variables.f_c_plasma_diamagnetic_hender
)
elif (
physics_variables.i_diamagnetic_current
== PlasmaDiamagneticCurrentModel.SCENE_FIT
):
current_drive_variables.f_c_plasma_diamagnetic = (
current_drive_variables.f_c_plasma_diamagnetic_scene
)

Comment on lines +978 to +992
Copy link

Copilot AI Apr 7, 2026

Choose a reason for hiding this comment

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

PlasmaDiamagneticCurrent.run() never assigns current_drive_variables.f_c_plasma_diamagnetic when i_diamagnetic_current is NONE (0) (or any unexpected value). In iterative/convergence runs this can leave a stale nonzero fraction from a previous iteration, affecting f_c_plasma_internal and downstream current balance. Set f_c_plasma_diamagnetic = 0.0 for NONE (and ideally handle/raise on unknown values) so the result is deterministic each run.

Copilot uses AI. Check for mistakes.
Comment on lines +963 to +992
Copy link

Copilot AI Apr 7, 2026

Choose a reason for hiding this comment

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

PlasmaDiamagneticCurrent.run() introduces selection logic based on physics_variables.i_diamagnetic_current, but there is no unit test asserting that f_c_plasma_diamagnetic is correctly set for each mode (including NONE/0 and that it does not retain stale values across runs). Adding a small test around run() would help prevent regressions in current fraction accounting.

Copilot uses AI. Check for mistakes.
def output(self):
po.oblnkl(self.outfile)
po.ovarrf(
self.outfile,
"Diamagnetic fraction (Hender)",
"(f_c_plasma_diamagnetic_hender)",
current_drive_variables.f_c_plasma_diamagnetic_hender,
"OP ",
)
po.ovarrf(
self.outfile,
"Diamagnetic fraction (SCENE)",
"(f_c_plasma_diamagnetic_scene)",
current_drive_variables.f_c_plasma_diamagnetic_scene,
"OP ",
)
po.oblnkl(self.outfile)

@staticmethod
@nb.njit(cache=True)
def diamagnetic_fraction_hender(beta: float) -> float:
"""Calculate the diamagnetic fraction based on the Hender fit.

Parameters
----------
beta :
the plasma beta value

Returns
-------
float
the diamagnetic fraction

"""
return beta / 2.8

@staticmethod
@nb.njit(cache=True)
def diamagnetic_fraction_scene(beta: float, q95: float, q0: float) -> float:
"""Calculate the diamagnetic fraction based on the SCENE fit by Tim Hender.
Comment on lines +1011 to +1032
Copy link

Copilot AI Apr 7, 2026

Choose a reason for hiding this comment

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

The new diamagnetic fraction helpers use @nb.jit(cache=True) without nopython=True/@nb.njit. With Numba defaults this can silently fall back to object mode, losing most of the performance benefit and potentially changing behavior between environments. Since these functions are simple numeric expressions, consider using @nb.njit(cache=True) (or @nb.jit(nopython=True, cache=True)) to guarantee nopython compilation.

Copilot uses AI. Check for mistakes.

Parameters
----------
beta :
the plasma beta value
q95 :
the normalized safety factor at 95% of the plasma radius
q0 :
the normalized safety factor at the magnetic axis

Returns
-------
float
the diamagnetic fraction
"""
return beta * (0.1 * q95 / q0 + 0.44) * 0.414
13 changes: 6 additions & 7 deletions tests/unit/test_physics.py
Original file line number Diff line number Diff line change
Expand Up @@ -31,13 +31,11 @@
PlasmaBeta,
PlasmaInductance,
calculate_cylindrical_safety_factor,
diamagnetic_fraction_hender,
diamagnetic_fraction_scene,
ps_fraction_scene,
res_diff_time,
rether,
)
from process.models.physics.plasma_current import PlasmaCurrent
from process.models.physics.plasma_current import PlasmaCurrent, PlasmaDiamagneticCurrent
from process.models.physics.plasma_fields import PlasmaFields
from process.models.physics.plasma_profiles import PlasmaProfile

Expand Down Expand Up @@ -68,6 +66,7 @@ def physics():
PlasmaConfinementTransition(),
PlasmaCurrent(),
PlasmaFields(),
plasma_dia_current=PlasmaDiamagneticCurrent(),
)


Expand All @@ -83,19 +82,19 @@ def test_res_diff_time():
assert t_plasma_res_diffusion == pytest.approx(4784.3, abs=0.1)


def test_diamagnetic_fraction_hender():
def test_diamagnetic_fraction_hender(physics):
"""Test diamagnetic_fraction_hender()."""
beta = 0.14
diacf = diamagnetic_fraction_hender(beta)
diacf = physics.dia_current.diamagnetic_fraction_hender(beta)
assert diacf == pytest.approx(0.05, abs=0.0001)


def test_diamagnetic_fraction_scene():
def test_diamagnetic_fraction_scene(physics):
"""Test diamagnetic_fraction_scene."""
beta = 0.15
q95 = 3.0
q0 = 1.0
diacf = diamagnetic_fraction_scene(beta, q95, q0)
diacf = physics.dia_current.diamagnetic_fraction_scene(beta, q95, q0)
assert diacf == pytest.approx(0.0460, abs=0.0001)


Expand Down
Loading
Loading