diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md index 5e4f16160..ff8cce68d 100644 --- a/documentation/source/development/add-vars.md +++ b/documentation/source/development/add-vars.md @@ -80,23 +80,22 @@ 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. -4. Add the new figure of merit description to the `OBJECTIVES_NAMES` dictionary in `objectives.py` - An example can be found below: ```python objective_function(): ... - match figure_of_merit: + try: + figure_of_merit = FiguresOfMerit(abs(minmax)) ... - case 1: + if figure_of_merit == FiguresOfMerit.MAJOR_RADIUS: objective_metric = 0.2 * physics_variables.rmajor ``` diff --git a/documentation/source/io/input-guide.md b/documentation/source/io/input-guide.md index 4176d5ba1..d929d44b6 100644 --- a/documentation/source/io/input-guide.md +++ b/documentation/source/io/input-guide.md @@ -118,9 +118,12 @@ 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) ``` +!!! Info "Figures of Merit" + The full list of valid `minmax` values is documented in the API reference for [`process.data_structure.numerics.FiguresOfMerit`](../../source/reference/process/data_structure/numerics/#process.data_structure.numerics.FiguresOfMerit). + In this case the user is choosing option `1`, which is major radius. For `minmax` * a **positive** value means **minimise** the figure of merit @@ -132,10 +135,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: diff --git a/process/core/init.py b/process/core/init.py index 1ecc36d60..e042e0cfd 100644 --- a/process/core/init.py +++ b/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}", diff --git a/process/core/io/plot/solutions.py b/process/core/io/plot/solutions.py index 1e3848e13..f4fce3a19 100644 --- a/process/core/io/plot/solutions.py +++ b/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: diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index a0180e00c..0c81ea02a 100644 --- a/process/core/io/plot/summary.py +++ b/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" diff --git a/process/core/scan.py b/process/core/scan.py index 87cf5cff6..80f22848f 100644 --- a/process/core/scan.py +++ b/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 diff --git a/process/core/solver/objectives.py b/process/core/solver/objectives.py index 716607e44..53820bf9b 100644 --- a/process/core/solver/objectives.py +++ b/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.P_TF_PLUS_P_PF: + objective_metric = (tfcoil_variables.tfcmw + 1e-3 * data.pf_power.srcktpm) / 10.0 + elif figure_of_merit == FiguresOfMerit.FUSION_GAIN_Q: + 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 diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 79f1feb67..f1e1c2bdb 100644 --- a/process/data_structure/numerics.py +++ b/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: