Add figure of merit enum

documentation/source/development/add-vars.md

@@ -80,9 +80,9 @@ ITERATION_VARIABLES = {
 
 New figures of merit are added to `PROCESS` in the following way:
 
-1. Increment the parameter `ipnfoms` in module `numerics` in source file `numerics.f90` to accommodate the new figure of merit.
+1. Increment the parameter `ipnfoms` in module `numerics` in source file `numerics.py` to accommodate the new figure of merit.
 
-2. Assign a description of the new figure of merit to the relevant element of array `lablmm` in module `numerics` in the source file `numerics.f90`.
+2. Assign the new integer value and description string of the new figure of merit to the `FiguresOfMerit` enumerator in `numerics.py`.
 
 3. Add the new figure of merit equation to `objective_function()` in `objectives.py`, following the method used in the existing examples. The value of figure of merit case should be of order unity, so select a reasonable scaling factor if necessary. 
 

documentation/source/io/input-guide.md

@@ -118,7 +118,7 @@ ioptimz  = 1 * for optimisation VMCON only
 The user can select the figure of merit to be used:
 
 ```
-minmax   = 1 * Switch for figure-of-merit (see lablmm for descriptions)
+minmax   = 1 * Switch for figure-of-merit (see `FiguresOfMerit` for descriptions)
 ```
 
 In this case the user is choosing option `1`, which is major radius. For `minmax`
@@ -132,10 +132,6 @@ The user can also input the allowed error tolerance on the solver solution:
 epsvmc   = 1.0e-8 * Error tolerance for vmcon
 ```
 
-!!! Info "Figure of Merit"  
-    A full list of figures of merit is given on the variable description page in the row labelled 
-    `lablmm` [here](../../source/reference/process/data_structure/numerics/#process.data_structure.numerics.lablmm).
-
 ## Input Variables
 
 One can enter an input into the `IN.DAT` by:

process/core/init.py

@@ -17,6 +17,7 @@
 from process.core.solver import iteration_variables
 from process.core.solver.constraints import ConstraintManager
 from process.data_structure.impurity_radiation_variables import N_IMPURITIES
+from process.data_structure.numerics import FiguresOfMerit
 from process.data_structure.physics_variables import (
     DivertorNumberModels,
     init_physics_module,
@@ -191,9 +192,7 @@ def run_summary():
                 minmax_string = "  -- maximise "
                 minmax_sign = "-"
 
-            fom_string = data_structure.numerics.lablmm[
-                abs(data_structure.numerics.minmax) - 1
-            ]
+            fom_string = FiguresOfMerit(abs(data_structure.numerics.minmax)).description
             process_output.ocmmnt(
                 outfile,
                 f"Figure of merit : {minmax_sign}{abs(data_structure.numerics.minmax)}{minmax_string}{fom_string}",

process/core/io/plot/solutions.py

@@ -22,6 +22,7 @@
 
 from process.core.io.mfile import MFile
 from process.data_structure import numerics
+from process.data_structure.numerics import FiguresOfMerit
 
 # Variables of interest in mfiles and subsequent dataframes
 # Be specific about exact names, patterns and regex
@@ -446,7 +447,7 @@ def _plot_solutions(
     else:
         numerics.init_numerics()
         objf_list = {
-            numerics.lablmm[int(abs(minmax)) - 1] for minmax in diffs_df["minmax"]
+            FiguresOfMerit(abs(int(minmax))).description for minmax in diffs_df["minmax"]
         }
 
     if len(objf_list) != 1:

process/core/io/plot/summary.py

@@ -20,8 +20,8 @@
 
 from process.core import constants
 from process.core.io.mfile import MFile, MFileErrorClass
-from process.core.solver.objectives import OBJECTIVE_NAMES
 from process.data_structure.impurity_radiation_variables import N_IMPURITIES
+from process.data_structure.numerics import FiguresOfMerit
 from process.data_structure.pfcoil_variables import NFIXMX
 from process.models.build import Build
 from process.models.geometry.blanket import (
@@ -7822,7 +7822,7 @@ def plot_header(axis: plt.Axes, mfile: MFile, scan: int):
         ("!Evaluation", "Run type", "")
         if isinstance(mfile.data["minmax"], MFileErrorClass)
         else (
-            f"!{OBJECTIVE_NAMES[abs(int(mfile.get('minmax', scan=-1)))]}",
+            f"!{FiguresOfMerit(abs(int(mfile.get('minmax', scan=-1)))).description}",
             "Optimising:",
             "",
         ),
@@ -11596,10 +11596,10 @@ def plot_cover_page(
         objective_text = ""
     elif minmax_switch >= 0:
         minmax_switch = int(minmax_switch)
-        objective_text = f"• Minimising {objf_name}"
+        objective_text = f"  -> Minimising: {objf_name}"
     else:
         minmax_switch = int(minmax_switch)
-        objective_text = f"• Maximising {objf_name}"
+        objective_text = f"  -> Maximising: {objf_name}"
 
     axis.text(
         0.1,
@@ -11661,9 +11661,9 @@ def plot_cover_page(
     settings_info = (
         f"• Optimisation Switch: {int(optmisation_switch)}\n"
         f"• Figure of Merit Switch (minmax): {minmax_switch}\n"
+        f"     {objective_text}\n"
         f"• Fail Status (ifail): {int(ifail)}\n"
         f"• Number of Iteration Variables: {int(nvars)}\n"
-        f"{objective_text}\n"
         f"• Constraint Residuals (sqrt sum sq): {sqsumsq}\n"
         f"• Convergence Parameter: {convergence_parameter}\n"
         f"• Solver Iterations: {nviter}\n"

process/core/scan.py

@@ -23,6 +23,7 @@
     scan_variables,
     tfcoil_variables,
 )
+from process.data_structure.numerics import FiguresOfMerit
 
 if TYPE_CHECKING:
     from process.core.model import DataStructure, Model
@@ -353,7 +354,7 @@ def post_optimise(self, ifail: int):
                 constants.NOUT, "Figure of merit switch", "(minmax)", numerics.minmax
             )
 
-            objf_name = f'"{numerics.lablmm[abs(numerics.minmax) - 1]}"'
+            objf_name = f'"{FiguresOfMerit(abs(numerics.minmax)).description}"'
 
             numerics.objf_name = objf_name
 

process/core/solver/objectives.py

@@ -6,25 +6,7 @@
     physics_variables,
     tfcoil_variables,
 )
-
-OBJECTIVE_NAMES = {
-    1: "Plasma major radius",
-    3: "neutron wall load",
-    4: "total TF + PF coil power",
-    5: "ratio fusion power:injection power",
-    6: "cost of electricity",
-    7: "constructed cost",
-    8: "aspect ratio",
-    9: "divertor heat load",
-    10: "toroidal field on axis",
-    11: "injection power",
-    14: "pulse length",
-    15: "plant availability factor",
-    16: "Major radius/burn time",
-    17: "net electrical output",
-    18: "NULL",
-    19: "Major radius/burn time",
-}
+from process.data_structure.numerics import FiguresOfMerit
 
 
 def objective_function(minmax: int, data: DataStructure) -> float:
@@ -56,60 +38,56 @@ def objective_function(minmax: int, data: DataStructure) -> float:
         data structure object for providing data to the
         objective function
     """
-    figure_of_merit = abs(minmax)
+    try:
+        figure_of_merit = FiguresOfMerit(abs(minmax))
+    except ValueError as err:
+        raise ProcessValueError(f"Invalid minmax value: {minmax}") from err
 
     # -1 = maximise
     # +1 = minimise
     objective_sign = np.sign(minmax)
 
-    match figure_of_merit:
-        case 1:
-            objective_metric = 0.2 * physics_variables.rmajor
-        case 3:
-            objective_metric = physics_variables.pflux_fw_neutron_mw
-        case 4:
-            objective_metric = (
-                tfcoil_variables.tfcmw + 1e-3 * data.pf_power.srcktpm
-            ) / 10.0
-        case 5:
-            objective_metric = physics_variables.p_fusion_total_mw / (
-                data.current_drive.p_hcd_injected_total_mw
-                + data.current_drive.p_beam_orbit_loss_mw
-                + physics_variables.p_plasma_ohmic_mw
-            )
-        case 6:
-            objective_metric = data.costs.coe / 100.0
-        case 7:
-            objective_metric = (
-                data.costs.cdirt / 1.0e3
-                if data.costs.ireactor == 0
-                else data.costs.concost / 1.0e4
-            )
-        case 8:
-            objective_metric = physics_variables.aspect
-        case 9:
-            objective_metric = data.divertor.pflux_div_heat_load_mw
-        case 10:
-            objective_metric = physics_variables.b_plasma_toroidal_on_axis
-        case 11:
-            objective_metric = data.current_drive.p_hcd_injected_total_mw
-        case 14:
-            objective_metric = data.times.t_plant_pulse_burn / 2.0e4
-        case 15:
-            if data.costs.i_plant_availability != 1:
-                raise ProcessValueError("minmax=15 requires i_plant_availability=1")
-            objective_metric = data.costs.f_t_plant_available
-        case 16:
-            objective_metric = 0.95 * (physics_variables.rmajor / 9.0) - 0.05 * (
-                data.times.t_plant_pulse_burn / 7200.0
-            )
-        case 17:
-            objective_metric = data.heat_transport.p_plant_electric_net_mw / 500.0
-        case 18:
-            objective_metric = 1.0
-        case 19:
-            objective_metric = -0.5 * (data.current_drive.big_q_plasma / 20.0) - 0.5 * (
-                data.times.t_plant_pulse_burn / 7200.0
-            )
+    if figure_of_merit == FiguresOfMerit.MAJOR_RADIUS:
+        objective_metric = 0.2 * physics_variables.rmajor
+    elif figure_of_merit == FiguresOfMerit.NEUTRON_WALL_LOAD:
+        objective_metric = physics_variables.pflux_fw_neutron_mw
+    elif figure_of_merit == FiguresOfMerit.TF_PF_COIL_POWER:
+        objective_metric = (tfcoil_variables.tfcmw + 1e-3 * data.pf_power.srcktpm) / 10.0
+    elif figure_of_merit == FiguresOfMerit.FUSION_GAIN:
+        objective_metric = data.current_drive.big_q_plasma
+    elif figure_of_merit == FiguresOfMerit.COST_OF_ELECTRICITY:
+        objective_metric = data.costs.coe / 100.0
+    elif figure_of_merit == FiguresOfMerit.CAPITAL_COST:
+        objective_metric = (
+            data.costs.cdirt / 1.0e3
+            if data.costs.ireactor == 0
+            else data.costs.concost / 1.0e4
+        )
+    elif figure_of_merit == FiguresOfMerit.ASPECT_RATIO:
+        objective_metric = physics_variables.aspect
+    elif figure_of_merit == FiguresOfMerit.DIVERTOR_HEAT_LOAD:
+        objective_metric = data.divertor.pflux_div_heat_load_mw
+    elif figure_of_merit == FiguresOfMerit.TOROIDAL_FIELD:
+        objective_metric = physics_variables.b_plasma_toroidal_on_axis
+    elif figure_of_merit == FiguresOfMerit.INJECTED_POWER:
+        objective_metric = data.current_drive.p_hcd_injected_total_mw
+    elif figure_of_merit == FiguresOfMerit.PULSE_LENGTH:
+        objective_metric = data.times.t_plant_pulse_burn / 2.0e4
+    elif figure_of_merit == FiguresOfMerit.PLANT_AVAILABILITY:
+        if data.costs.i_plant_availability != 1:
+            raise ProcessValueError("minmax=15 requires i_plant_availability=1")
+        objective_metric = data.costs.f_t_plant_available
+    elif figure_of_merit == FiguresOfMerit.MAJOR_RADIUS_BURN_TIME:
+        objective_metric = 0.95 * (physics_variables.rmajor / 9.0) - 0.05 * (
+            data.times.t_plant_pulse_burn / 7200.0
+        )
+    elif figure_of_merit == FiguresOfMerit.NET_ELECTRICAL_OUTPUT:
+        objective_metric = data.heat_transport.p_plant_electric_net_mw / 500.0
+    elif figure_of_merit == FiguresOfMerit.NULL:
+        objective_metric = 1.0
+    elif figure_of_merit == FiguresOfMerit.FUSION_GAIN_BURN_TIME:
+        objective_metric = -0.5 * (data.current_drive.big_q_plasma / 20.0) - 0.5 * (
+            data.times.t_plant_pulse_burn / 7200.0
+        )
 
     return objective_sign * objective_metric

process/data_structure/numerics.py

@@ -1,5 +1,59 @@
+from enum import IntEnum
+from types import DynamicClassAttribute
+
 import numpy as np
 
+
+class FiguresOfMerit(IntEnum):
+    """Enumeration of the available figures of merit (FoM) that can be used as
+    objective functions for optimisation in PROCESS.
+    """
+
+    MAJOR_RADIUS = (1, "Plasma major radius (R₀)")
+    NEUTRON_WALL_LOAD = (3, "Neutron wall load")
+    P_TF_PLUS_P_PF = (4, "TF & PF coil power")
+    FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)")
+    COST_OF_ELECTRICITY = (6, "Cost of electricity")
+    CAPITAL_COST = (7, "Plant capital cost")
+    ASPECT_RATIO = (8, "Plasma aspect ratio")
+    DIVERTOR_HEAT_LOAD = (9, "Divertor heat load")
+    TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)")
+    TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)")
+    PULSE_LENGTH = (14, "Pulse length")
+    PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor")
+    MIN_R0_MAX_TAU_BURN = (
+        16,
+        "Linear combination of major radius (minimised) and pulse length (maximised)",
+    )
+    NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output")
+    NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit")
+    MAX_Q_MAX_T_PLANT_PULSE_BURN = (
+        19,
+        "Linear combination of big Q and pulse length (maximised)",
+    )
+
+    def __new__(cls, value: int, description: str):
+        """Create a new FiguresOfMerit enum member with description.
+
+        Args:
+            value: The integer value of the enum member.
+            description: The description for this figure of merit.
+
+        Returns
+        -------
+            The new enum member with attached description.
+        """
+        obj = int.__new__(cls, value)
+        obj._value_ = value
+        obj._description_ = description
+        return obj
+
+    @DynamicClassAttribute
+    def description(self):
+        """Return the description for this figure of merit."""
+        return self._description_
+
+
 ipnvars: int = 177
 """total number of variables available for iteration"""
 
@@ -21,37 +75,10 @@
 
 minmax: int = None
 """
-Switch for figure-of-merit (see lablmm for descriptions)
+Switch for figure-of-merit (see `FiguresOfMerit` for descriptions)
 negative => maximise, positive => minimise
 """
 
-lablmm: list[str] = None
-"""lablmm(ipnfoms) : labels describing figures of merit:<UL>
-* ( 1) major radius
-* ( 2) not used
-* ( 3) neutron wall load
-* ( 4) P_tf + P_pf
-* ( 5) fusion gain Q
-* ( 6) cost of electricity
-* ( 7) capital cost (direct cost if ireactor=0,
-constructed cost otherwise)
-* ( 8) aspect ratio
-* ( 9) divertor heat load
-* (10) toroidal field
-* (11) total injected power
-* (12) hydrogen plant capital cost OBSOLETE
-* (13) hydrogen production rate OBSOLETE
-* (14) pulse length
-* (15) plant availability factor (N.B. requires
-i_plant_availability=1 to be set)
-* (16) linear combination of major radius (minimised) and pulse length (maximised)
-note: FoM should be minimised only!
-* (17) net electrical output
-* (18) Null Figure of Merit
-* (19) linear combination of big Q and pulse length (maximised)
-note: FoM should be minimised only!
-"""
-
 n_constraints: int = None
 """Total number of constraints (neqns + nineqns)"""
 
@@ -454,7 +481,6 @@ def init_numerics():
     global \
         ioptimz, \
         minmax, \
-        lablmm, \
         n_constraints, \
         ncalls, \
         neqns, \
@@ -492,27 +518,6 @@ def init_numerics():
     """Initialise module variables"""
     ioptimz = 1
     minmax = 7
-    lablmm = [
-        "major radius          ",
-        "not used              ",
-        "neutron wall load     ",
-        "P_tf + P_pf           ",
-        "fusion gain           ",
-        "cost of electricity   ",
-        "capital cost          ",
-        "aspect ratio          ",
-        "divertor heat load    ",
-        "toroidal field        ",
-        "total injected power  ",
-        "H plant capital cost  ",
-        "H production rate     ",
-        "pulse length          ",
-        "plant availability    ",
-        "min R0, max tau_burn  ",
-        "net electrical output ",
-        "Null figure of merit  ",
-        "max Q, max t_plant_pulse_burn     ",
-    ]
 
     ncalls = 0
     neqns = -1