diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md index e925ea093d..703a648e52 100644 --- a/documentation/source/development/add-vars.md +++ b/documentation/source/development/add-vars.md @@ -1,6 +1,6 @@ # Guide for adding variables and constraints -A guide for how to add new inputs, constraints, figures of merit, and scan variables. +A guide for how to add new inputs, constraints, figures of merit, and scan variables. **At all times ensure new variables and functions adhere to the [`PROCESS` style guide](../development/standards.md).** @@ -23,7 +23,7 @@ It is advised that a variable is either used to define a particular PROCESS run ## Add a new variable You may need to add a variable to PROCESS when changing or creating models. In most cases, you will want to add your variable to an existing data structure. Creating an entirely new data structure is beyond the scope of this guide, so please seek support from the PROCESS maintainers. -For example, if you are adding a new variable that relates to the blanket model, you would add the variable to `process/data_structure/blanket_variables.py` as part of the `BlanketData` dataclass. +For example, if you are adding a new variable that relates to the blanket model, you would add the variable to `process/data_structure/blanket_variables.py` as part of the `BlanketData` dataclass. ```python @dataclass(slots=True) @@ -116,7 +116,7 @@ class FiguresOfMerit(IntEnum): ... BLANKET_FIGURE_OF_MERIT = (20, "my FOM description") ``` -Here `20` will be the identifier of the figure of merit, and **must** be unique. +Here `20` will be the identifier of the figure of merit, and **must** be unique. Finally, add the equation to `process/core/solver/objectives.py`: ```python @@ -136,26 +136,22 @@ Remember, setting `minmax = -20` would minimise instead of maximise our new vari ## Add a scan variable -After following the instructions to add an input variable, you can then make a scan variable. +After following the instruction to add an input variable, you can make the variable a scan variable by following these steps: -First, add the variable to the `ScanVariables` enum in `process/core/scan.py`. -```python -class ScanVariables(Enum): - ... - blanket_scan_variable = ScanVariable( - "my_new_blanket_variable", "A blanket variable", 82 - ) -``` -Here, `82` is the identifier of the scan variable and must be unique. +1. Update the `ScanVariables` enum in the `scan.py` file by adding a new case statement connecting the variable to the scan integer switch, the variable name and a short description. -Next, increment the parameter `IPNSCNV` in `process/data_structure/scan_variables.py` and be sure to add a description of the scan variable in the docstring of the `nsweep` variable. +2. Add a comment in the corresponding variable file in the data_structure directory, eg, `data_structure/[XX]_variables.py`, to add the variable description indicating the scan switch number. -Finally, in `process/core/scan.py`, add the scan variable to the `Scan.scan_select()` method. +`SCAN_VARIABLES` case example: + +Add the variable to the `ScanVariables` enum in `process/core/scan.py`. ```python -match nwp: - ... - case 82: - self.data.tfcoil.my_new_blanket_variable = swp[iscn - 1] + class ScanVariables(ScanVariable, Enum): + aspect = (1, SVE.P), + pflux_div_heat_load_max_mw = (2, SVE.D) + ... + b_crit_upper_nbti = (54, SVE.T) + dr_shld_inboard = (55, SVE.B), ``` Please see the [scan documentation](../io/input-guide.md#input-guide) for how to set up a scan `IN.DAT` diff --git a/process/core/caller.py b/process/core/caller.py index 149291bc2a..0cf0c486ce 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -7,13 +7,14 @@ from tabulate import tabulate from process.core import constants -from process.core.final import finalise +from process.core import process_output as po from process.core.io.mfile import MFile from process.core.process_output import OutputFileManager, ovarre from process.core.solver import constraints from process.core.solver.iteration_variables import set_scaled_iteration_variable from process.core.solver.objectives import objective_function from process.data_structure.blanket_variables import BlktModelTypes +from process.data_structure.numerics import PROCESSRunMode from process.models.tfcoil.base import TFConductorModel from process.models.tfcoil.superconducting import SuperconductingTFTurnType @@ -394,6 +395,103 @@ def _call_models_once(self, xc: np.ndarray): # FISPACT and LOCA model (not used)- removed +def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None): + """Routine to print out the final point in the scan. + + Writes to OUT.DAT and MFILE.DAT. + + Parameters + ---------- + models : process.main.Models + physics and engineering model objects + data: DataStructure + data structure object to provide data to evaluate the constraints + ifail : int + error flag + non_idempotent_msg : None | str, optional + warning about non-idempotent variables, defaults to None + """ + if ifail == 1: + po.oheadr(constants.NOUT, "Final Feasible Point") + else: + po.oheadr(constants.NOUT, "Final UNFEASIBLE Point") + + # Output relevant to no optimisation + if data.numerics.ioptimz == PROCESSRunMode.EVALUATION: + output_evaluation(data) + + # Print non-idempotence warning to OUT.DAT only + if non_idempotent_msg: + po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES") + po.ocmmnt(constants.NOUT, non_idempotent_msg) + + # Write output to OUT.DAT and MFILE.DAT + models.write(data, constants.NOUT) + + +def output_evaluation(data): + """Write output for an evaluation run of PROCESS + + Parameters + ---------- + data: DataStructure + data structure object to provide data to evaluate the constraints + """ + po.oheadr(constants.NOUT, "Numerics") + po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.") + po.oblnkl(constants.NOUT) + + # Evaluate objective function + norm_objf = objective_function(data.numerics.minmax, data) + po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf) + + # Print the residuals of the constraint equations + + residual_error, value, residual, symbols, units = constraints.constraint_eqns( + data.numerics.neqns + data.numerics.nineqns, -1, data + ) + + labels = [ + data.numerics.lablcc[j - 1] + for j in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] + ] + + def _fmt(a, units): + return [f"{c} {u}" for c, u in zip(a, units, strict=False)] + + po.write( + constants.NOUT, + tabulate( + { + "Constraint Name": labels, + "Constraint Type": symbols, + "Physical constraint": _fmt(value, units), + "Constraint residual": _fmt(residual, units), + "Normalised residual": residual_error, + }, + headers="keys", + ), + ) + + for i in range(data.numerics.neqns): + constraint_id = data.numerics.icc[i] + po.ovarre( + constants.MFILE, + f"{labels[i]} normalised residue", + f"(eq_con{constraint_id:03d})", + residual_error[i], + ) + + for i in range(data.numerics.nineqns): + constraint_id = data.numerics.icc[data.numerics.neqns + i] + po.ovarre( + constants.MFILE, + f"{labels[data.numerics.neqns + i]}", + f"(ineq_con{constraint_id:03d})", + residual_error[data.numerics.neqns + i], + ) + + def write_output_files( models: Models, data: DataStructure, ifail: int, *, runtime: float | None = None ): diff --git a/process/core/final.py b/process/core/final.py index d5da7593bc..e69de29bb2 100644 --- a/process/core/final.py +++ b/process/core/final.py @@ -1,105 +0,0 @@ -"""Final output at the end of a scan.""" - -from tabulate import tabulate - -from process.core import constants -from process.core import output as op -from process.core import process_output as po -from process.core.solver import constraints -from process.core.solver.objectives import objective_function -from process.data_structure.numerics import PROCESSRunMode - - -def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None): - """Routine to print out the final point in the scan. - - Writes to OUT.DAT and MFILE.DAT. - - Parameters - ---------- - models : process.main.Models - physics and engineering model objects - data: DataStructure - data structure object to provide data to evaluate the constraints - ifail : int - error flag - non_idempotent_msg : None | str, optional - warning about non-idempotent variables, defaults to None - """ - if ifail == 1: - po.oheadr(constants.NOUT, "Final Feasible Point") - else: - po.oheadr(constants.NOUT, "Final UNFEASIBLE Point") - - # Output relevant to no optimisation - if data.numerics.ioptimz == PROCESSRunMode.EVALUATION: - output_evaluation(data) - - # Print non-idempotence warning to OUT.DAT only - if non_idempotent_msg: - po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES") - po.ocmmnt(constants.NOUT, non_idempotent_msg) - - # Write output to OUT.DAT and MFILE.DAT - op.write(models, data, constants.NOUT) - - -def output_evaluation(data): - """Write output for an evaluation run of PROCESS - - Parameters - ---------- - data: DataStructure - data structure object to provide data to evaluate the constraints - """ - po.oheadr(constants.NOUT, "Numerics") - po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.") - po.oblnkl(constants.NOUT) - - # Evaluate objective function - norm_objf = objective_function(data.numerics.minmax, data) - po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf) - - # Print the residuals of the constraint equations - - residual_error, value, residual, symbols, units = constraints.constraint_eqns( - data.numerics.neqns + data.numerics.nineqns, -1, data - ) - - labels = [ - data.numerics.lablcc[j] - for j in [ - i - 1 - for i in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] - ] - ] - physical_constraint = [f"{c} {u}" for c, u in zip(value, units, strict=False)] - physical_residual = [f"{c} {u}" for c, u in zip(residual, units, strict=False)] - - table_data = { - "Constraint Name": labels, - "Constraint Type": symbols, - "Physical constraint": physical_constraint, - "Constraint residual": physical_residual, - "Normalised residual": residual_error, - } - - po.write(constants.NOUT, tabulate(table_data, headers="keys")) - - for i in range(data.numerics.neqns): - constraint_id = data.numerics.icc[i] - po.ovarre( - constants.MFILE, - f"{labels[i]} normalised residue", - f"(eq_con{constraint_id:03d})", - residual_error[i], - ) - - for i in range(data.numerics.nineqns): - constraint_id = data.numerics.icc[data.numerics.neqns + i] - po.ovarre( - constants.MFILE, - f"{labels[data.numerics.neqns + i]}", - f"(ineq_con{constraint_id:03d})", - residual_error[data.numerics.neqns + i], - ) diff --git a/process/core/input.py b/process/core/input.py index 1d7c132b6a..57a0f2757b 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1138,24 +1138,15 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]: "scan", int, choices=range(IPNSCNS + 1), + array=True, ), "nsweep": InputVariable( "scan", int, choices=range(1, IPNSCNV + 1), - ), - "isweep_2": InputVariable( - "scan", - int, - choices=range(IPNSCNS + 1), - ), - "nsweep_2": InputVariable( - "scan", - int, - choices=range(1, IPNSCNV + 1), + array=True, ), "sweep": InputVariable("scan", float, array=True), - "sweep_2": InputVariable("scan", float, array=True), "impvardiv": InputVariable( "reinke", int, diff --git a/process/core/io/plot/cli.py b/process/core/io/plot/cli.py index 804b969a7d..9a8dd6bff6 100644 --- a/process/core/io/plot/cli.py +++ b/process/core/io/plot/cli.py @@ -241,10 +241,10 @@ def plot_scans_cli( list(map(float, filter(None, x_axis_max))), list(map(float, filter(None, x_axis_range))), y_axis_percent, - y_axis_percent2, list(map(float, filter(None, y_axis_max))), - list(map(float, filter(None, y_axis2_max))), list(map(float, filter(None, y_axis_range))), + y_axis_percent2, + list(map(float, filter(None, y_axis2_max))), list(map(float, filter(None, y_axis_range2))), label_name, twod_contour, diff --git a/process/core/io/plot/scans.py b/process/core/io/plot/scans.py index 80ac996f60..e723e8bd2f 100644 --- a/process/core/io/plot/scans.py +++ b/process/core/io/plot/scans.py @@ -22,17 +22,133 @@ - If the file is a folder, the contained MFILE is used as an input. """ +from __future__ import annotations + import math -import sys from collections.abc import Iterable, Sequence +from dataclasses import dataclass +from enum import Enum, auto from pathlib import Path +from typing import TYPE_CHECKING import matplotlib.pyplot as plt -import matplotlib.ticker as mtick import numpy as np +from matplotlib.ticker import MultipleLocator, PercentFormatter from process.core.io.mfile import MFile +from process.core.io.mfile.cli import mfile from process.core.io.variable_metadata import var_dicts as meta +from process.core.scan import ScanVariables + +if TYPE_CHECKING: + from matproplib.axes import Axes + + +@dataclass +class AxisData: + name: str + percent: bool + max_: Sequence[float] + range_: Sequence[float] + tick_size: float + font_size: float + legend_size: float = 12 + + +class AxisChoice(Enum): + X = auto() + Y = auto() + + def axis(self, ax): + return getattr(ax, f"{self.name.lower()}axis") + + def set_lim(self, ax, lower, upper): + getattr(ax, f"set_{self.name.lower()}lim")(lower, upper) + + +def get_list_padded(inp, names): + target_len = len(names) + inp_array = np.array(inp, dtype=float) + if (i_len := len(inp_array)) < target_len: + if i_len == 0: + return [None] * target_len + + return np.concatenate(( + inp_array, + np.full(target_len - i_len, inp_array[-1], dtype=float), + )) + return inp_array[:target_len] + + +def value_checks( + scan_var: ScanVariables, + scan_2_var: ScanVariables | None, + m_file: MFile, + input_files: Sequence[Path], +): + ve_string = ( + "`{}` does not exist in PROCESS dicts\n" + " The scan variable is probably an upper/lower boundary\n" + " Please modify 'nsweep_dict' dict with the constrained var" + ) + # Check if the scan variable is present + if scan_var.name not in m_file.data: + raise ValueError(ve_string.format(scan_var.name)) + + # Check if the second scan variable is present + if scan_2_var is not None: + if scan_2_var.name not in m_file.data: + raise ValueError(ve_string.format(scan_2_var.name)) + + if len(input_files) > 1: + raise ValueError("Only one input file can be used for 2D scans") + + +def array_check(output_name: str, m_file: MFile) -> bool: + # Check if the output variable exists in the MFILE + if output_name not in m_file.data: + print( + f"Warning : `{output_name}` does not exist in PROCESS dicts\n" + f"Warning : `{output_name}` will not be output" + ) + return False + return True + + +def create_o_array( + n_scan: int, m_file: MFile, output_name: str, conv_i: list[int] +) -> np.ndarray: + return np.array([m_file.get(output_name, scan=conv_i[ii]) for ii in range(n_scan)]) + + +def get_label(name: str) -> str: + return meta[name].latex if name in meta else f"{name}" + + +def axis_manipulation(ax: Axes, axis: AxisData, index: int, contour: np.ndarray): + + an = AxisChoice[axis.name.upper()] + + if len(axis.range_) > 0: + divisions = (axis.range_[1] - axis.range_[0]) / 10 + if axis.percent: + if axis.max_[index] is None: + axis.max_[index] = max(np.abs(contour)) + ticks = PercentFormatter(axis.max_[index]) + if len(axis.range_) > 0: + scale = axis.max_[index] / 100 + divisions = 5 * math.ceil(divisions / 5) * scale + range_ = (axis.range_[0] * scale, axis.range_[1] * scale) + an.axis(ax).set_major_formatter(ticks) + + if len(axis.range_) > 0: + if axis.percent is False: + range_ = axis.range_ + an.set_lim(ax, range_[0], range_[1]) + an.axis(ax).set_major_locator(MultipleLocator(divisions)) + + ax.figure.tight_layout() + ax.tick_params(axis=an.name.lower(), labelsize=axis.tick_size) def plot_scan( @@ -48,10 +164,10 @@ def plot_scan( x_axis_max: Sequence[float] = (), x_axis_range: Sequence[float] = (), y_axis_percent: bool = False, - y_axis_percent2: bool = False, y_axis_max: Sequence[float] = (), - y_axis2_max: Sequence[float] = (), y_axis_range: Sequence[float] = (), + y_axis_percent2: bool = False, + y_axis2_max: Sequence[float] = (), y_axis_range2: Sequence[float] = (), label_name: Sequence[str] = (), twod_contour: bool = False, @@ -60,908 +176,445 @@ def plot_scan( """Main plot scans script.""" outputdir = outputdir or Path.cwd() input_files = mfiles if isinstance(mfiles, Iterable) else [mfiles] - x_max_input = x_axis_max - - y_max_input = y_axis_max - y_max2_input = y_axis2_max # If the input file is a directory, add MFILE.DAT for ii, if_ in enumerate(input_files): if if_.is_dir(): input_files[ii] = if_ / "MFILE.DAT" - # nsweep variable dict - # ------------------- - # TODO WOULD BE GREAT TO HAVE IT AUTOMATICALLY GENERATED ON THE PROCESS CMAKE! - # THE SAME WAY THE DICTS ARE - # This needs to be kept in sync automatically; this will break frequently - # otherwise - # Rem : Some variables are not in the MFILE, making the definition rather tricky... - nsweep_dict = { - 1: "aspect", - 2: "pflux_div_heat_load_max_mw", - 3: "p_plant_electric_net_mw", - 4: "hfact", - 5: "j_tf_coil_full_area", - 6: "pflux_fw_neutron_max_mw", - 7: "beamfus0", - 8: "Obsolete", # OBSOLETE - 9: "temp_plasma_electron_vol_avg_kev", - 10: "boundu(15)", - 11: "beta_norm_max", - 12: "f_c_plasma_bootstrap_max", - 13: "boundu(10)", - 14: "f_j_tf_wp_critical_max", - 16: "rmajor", - # b_tf_inboard_max the maximum T field upper limit is the scan variable - 17: "b_tf_inboard_peak_symmetric", - 18: "eta_cd_norm_hcd_primary_max", - 19: "boundl(16)", - 20: "cnstv.t_burn_min", - 21: "", - 22: "f_t_plant_available", - 23: "boundu(72)", - 24: "p_fusion_total_max_mw", - 25: "kappa", - 26: "triang", - 27: "tbrmin", - 28: "b_plasma_toroidal_on_axis", - 29: "radius_plasma_core_norm", - 30: "", # OBSOLETE - 31: "f_t_alpha_energy_confinement_min", - 32: "epsvmc", - 33: "ttarget", - 34: "qtargettotal", - 35: "lambda_q_omp", - 36: "lambda_target", - 37: "lcon_factor", - 38: "boundu(129)", - 39: "boundu(131)", - 40: "boundu(135)", - 41: "dr_blkt_outboard", - 42: "f_nd_impurity_electrons(9)", - 43: "Obsolete", # OBSOLETE - 44: "alstrtf", - 45: "temp_tf_superconductor_margin_min", - 46: "boundu(152)", - 47: "impurity_enrichment(9)", - 48: "n_tf_wp_pancakes", - 49: "n_tf_wp_layers", - 50: "f_nd_impurity_electrons(13)", - 51: "f_p_div_lower", - 52: "rad_fraction_sol", - 53: "Obsolete", # OBSOLETE - 54: "b_crit_upper_nbti", - 55: "dr_shld_inboard", - 56: "p_cryo_plant_electric_max_mw", - # Genuinely b_plasma_toroidal_on_axis lower bound - 57: "b_plasma_toroidal_on_axis", - 58: "dr_fw_plasma_gap_inboard", - 59: "dr_fw_plasma_gap_outboard", - 60: "sig_tf_wp_max", - 61: "copperaoh_m2_max", - 62: "j_cs_flat_top_end", - 63: "dr_cs", - 64: "f_z_cs_tf_internal", - 65: "n_cycle_min", - 66: "f_a_cs_turn_steel", - 67: "t_crack_vertical", - 68: "inlet_temp_liq", - 69: "outlet_temp_liq", - 70: "blpressure_liq", - 71: "n_liq_recirc", - 72: "bz_channel_conduct_liq", - 73: "pnuc_fw_ratio_dcll", - 74: "f_nuc_pow_bz_struct", - 75: "dx_fw_module", - 76: "eta_turbine", - 77: "startupratio", - 78: "fkind", - 79: "eta_ecrh_injector_wall_plug", - 80: "fcoolcp", - 81: "n_tf_coil_turns", - } - # ------------------- - # Getting the scanned variable name m_file = MFile(filename=input_files[-1]) - nsweep_ref = int(m_file.data["nsweep"].get_scan(-1)) - scan_var_name = nsweep_dict[nsweep_ref] - # Get the eventual second scan variable - nsweep_2_ref = 0 - is_2D_scan = False - scan_2_var_name = "" - if "nsweep_2" in m_file.data: - is_2D_scan = True - nsweep_2_ref = int(m_file.data["nsweep_2"].get_scan(-1)) - scan_2_var_name = nsweep_dict[nsweep_2_ref] - - # Checks - # ------ - # Check if the nsweep dict has been updated - if nsweep_ref > len(nsweep_dict) + 1: - print( - f"ERROR : nsweep = {nsweep_ref} not supported by the utility\n" - "ERROR : Please update the 'nsweep_dict' dict" - ) - sys.exit() - - # Check if the scan variable is present in the - if scan_var_name not in m_file.data: - print( - f"ERROR : `{scan_var_name}` does not exist in PROCESS dicts\n" - "ERROR : The scan variable is probably an upper/lower boundary\n" - "ERROR : Please modify 'nsweep_dict' dict with the constrained var" - ) - sys.exit() + nsweep_ref = int(m_file.get("nsweep", scan=-1)) + scan_var = ScanVariables(nsweep_ref) - # Check if the second scan variable is present in the - if is_2D_scan and (scan_2_var_name not in m_file.data): - print( - f"ERROR : `{scan_2_var_name}` does not exist in PROCESS dicts\n" - "ERROR : The scan variable is probably an upper/lower boundary\n" - "ERROR : Please modify 'nsweep_dict' dict with the constrained var" - ) - sys.exit() - - # Only one input must be used for a 2D scan - if is_2D_scan and len(input_files) > 1: - print("ERROR : Only one input file can be used for 2D scans\nERROR : Exiting") - sys.exit() - # ------ - - # Plot settings - # ------------- - # Plot cosmetic settings - def _format_lists(inp, output_names): - x_max = [] - if len(inp) > 0: - for i in range(len(output_names)): - j = 0 - try: - x_max += [float(inp[i])] - j += 1 - except IndexError: - x_max += [float(inp[j])] - else: - x_max = [None] * len(output_names) + # Get the eventual second scan variable + scan_2_var = ( + ScanVariables(int(m_file.get("nsweep_2", scan=-1))) + if "nsweep_2" in m_file.data + else None + ) - return x_max + value_checks(scan_var, scan_2_var, m_file, input_files) - legend_size = 12 - x_max = ( - _format_lists(x_max_input, output_names) - if len(x_max_input) != len(output_names) - else np.float64(x_max_input) + x_max = get_list_padded(x_axis_max, output_names) + x_axis = AxisData( + "x", x_axis_percent, x_max, x_axis_range, axis_tick_size, axis_font_size ) - y_max = ( - _format_lists(y_max_input, output_names) - if len(y_max_input) != len(output_names) - else np.float64(y_max_input) + + y_max = get_list_padded(y_axis_max, output_names) + y_axis = AxisData( + "y", y_axis_percent, y_max, y_axis_range, axis_tick_size, axis_font_size ) - if len(output_names2) > 0: - y_max2 = ( - _format_lists(y_max2_input, output_names) - if len(y_max2_input) != len(output_names) - else np.float64(y_max2_input) + if scan_2_var is None: + y_axis2 = AxisData( + "y", + y_axis_percent2, + ( + get_list_padded(y_axis2_max, output_names) + if len(output_names2) > 0 + else y_axis2_max + ), + y_axis_range2, + axis_tick_size, + axis_font_size, + ) + + scan_var_array, output_arrays, output_arrays2 = oned_scan( + input_files, + nsweep_ref, + scan_var, + output_names, + output_names2, + term_output=term_output, + ) + plot_1d_scan( + input_files, + m_file, + output_names, + output_names2, + scan_var, + outputdir, + save_format, + label_name, + scan_var_array, + output_arrays, + output_arrays2, + x_axis, + y_axis, + y_axis2, + stack_plots=stack_plots, ) else: - y_max2 = y_max2_input - # ------------- - - # Case of a set of 1D scans - if not is_2D_scan: - # Loop over the MFILEs - output_arrays = {} - output_arrays2 = {} - scan_var_array = {} - for input_file in input_files: - # Opening the MFILE.DAT - m_file = MFile(filename=input_file) - - # Check if the the scan variable is the same for all inputs - # --- - # Same scan var - nsweep = int(m_file.data["nsweep"].get_scan(-1)) - if nsweep != nsweep_ref: - print( - "ERROR : You must use inputs files with the same scan variables\n" - "ERROR : Exiting" - ) - sys.exit() - - # No D scans - if "nsweep_2" in m_file.data: - print("ERROR : You cannot mix 1D with 2D scans\nERROR : Exiting") - sys.exit() - # --- - - # Only selecting the scans that has converged - # --- - # Number of scan points - n_scan = int(m_file.data["isweep"].get_scan(-1)) - - # Converged indexes - conv_i = [] - for ii in range(n_scan): - ifail = m_file.data["ifail"].get_scan(ii + 1) - if ifail == 1: - conv_i.append(ii + 1) - else: - failed_value = m_file.data[scan_var_name].get_scan(ii + 1) - print( - "Warning : Non-convergent scan point : " - f"{scan_var_name} = {failed_value}\n" - "Warning : This point will not be shown." - ) + twod_scan( + input_files, + scan_var, + scan_2_var, + output_names, + outputdir, + save_format, + x_axis, + y_axis, + twod_contour=twod_contour, + ) - # Updating the number of scans - n_scan = len(conv_i) - # --- - # Scanned variable - scan_var_array[input_file] = np.zeros(n_scan) - for ii in range(n_scan): - scan_var_array[input_file][ii] = m_file.data[scan_var_name].get_scan( - conv_i[ii] +def oned_scan( + input_files: Sequence[Path], + nsweep_ref: int, + scan_var: ScanVariables, + output_names: Sequence[str], + output_names2: Sequence[str], + *, + term_output: bool, +) -> tuple[np.ndarray, ...]: + # input file, output_name, scan + output_arrays = {} + # input file, output_name2, scan + output_arrays2 = {} + # input_file, scan + scan_var_array = {} + for input_file in input_files: + # Opening the MFILE.DAT + m_file = MFile(filename=input_file) + n_scan = int(m_file.get("isweep", scan=-1)) + + # Check if the the scan variable is the same for all inputs + # Same scan var + nsweep = int(m_file.get("nsweep", scan=-1)) + if nsweep != nsweep_ref: + raise ValueError("You must use inputs files with the same scan variables\n") + + if "nsweep_2" in m_file.data: + raise ValueError("You cannot mix 1D with 2D scans\nERROR : Exiting") + + # Only selecting the scans that has converged + # Converged indexes + conv_i = [] + for ii in range(n_scan): + ifail = m_file.get("ifail", scan=ii + 1) + if ifail == 1: + conv_i.append(ii + 1) + else: + failed_value = scan_var.get_val(m_file, scan=ii + 1) + print( + f"Warning : Non-convergent scan point : {scan_var.name} = {failed_value}\n" + "Warning : This point will not be shown." ) - # output list declaration - output_arrays[input_file] = {} - output_arrays2[input_file] = {} - # First variable scan - for output_name in output_names: - ouput_array = np.zeros(n_scan) - - # Check if the output variable exists in the MFILE - if output_name not in m_file.data: - print( - f"Warning : `{output_name}` does not exist in PROCESS dicts\n" - f"Warning : `{output_name}` will not be output" - ) - continue - - for ii in range(n_scan): - ouput_array[ii] = m_file.data[output_name].get_scan(conv_i[ii]) - output_arrays[input_file][output_name] = ouput_array - # Second variable scan - for output_name2 in output_names2: - ouput_array2 = np.zeros(n_scan) - - # Check if the output variable exists in the MFILE - if output_name2 not in m_file.data: - print( - f"Warning : `{output_name2}` does not exist in PROCESS dicts\n" - f"Warning : `{output_name2}` will not be output" - ) - continue - for ii in range(n_scan): - ouput_array2[ii] = m_file.data[output_name2].get_scan(conv_i[ii]) - output_arrays2[input_file][output_name2] = ouput_array2 - # Terminal output - if term_output: + # Updating the number of scans + n_scan = len(conv_i) + scan_var_array[input_file] = np.array([ + scan_var.get_val(mfile, scan=conv_i[ii]) for ii in range(n_scan) + ]) + output_arrays[input_file] = { + output_name: create_o_array(n_scan, m_file, output_name, conv_i) + for output_name in output_names + if array_check(output_name, m_file) + } + output_arrays2[input_file] = { + output_name2: create_o_array(n_scan, m_file, output_name2, conv_i) + for output_name2 in output_names2 + if array_check(output_name2, m_file) + } + # Terminal output + if term_output: + print( + f"\n{input_file} scan output\n\nX-axis:\n" + f"scan var {scan_var.name} : {scan_var_array[input_file]}\n\nY-axis:" + + "\n".join( + f"{output_name} : {output_arrays[input_file][output_name]}" + for output_name in output_names + if output_name in m_file.data + ) + + "\n" + ) + if len(output_names2) > 0: + last_name = output_names2[-1] print( - f"\n{input_file} scan output\n\nX-axis:\n" - f"scan var {scan_var_name} : {scan_var_array[input_file]}\n\nY-axis:" + f"Y2-Axis\n {last_name} : {output_arrays2[input_file][last_name]}\n" ) - for output_name in output_names: - # Check if the output variable exists in the MFILE - if output_name not in m_file.data: - continue - - print(f"{output_name} : {output_arrays[input_file][output_name]}") - print() - if len(output_names2) > 0: - print( - f"Y2-Axis\n {output_name2} : " - f"{output_arrays2[input_file][output_name2]}\n" - ) - # Plot section - # ----------- - for index, output_name in enumerate(output_names): - if not stack_plots: - fig, ax = plt.subplots() - if len(output_names2) > 0: - ax2 = ax.twinx() - # reset counter for label_name - kk = 0 - - # Check if the output variable exists in the MFILE - if output_name not in m_file.data: - continue - - # Loop over inputs - for input_file in input_files: - # Legend label formating - if len(label_name) == 0: - labl = input_file.name - else: - labl = label_name[kk] - kk += 1 - - # Plot the graph - if len(output_names2) > 0 and not stack_plots: - ax.plot( - scan_var_array[input_file], - output_arrays[input_file][output_name], - "--o", - color="blue" if len(input_files) == 1 else None, - label=labl, - ) - if len(y_axis_range) > 0: - y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10 - if y_axis_percent: - if y_max[index] is None: - y_max[index] = max( - np.abs(output_arrays[input_file][output_name]) - ) - yticks = mtick.PercentFormatter(y_max[index]) - if len(y_axis_range) > 0: - y_divisions = ( - 5 * math.ceil(y_divisions / 5) * y_max[index] / 100 - ) - y_range = ( - y_axis_range[0] * y_max[index] / 100, - y_axis_range[1] * y_max[index] / 100, - ) - ax.yaxis.set_major_formatter(yticks) - if len(y_axis_range) > 0: - if y_axis_percent is False: - y_range = y_axis_range - ax.set_ylim(y_range[0], y_range[1]) - ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions)) - if len(x_axis_range) > 0: - x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10 - if x_axis_percent: - if x_max[index] is None: - x_max[index] = max(np.abs(scan_var_array[input_file])) - xticks = mtick.PercentFormatter(x_max[index]) - ax.xaxis.set_major_formatter(xticks) - if len(x_axis_range) > 0: - x_divisions = ( - 5 * math.ceil(x_divisions / 5) * x_max[index] / 100 - ) - x_range = ( - x_axis_range[0] * x_max[index] / 100, - x_axis_range[1] * x_max[index] / 100, - ) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - if len(x_axis_range) > 0: - if x_axis_percent is False: - x_range = x_axis_range - plt.xlim(x_range[0], x_range[1]) - ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions)) - plt.tight_layout() - elif stack_plots: - # check stack plots will work - if len(output_names) <= 1: - raise ValueError( - "stack_plots requires at least two output variables" - ) - # Create subplots only once for the first output - if index == 0: - fig, axs = plt.subplots( - len(output_names), - 1, - figsize=(8.0, (3.5 + (1 * len(output_names)))), - sharex=True, - ) - fig.subplots_adjust(hspace=0.0) - - axs[index].plot( - scan_var_array[input_file], - output_arrays[input_file][output_name], - "--o", - color="blue" if len(output_names2) > 0 else None, - label=labl, - ) - if len(y_axis_range) > 0: - y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10 - if y_axis_percent: - if y_max[index] is None: - y_max[index] = max( - np.abs(output_arrays[input_file][output_name]) - ) - yticks = mtick.PercentFormatter(y_max[index]) - if len(y_axis_range) > 0: - y_divisions = ( - 5 * math.ceil(y_divisions / 5) * y_max[index] / 100 - ) - y_range = ( - y_axis_range[0] * y_max[index] / 100, - y_axis_range[1] * y_max[index] / 100, - ) - axs[output_names.index(output_name)].yaxis.set_major_formatter( - yticks - ) - if len(y_axis_range) > 0: - if y_axis_percent is False: - y_range = y_axis_range - axs[output_names.index(output_name)].set_ylim( - y_range[0], y_range[1] - ) - axs[output_names.index(output_name)].yaxis.set_major_locator( - mtick.MultipleLocator(y_divisions) - ) - if len(x_axis_range) > 0: - x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10 - if x_axis_percent: - if x_max[index] is None: - x_max[index] = max(np.abs(scan_var_array[input_file])) - xticks = mtick.PercentFormatter(x_max[index]) - if len(x_axis_range) > 0: - x_divisions = ( - 5 * math.ceil(x_divisions / 5) * x_max[index] / 100 - ) - x_range = ( - x_axis_range[0] * x_max[index] / 100, - x_axis_range[1] * x_max[index] / 100, - ) - axs[output_names.index(output_name)].xaxis.set_major_formatter( - xticks - ) - if len(x_axis_range) > 0: - if x_axis_percent is False: - x_range = x_axis_range - plt.xlim(x_range[0], x_range[1]) - axs[output_names.index(output_name)].xaxis.set_major_locator( - mtick.MultipleLocator(x_divisions) - ) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - plt.tight_layout() - else: - ax.plot( - scan_var_array[input_file], - output_arrays[input_file][output_name], - "--o", - color="blue" if len(output_names2) > 0 else None, - label=labl, - ) - if len(y_axis_range) > 0: - y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10 - if y_axis_percent: - if y_max[index] is None: - y_max[index] = max( - np.abs(output_arrays[input_file][output_name]) - ) - yticks = mtick.PercentFormatter(y_max[index]) - if len(y_axis_range) > 0: - y_divisions = ( - 5 * math.ceil(y_divisions / 5) * y_max[index] / 100 - ) - y_range = ( - y_axis_range[0] * y_max[index] / 100, - y_axis_range[1] * y_max[index] / 100, - ) - ax.yaxis.set_major_formatter(yticks) - if len(y_axis_range) > 0: - if y_axis_percent is False: - y_range = y_axis_range - ax.set_ylim(y_range[0], y_range[1]) - ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions)) - if len(x_axis_range) > 0: - x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10 - if x_axis_percent: - if x_max[index] is None: - x_max[index] = max(np.abs(scan_var_array[input_file])) - xticks = mtick.PercentFormatter(x_max[index]) - if len(x_axis_range) > 0: - x_divisions = ( - 5 * math.ceil(x_divisions / 5) * x_max[index] / 100 - ) - x_range = ( - x_axis_range[0] * x_max[index] / 100, - x_axis_range[1] * x_max[index] / 100, - ) - ax.xaxis.set_major_formatter(xticks) - if len(x_axis_range) > 0: - if x_axis_percent is False: - x_range = x_axis_range - plt.xlim(x_range[0], x_range[1]) - ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions)) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - plt.tight_layout() - if len(output_names2) > 0: + return scan_var_array, output_arrays, output_arrays2 + + +def plot_1d_scan( + input_files: Sequence[Path], + m_file: MFile, + output_names: Sequence[str], + output_names2: Sequence[str], + scan_var: ScanVariables, + outputdir: Path, + save_format: str, + label_name: Sequence[str], + scan_var_array: np.ndarray, + output_arrays: np.ndarray, + output_arrays2: np.ndarray, + x_axis: AxisData, + y_axis: AxisData, + y_axis2: AxisData, + *, + stack_plots: bool, +): + + if stack_plots: + # check stack plots will work + if len(output_names) <= 1: + raise ValueError("stack_plots requires at least two output variables") + # Create subplots only once for the first output + fig, axs = plt.subplots( + len(output_names), + 1, + figsize=(8.0, (3.5 + (1 * len(output_names)))), + sharex=True, + ) + fig.subplots_adjust(hspace=0.0) + + colour = ( # be careful changing this + ("blue" if len(output_names2) > 0 else None) + if len(output_names2) <= 0 or stack_plots + else ("blue" if len(input_files) == 1 else None) + ) + + for index, output_name in enumerate(output_names): + # reset counter for label_name + kk = 0 + + if output_name not in m_file.data: + continue + + if stack_plots: + ax_ = axs[index] + ax = axs[output_names.index(output_name)] + else: + fig, ax = plt.subplots() + if len(output_names2) > 0: + ax2 = ax.twinx() + ax_ = ax + + for input_file in input_files: + if len(label_name) == 0: + labl = input_file.name + else: + labl = label_name[kk] + kk += 1 + + ax_.plot( + scan_var_array[input_file], + output_arrays[input_file][output_name], + "--o", + color=colour, + label=labl, + ) + + axis_manipulation( + ax, axis=x_axis, index=index, contour=scan_var_array[input_file] + ) + axis_manipulation( + ax, + axis=y_axis, + index=index, + contour=output_arrays[input_file][output_name], + ) + + if len(output_names2) > 0: + for output_name2 in output_names2: + yval = output_arrays2[input_file][output_name2] + colour = "red" if len(input_files) == 1 else None ax2.plot( scan_var_array[input_file], - output_arrays2[input_file][output_name2], + yval, "--o", - color="red" if len(input_files) == 1 else None, + color=colour, label=labl, ) ax2.set_ylabel( - ( - meta[output_name2].latex - if output_name2 in meta - else f"{output_name2}" - ), - fontsize=axis_font_size, - color="red" if len(input_files) == 1 else "black", + get_label(output_name2), + fontsize=y_axis.font_size, + color=colour or "black", ) - if len(y_axis_range2) > 0: - y_divisions2 = (y_axis_range2[1] - y_axis_range2[0]) / 10 - if y_axis_percent2: - if y_max2[index] is None: - y_max2[index] = max( - np.abs(output_arrays2[input_file][output_name2]) - ) - yticks2 = mtick.PercentFormatter(y_max2[index]) - if len(y_axis_range2) > 0: - y_divisions2 = ( - 5 * math.ceil(y_divisions2 / 5) * y_max2[index] / 100 - ) - y_range2 = ( - y_axis_range2[0] * y_max2[index] / 100, - y_axis_range2[1] * y_max2[index] / 100, - ) - ax2.yaxis.set_major_formatter(yticks2) - if len(y_axis_range2) > 0: - if y_axis_percent2 is False: - y_range2 = y_axis_range2 - ax2.set_ylim(y_range2[0], y_range2[1]) - ax2.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions2)) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - plt.tight_layout() - if len(output_names2) > 0: - ax2.yaxis.grid(True) - ax.xaxis.grid(True) - ax.set_ylabel( - ( - meta[output_name].latex - if output_name in meta - else f"{output_name}" - ), - fontsize=axis_font_size, - color="blue" if len(input_files) == 1 else "black", - ) - ax.set_xlabel( - ( - meta[scan_var_name].latex - if scan_var_name in meta - else f"{scan_var_name}" - ), - fontsize=axis_font_size, - ) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - if len(input_files) != 1: - plt.legend(loc="best", fontsize=legend_size) - plt.tight_layout() - elif stack_plots: - axs[output_names.index(output_name)].minorticks_on() - axs[output_names.index(output_name)].grid(True) - axs[output_names.index(output_name)].set_ylabel( - ( - meta[output_name].latex - if output_name in meta - else f"{output_name}" - ), + axis_manipulation(ax2, axis=y_axis2, index=index, contour=yval) + if len(output_names2) > 0: + ax2.yaxis.grid(True) + ax.xaxis.grid(True) + ax.set_ylabel( + get_label(output_name), + fontsize=y_axis.font_size, + color="blue" if len(input_files) == 1 else "black", + ) + ax.set_xlabel( + get_label(scan_var.name), + fontsize=x_axis.font_size, + ) + if len(input_files) != 1: + fig.legend(loc="best", fontsize=x_axis.legend_size) + elif stack_plots: + ax.minorticks_on() + ax.grid(True) + ax.set_ylabel(get_label(output_name)) + ax.set_xlabel(get_label(scan_var.name), fontsize=x_axis.font_size) + + ymin, ymax = ax.get_ylim() + if ymin < 0 and ymax > 0: + mod_min = ymin * 1.1 + mod_max = ymax * 1.1 + elif ymin >= 0: + mod_min = ymin * 0.9 + mod_max = ymax * 1.1 + else: + mod_min = ymin * 1.1 + mod_max = ymax * 0.9 + ax.set_ylim(mod_min, mod_max) + + if len(input_files) > 1: + fig.legend( + loc="lower center", + fontsize=x_axis.legend_size, + bbox_to_anchor=(0.5, -0.5 - (0.1 * len(output_names))), + fancybox=True, + shadow=False, + ncol=len(input_files), + columnspacing=0.8, ) - plt.xlabel( - ( - meta[scan_var_name].latex - if scan_var_name in meta - else f"{scan_var_name}" - ), - fontsize=axis_font_size, - ) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - if len(input_files) > 1: - plt.legend( - loc="lower center", - fontsize=legend_size, - bbox_to_anchor=(0.5, -0.5 - (0.1 * len(output_names))), - fancybox=True, - shadow=False, - ncol=len(input_files), - columnspacing=0.8, - ) - plt.tight_layout() - ymin, ymax = axs[output_names.index(output_name)].get_ylim() - if ymin < 0 and ymax > 0: - axs[output_names.index(output_name)].set_ylim(ymin * 1.1, ymax * 1.1) - elif ymin >= 0: - axs[output_names.index(output_name)].set_ylim(ymin * 0.9, ymax * 1.1) - else: - axs[output_names.index(output_name)].set_ylim(ymin * 1.1, ymax * 0.9) + else: + ax.grid(True) + ax.set_ylabel( + get_label(output_name), + fontsize=x_axis.font_size, + color="red" if len(output_names2) > 0 else "black", + ) + ax.set_xlabel(get_label(scan_var.name), fontsize=x_axis.font_size) + + fig.title( + f"{get_label(output_name)} vs {get_label(scan_var.name)}", + fontsize=x_axis.font_size, + ) + if len(input_files) != 1: + fig.legend(loc="best", fontsize=x_axis.legend_size) + + ax.tick_params(axis=x_axis.name.lower(), labelsize=x_axis.tick_size) + ax.tick_params(axis=y_axis.name.lower(), labelsize=y_axis.tick_size) + fig.tight_layout() + + # Output file naming + # This uses exclusively the last output_name2 defined in an earlier loop ignoring all other output_name2s... + if output_name == "plasma_current_MA": + extra_str = f"plasma_current{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}" + elif stack_plots and output_names[-1] == output_name: + extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else '_vs_'.join(output_names)}" + else: + extra_str = ( + f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}" + ) + + if (not stack_plots) or (stack_plots and output_names[-1] == output_name): + fig.savefig( + f"{outputdir}/scan_{scan_var.name}_vs_{extra_str}.{save_format}", + dpi=300, + ) + plt.show() + + +def twod_scan( + input_files: Sequence[Path], + scan_var: ScanVariables, + scan_2_var: ScanVariables, + output_names: Sequence[str], + outputdir: Path, + save_format: str, + x_axis: AxisData, + y_axis: AxisData, + *, + twod_contour: bool, +): + m_file = MFile(filename=input_files[0]) + + # Number of scan points + n_scan_1 = int(m_file.get("isweep", scan=-1)) + n_scan_2 = int(m_file.get("isweep_2", scan=-1)) + # Selecting the converged runs only + contour_conv_ij = [] # List of non-converged scan point numbers + conv_ij = [] # 2D array of converged scan point numbers (sweep = rows, sweep_2 = columns) + ii_jj = 0 + for ii in range(n_scan_1): + conv_ij.append([]) + for _jj in range(n_scan_2): + ii_jj += 1 # Represents the scan point number in the MFILE + ifail = m_file.get("ifail", scan=ii_jj) + if ifail == 1: + conv_ij[ii].append(ii_jj) # Only appends scan number if scan converged + contour_conv_ij.append(ii_jj) else: - plt.grid(True) - plt.ylabel( - ( - meta[output_name].latex - if output_name in meta - else f"{output_name}" - ), - fontsize=axis_font_size, - color="red" if len(output_names2) > 0 else "black", - ) - plt.xlabel( - ( - meta[scan_var_name].latex - if scan_var_name in meta - else f"{scan_var_name}" - ), - fontsize=axis_font_size, - ) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - plt.title( - f"{meta[output_name].latex if output_name in meta else {output_name}} vs " # noqa: E501 - f"{meta[scan_var_name].latex if scan_var_name in meta else {scan_var_name}}", # noqa: E501 - fontsize=axis_font_size, + failed_value_1 = scan_var.get_val(m_file, scan=ii_jj) + failed_value_2 = scan_2_var.get_val(m_file, scan=ii_jj) + print( + f"Warning : Non-convergent scan point : ({scan_var.name},{scan_2_var.name}) " + f"= ({failed_value_1},{failed_value_2})\n" + "Warning : This point will not be shown." ) - plt.tight_layout() - if len(input_files) != 1: - plt.legend(loc="best", fontsize=legend_size) - plt.tight_layout() - - # Output file naming - if output_name == "plasma_current_MA": - extra_str = f"plasma_current{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}" # noqa: E501 - elif stack_plots and output_names[-1] == output_name: - extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else '_vs_'.join(output_names)}" # noqa: E501 - else: - extra_str = f"{output_name}{f'_vs_{output_name2}' if len(output_names2) > 0 else ''}" # noqa: E501 - if (not stack_plots) or (stack_plots and output_names[-1] == output_name): - plt.savefig( - f"{outputdir}/scan_{scan_var_name}_vs_{extra_str}.{save_format}", - dpi=300, + for index, output_name in enumerate(output_names): + if output_name not in m_file.data: + print( + f"Warning : `{output_name}` does not exist in PROCESS dicts\n" + f"Warning : `{output_name}` will not be output" + ) + continue + + fig, ax = plt.subplots() + x_contour = [scan_2_var.get_val(m_file, scan=i + 1) for i in range(n_scan_2)] + + if twod_contour: + y_contour = [ + scan_var.get_val(m_file, scan=i + 1) + for i in range(1, n_scan_1 * n_scan_2, n_scan_2) + ] + + output_contour_z = np.zeros((n_scan_1, n_scan_2)) + for i in contour_conv_ij: + ind1 = (i - 1) // n_scan_2 + ind2 = (i - 1) % n_scan_2 + output_contour_z[ind1][(ind2 if ind1 % 2 == 0 else (-ind2 - 1))] = ( + m_file.get(output_name, scan=i) ) - plt.show() - plt.clf() - plt.close() - # ------------ + flat_output_z = output_contour_z.flatten() + flat_output_z.sort() - # In case of a 2D scan - else: - # Opening the MFILE.DAT - m_file = MFile(filename=input_files[0]) - - # Number of scan points - n_scan_1 = int(m_file.data["isweep"].get_scan(-1)) - n_scan_2 = int(m_file.data["isweep_2"].get_scan(-1)) - # Selecting the converged runs only - contour_conv_ij = [] # List of non-converged scan point numbers - # 2D array of converged scan point numbers (sweep = rows, sweep_2 = columns) - conv_ij = [] - ii_jj = 0 - for ii in range(n_scan_1): - conv_ij.append([]) - for _jj in range(n_scan_2): - ii_jj += 1 # Represents the scan point number in the MFILE - ifail = m_file.data["ifail"].get_scan(ii_jj) - if ifail == 1: - conv_ij[ii].append( - ii_jj - ) # Only appends scan number if scan converged - contour_conv_ij.append(ii_jj) - else: - failed_value_1 = m_file.data[scan_var_name].get_scan(ii_jj) - failed_value_2 = m_file.data[scan_2_var_name].get_scan(ii_jj) - print( - "Warning : Non-convergent scan point : " - f"({scan_var_name},{scan_2_var_name}) " - f"= ({failed_value_1},{failed_value_2})\n" - "Warning : This point will not be shown." - ) + levels = np.linspace( + next(filter(lambda i: i > 0.0, flat_output_z)), flat_output_z.max(), 50 + ) + contour = ax.contourf(x_contour, y_contour, output_contour_z, levels=levels) - # Looping over requested outputs - for index, output_name in enumerate(output_names): - # Check if the output variable exists in the MFILE - if output_name not in m_file.data: - print( - f"Warning : `{output_name}` does not exist in PROCESS dicts\n" - f"Warning : `{output_name}` will not be output" - ) - continue - - # Declaring the outputs - output_arrays = [] - - if twod_contour: - output_contour_z = np.zeros((n_scan_1, n_scan_2)) - x_contour = [ - m_file.data[scan_2_var_name].get_scan(i + 1) for i in range(n_scan_2) - ] - y_contour = [ - m_file.data[scan_var_name].get_scan(i + 1) - for i in range(1, n_scan_1 * n_scan_2, n_scan_2) - ] - for i in contour_conv_ij: - output_contour_z[((i - 1) // n_scan_2)][ - ( - ((i - 1) % n_scan_2) - if ((i - 1) // n_scan_2) % 2 == 0 - else (-((i - 1) % n_scan_2) - 1) - ) - ] = m_file.data[output_name].get_scan(i) - - flat_output_z = output_contour_z.flatten() - flat_output_z.sort() - fig, ax = plt.subplots() - levels = np.linspace( - next(filter(lambda i: i > 0.0, flat_output_z)), - flat_output_z.max(), - 50, - ) - contour = ax.contourf( - x_contour, - y_contour, - output_contour_z, - levels=levels, - ) + fig.colorbar(contour).set_label( + label=get_label(output_name), size=y_axis.font_size + ) + ax.set_ylabel(get_label(scan_var.name), fontsize=y_axis.font_size) - fig.colorbar(contour).set_label( - label=( - meta[output_name].latex - if output_name in meta - else f"{output_name}" - ), - size=axis_font_size, - ) - plt.ylabel( - ( - meta[scan_var_name].latex - if scan_var_name in meta - else f"{scan_var_name}" - ), - fontsize=axis_font_size, - ) - plt.xlabel( - ( - meta[scan_2_var_name].latex - if scan_2_var_name in meta - else f"{scan_2_var_name}" - ), - fontsize=axis_font_size, - ) - if len(y_axis_range) > 0: - y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10 - if y_axis_percent: - if y_max[index] is None: - y_max[index] = max(np.abs(y_contour)) - yticks = mtick.PercentFormatter(y_max[index]) - if len(y_axis_range) > 0: - y_divisions = 5 * math.ceil(y_divisions / 5) * y_max[index] / 100 - y_range = ( - y_axis_range[0] * y_max[index] / 100, - y_axis_range[1] * y_max[index] / 100, - ) - ax.yaxis.set_major_formatter(yticks) - if len(y_axis_range) > 0: - if y_axis_percent is False: - y_range = y_axis_range - ax.set_ylim(y_range[0], y_range[1]) - ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions)) - if len(x_axis_range) > 0: - x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10 - if x_axis_percent: - if x_max[index] is None: - x_max[index] = max(np.abs(x_contour)) - xticks = mtick.PercentFormatter(x_max[index]) - if len(x_axis_range) > 0: - x_divisions = 5 * math.ceil(x_divisions / 5) * x_max[index] / 100 - x_range = ( - x_axis_range[0] * x_max[index] / 100, - x_axis_range[1] * x_max[index] / 100, - ) - ax.xaxis.set_major_formatter(xticks) - if len(x_axis_range) > 0: - if x_axis_percent is False: - x_range = x_axis_range - plt.xlim(x_range[0], x_range[1]) - ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions)) - plt.rc("xtick", labelsize=axis_tick_size) - plt.rc("ytick", labelsize=axis_tick_size) - plt.tight_layout() - plt.savefig( - outputdir - / f"scan_{output_name}_vs_{scan_var_name}_{scan_2_var_name}.{save_format}" # noqa: E501 - ) - plt.grid(True) - plt.show() - plt.clf() + else: + y_contour = [m_file.get(output_name, scan=i + 1) for i in range(n_scan_2)] - else: - # Converged indexes, for normal 2D line plot - fig, ax = plt.subplots() - for conv_j in ( - conv_ij - ): # conv_j is an array element containing the converged scan numbers - # Scanned variables - scan_1_var_array = np.zeros(len(conv_j)) - scan_2_var_array = np.zeros(len(conv_j)) - output_array = np.zeros(len(conv_j)) - for jj in range(len(conv_j)): - scan_1_var_array[jj] = m_file.data[scan_var_name].get_scan( - conv_j[jj] - ) - scan_2_var_array[jj] = m_file.data[scan_2_var_name].get_scan( - conv_j[jj] - ) - output_array[jj] = m_file.data[output_name].get_scan(conv_j[jj]) - - # Label formating - labl = f"{meta[scan_var_name].latex if scan_var_name in meta else {scan_var_name}} = {scan_1_var_array[0]}" # noqa: E501 - - # Plot the graph - ax.plot(scan_2_var_array, output_array, "--o", label=labl) - - plt.grid(True) - plt.ylabel( + # conv_j is an array element containing the converged scan numbers + for conv_j in conv_ij: + # Scanned variables + scan_1_var_array, scan_2_var_array, output_array = np.array([ ( - meta[output_name].latex - if output_name in meta - else f"{output_name}" - ), - fontsize=axis_font_size, - ) - plt.xlabel( - ( - meta[scan_2_var_name].latex - if scan_2_var_name in meta - else f"{scan_2_var_name}" - ), - fontsize=axis_font_size, - ) - plt.legend(loc="best", fontsize=legend_size) - y_data = [ - m_file.data[output_name].get_scan(i + 1) for i in range(n_scan_2) - ] - if len(y_axis_range) > 0: - y_divisions = (y_axis_range[1] - y_axis_range[0]) / 10 - if y_axis_percent: - if y_max[index] is None: - y_max[index] = max(np.abs(y_data)) - yticks = mtick.PercentFormatter(y_max[index]) - if len(y_axis_range) > 0: - y_divisions = 5 * math.ceil(y_divisions / 5) * y_max[index] / 100 - y_range = ( - y_axis_range[0] * y_max[index] / 100, - y_axis_range[1] * y_max[index] / 100, - ) - ax.yaxis.set_major_formatter(yticks) - if len(y_axis_range) > 0: - if y_axis_percent is False: - y_range = y_axis_range - ax.set_ylim(y_range[0], y_range[1]) - ax.yaxis.set_major_locator(mtick.MultipleLocator(y_divisions)) - x_data = [ - m_file.data[scan_2_var_name].get_scan(i + 1) for i in range(n_scan_2) - ] - if len(x_axis_range) > 0: - x_divisions = (x_axis_range[1] - x_axis_range[0]) / 10 - if x_axis_percent: - if x_max[index] is None: - x_max[index] = max(np.abs(x_data)) - xticks = mtick.PercentFormatter(x_max[index]) - if len(x_axis_range) > 0: - x_divisions = 5 * math.ceil(x_divisions / 5) * x_max[index] / 100 - x_range = ( - x_axis_range[0] * x_max[index] / 100, - x_axis_range[1] * x_max[index] / 100, - ) - ax.xaxis.set_major_formatter(xticks) - if len(x_axis_range) > 0: - if x_axis_percent is False: - x_range = x_axis_range - plt.xlim(x_range[0], x_range[1]) - ax.xaxis.set_major_locator(mtick.MultipleLocator(x_divisions)) - plt.rc("xtick", labelsize=8) - plt.rc("ytick", labelsize=8) - plt.tight_layout() - plt.savefig( - outputdir - / f"scan_{output_name}_vs_{scan_var_name}_{scan_2_var_name}.{save_format}" # noqa: E501 - ) + scan_var.get_val(m_file, scan=conv_j[jj]), + scan_2_var.get_val(m_file, scan=conv_j[jj]), + m_file.get(output_name, scan=conv_j[jj]), + ) + for jj in range(len(conv_j)) + ]).T + + label = f"{get_label(scan_var.name)} = {scan_1_var_array[0]}" + ax.plot(scan_2_var_array, output_array, "--o", label=label) + + ax.set_ylabel(get_label(output_name), fontsize=y_axis.font_size) + fig.legend(loc="best", fontsize=x_axis.legend_size) + + ax.set_xlabel(get_label(scan_2_var.name), fontsize=x_axis.font_size) - # Display plot (used in Jupyter notebooks) - plt.show() - plt.clf() + axis_manipulation(ax, axis=x_axis, index=index, contour=x_contour) + axis_manipulation(ax, axis=y_axis, index=index, contour=y_contour) + ax.grid(True) + fname = f"scan_{output_name}_vs_{scan_var.name}_{scan_2_var.name}.{save_format}" + fig.savefig(outputdir / fname) + plt.show() diff --git a/process/core/output.py b/process/core/output.py index 82c070b44b..8b13789179 100644 --- a/process/core/output.py +++ b/process/core/output.py @@ -1,147 +1 @@ -from process.core.log import logging_model_handler -from process.data_structure.blanket_variables import BlktModelTypes -from process.models.tfcoil.base import TFConductorModel -from process.models.tfcoil.superconducting import ( - SuperconductingTFTurnType, -) - -def write(models, data, _outfile): - """Write the results to the main output file (OUT.DAT). - - Write the program results to a file, in a tidy format. - - Parameters - ---------- - models : process.main.Models - physics and engineering model objects - _outfile : int - Fortran output unit identifier - - """ - # ensure we are capturing warnings that occur in the 'output' stage - # as these are warnings that occur at our solution point. - # So we clear existing warnings - logging_model_handler.start_capturing() - logging_model_handler.clear_logs() - - # Call stellarator output routine instead if relevant - if data.stellarator.istell != 0: - models.stellarator.output() - return - - # Call IFE output routine instead if relevant - if data.ife.ife != 0: - models.ife.output() - return - - # Costs model - # Cost switch values - # No. | model - # ---- | ------ - # 0 | 1990 costs model - # 1 | 2015 Kovari model - # 2 | Custom model - models.costs.output() - - # Availability model - models.availability.output() - - # Physics model - models.physics.output() - - # Detailed physics, currently only done at final point as values are not used - # by any other functions - models.physics_detailed.output() - - # TODO what is this? Not in caller.py? - models.current_drive.output() - - # Pulsed reactor model - models.pulse.output() - - models.divertor.output() - - # Machine Build Model - models.build.output() - - # Cryostat build - models.cryostat.output() - - # Toroidal field coil copper model - if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER: - models.copper_tf_coil.output() - - # Toroidal field coil superconductor model - if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: - tf_turn_type = SuperconductingTFTurnType( - data.superconducting_tfcoil.i_tf_turn_type - ) - if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT: - models.cicc_sctfcoil.output() - elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR: - models.croco_sctfcoil.output() - else: - raise ValueError( - "Unsupported superconducting TF turn type: " - f"{data.superconducting_tfcoil.i_tf_turn_type}" - ) - - # Toroidal field coil aluminium model - if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM: - models.aluminium_tf_coil.output() - - # Tight aspect ratio machine model - if ( - data.physics.itart == 1 - and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING - ): - models.tfcoil.output() - - # Poloidal field coil model - models.pfcoil.output() - - # Structure Model - models.structure.output() - - # Blanket model - # Blanket switch values - # No. | model - # ---- | ------ - # 1 | CCFE HCPB model - # 2 | KIT HCPB model - # 3 | CCFE HCPB model with Tritium Breeding Ratio calculation - # 4 | KIT HCLL model - # 5 | DCLL model - - models.shield.output() - models.vacuum_vessel.output() - - # First wall geometry - models.fw.output() - - if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB: - # CCFE HCPB model - models.ccfe_hcpb.output() - - elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL: - # DCLL model - models.dcll.output() - - # FISPACT and LOCA model (not used)- removed - - # Power model - models.power.output() - - # Vacuum model - models.vacuum.output() - - # Buildings model - models.buildings.output() - - # Water usage in secondary cooling system - models.water_use.output() - - # stop capturing warnings so that Outfile does not end up with - # a lot of non-model logs - logging_model_handler.stop_capturing() diff --git a/process/core/scan.py b/process/core/scan.py index 535eae7f36..2e0aaf8a98 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -2,8 +2,9 @@ import logging import time -from dataclasses import astuple, dataclass +from dataclasses import dataclass, field from enum import Enum +from types import DynamicClassAttribute from typing import TYPE_CHECKING import numpy as np @@ -15,1248 +16,368 @@ from process.core.log import logging_model_handler, show_errors from process.core.solver import constraints from process.core.solver.solver_handler import SolverHandler -from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode -from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData if TYPE_CHECKING: - from process.core.model import DataStructure, Model + from process.core.model import DataStructure + from process.main import Models logger = logging.getLogger(__name__) @dataclass class ScanVariable: - variable_name: str - variable_description: str - variable_num: int - - def __iter__(self): - return iter(astuple(self)[:2]) - - -class ScanVariables(Enum): + number: int + area: SVE = field(repr=False) + + +class SVE(Enum): + P = "physics" + D = "divertor" + C = "constraints" + T = "tfcoil" + TR = "rebco" + CD = "current_drive" + NUM = "numerics" + CST = "costs" + IR = "impurity_radiation" + B = "build" + HT = "heat_transport" + PF = "pf_coil" + CS = "cs_fatigue" + FWBS = "fwbs" + + +class ScanVariables(ScanVariable, Enum): @classmethod - def _missing_(cls, var): - if isinstance(var, int): + def _missing_(cls, value): + if isinstance(value, int): for sv in cls: - if sv.value.variable_num == var: + if sv.number == value: return sv - return super()._missing_(var) - - aspect = ScanVariable("aspect", "Aspect_ratio", 1) - pflux_div_heat_load_max_mw = ScanVariable( - "pflux_div_heat_load_max_mw", "Div_heat_limit_(MW/m2)", 2 - ) - p_plant_electric_net_required_mw = ScanVariable( - "p_plant_electric_net_required_mw", "Net_electric_power_(MW)", 3 - ) - hfact = ScanVariable("hfact", "Confinement_H_factor", 4) - j_tf_coil_full_area = ScanVariable( - "j_tf_coil_full_area", "TF_inboard_leg_J_(MA/m2)", 5 - ) - pflux_fw_neutron_max_mw = ScanVariable( - "pflux_fw_neutron_max_mw", "Allow._wall_load_(MW/m2)", 6 - ) - beamfus0 = ScanVariable("beamfus0", "Beam_bkgrd_multiplier", 7) - temp_plasma_electron_vol_avg_kev = ScanVariable( - "temp_plasma_electron_vol_avg_kev", "Electron_temperature_keV", 9 - ) - boundu15 = ScanVariable("boundu(15)", "Volt-second_upper_bound", 10) - beta_norm_max = ScanVariable("beta_norm_max", "Beta_coefficient", 11) - f_c_plasma_bootstrap_max = ScanVariable( - "f_c_plasma_bootstrap_max", "Bootstrap_fraction", 12 - ) - boundu10 = ScanVariable("boundu(10)", "H_factor_upper_bound", 13) - f_j_tf_wp_critical_max = ScanVariable( - "f_j_tf_wp_critical_max", "TFC_Iop_/_Icrit_margin", 14 - ) - rmajor = ScanVariable("rmajor", "Plasma_major_radius_(m)", 16) - b_tf_inboard_max = ScanVariable("b_tf_inboard_max", "Max_toroidal_field_(T)", 17) - eta_cd_norm_hcd_primary_max = ScanVariable( - "eta_cd_norm_hcd_primary_max", "Maximum_CD_gamma", 18 - ) - boundl16 = ScanVariable("boundl(16)", "CS_thickness_lower_bound", 19) - t_burn_min = ScanVariable("t_burn_min", "Minimum_burn_time_(s)", 20) - f_t_plant_available = ScanVariable( - "f_t_plant_available", "Plant_availability_factor", 22 - ) - p_fusion_total_max_mw = ScanVariable( - "p_fusion_total_max_mw", "Fusion_power_limit_(MW)", 24 - ) - kappa = ScanVariable("kappa", "Plasma_elongation", 25) - triang = ScanVariable("triang", "Plasma_triangularity", 26) - tbrmin = ScanVariable("tbrmin", "Min_tritium_breed._ratio", 27) - b_plasma_toroidal_on_axis = ScanVariable( - "b_plasma_toroidal_on_axis", "Tor._field_on_axis_(T)", 28 - ) - coreradius = ScanVariable("coreradius", "Core_radius", 29) - f_t_alpha_energy_confinement_min = ScanVariable( - "f_t_alpha_energy_confinement_min", "t_alpha_confinement/taueff_lower_limit", 31 - ) - epsvmc = ScanVariable("epsvmc", "VMCON error tolerance", 32) - boundu129 = ScanVariable("boundu(129)", " Neon upper limit", 38) - boundu131 = ScanVariable("boundu(131)", " Argon upper limit", 39) - boundu135 = ScanVariable("boundu(135)", " Xenon upper limit", 40) - dr_blkt_outboard = ScanVariable("dr_blkt_outboard", "Outboard blanket thick.", 41) - f_nd_impurity_electrons9 = ScanVariable( - "f_nd_impurity_electrons(9)", "Argon fraction", 42 - ) - sig_tf_case_max = ScanVariable( - "sig_tf_case_max", "Allowable_stress_in_tf_coil_case_Tresca_(pa)", 44 - ) - temp_tf_superconductor_margin_min = ScanVariable( - "temp_tf_superconductor_margin_min", "Minimum_allowable_temperature_margin", 45 - ) - boundu152 = ScanVariable( - "boundu(152)", "Max allowable f_nd_plasma_separatrix_greenwald", 46 - ) - n_tf_wp_pancakes = ScanVariable("n_tf_wp_pancakes", "TF Coil - n_tf_wp_pancakes", 48) - n_tf_wp_layers = ScanVariable("n_tf_wp_layers", "TF Coil - n_tf_wp_layers", 49) - f_nd_impurity_electrons13 = ScanVariable( - "f_nd_impurity_electrons(13)", "Xenon fraction", 50 - ) - f_p_div_lower = ScanVariable("f_p_div_lower", "lower_divertor_power_fraction", 51) - rad_fraction_sol = ScanVariable("rad_fraction_sol", "SoL radiation fraction", 52) - boundu157 = ScanVariable("boundu(157)", "Max allowable fvssu", 53) - Bc2_0K = ScanVariable("Bc2(0K)", "GL_NbTi Bc2(0K)", 54) - dr_shld_inboard = ScanVariable("dr_shld_inboard", "Inboard neutronic shield", 55) - p_cryo_plant_electric_max_mw = ScanVariable( - "p_cryo_plant_electric_max_mw", "max allowable p_cryo_plant_electric_mw", 56 - ) - boundl2 = ScanVariable("boundl(2)", "b_plasma_toroidal_on_axis minimum", 57) - dr_fw_plasma_gap_inboard = ScanVariable( - "dr_fw_plasma_gap_inboard", "Inboard FW-plasma sep gap", 58 - ) - dr_fw_plasma_gap_outboard = ScanVariable( - "dr_fw_plasma_gap_outboard", "Outboard FW-plasma sep gap", 59 - ) - sig_tf_wp_max = ScanVariable( - "sig_tf_wp_max", "Allowable_stress_in_tf_coil_conduit_Tresca_(pa)", 60 - ) - copperaoh_m2_max = ScanVariable( - "copperaoh_m2_max", "Max CS coil current / copper area", 61 - ) - coheof = ScanVariable("coheof", "CS coil current density at EOF (A/m2)", 62) - dr_cs = ScanVariable("dr_cs", "CS coil thickness (m)", 63) - ohhghf = ScanVariable("ohhghf", "CS height (m)", 64) - n_cycle_min = ScanVariable("n_cycle_min", "CS stress cycles min", 65) - oh_steel_frac = ScanVariable("oh_steel_frac", "CS steel fraction", 66) - t_crack_vertical = ScanVariable( - "t_crack_vertical", "Initial crack vertical size (m)", 67 - ) - inlet_temp_liq = ScanVariable( - "inlet_temp_liq", "Inlet Temperature Liquid Metal Breeder/Coolant (K)", 68 - ) - outlet_temp_liq = ScanVariable( - "outlet_temp_liq", "Outlet Temperature Liquid Metal Breeder/Coolant (K)", 69 - ) - blpressure_liq = ScanVariable( - "blpressure_liq", "Blanket liquid metal breeder/coolant pressure (Pa)", 70 - ) - n_liq_recirc = ScanVariable( - "n_liq_recirc", - "Selected number of liquid metal breeder recirculations per day", - 71, - ) - bz_channel_conduct_liq = ScanVariable( - "bz_channel_conduct_liq", - "Conductance of liquid metal breeder duct walls (A V-1 m-1)", - 72, - ) - pnuc_fw_ratio_dcll = ScanVariable( - "pnuc_fw_ratio_dcll", - "Ratio of FW nuclear power as fraction of total (FW+BB)", - 73, - ) - f_nuc_pow_bz_struct = ScanVariable( - "f_nuc_pow_bz_struct", - "Fraction of BZ power cooled by primary coolant for dual-coolant blanket", - 74, - ) - dx_fw_module = ScanVariable( - "dx_fw_module", "dx_fw_module of first wall cooling channels (m)", 75 - ) - eta_turbine = ScanVariable("eta_turbine", "Thermal conversion eff.", 76) - startupratio = ScanVariable("startupratio", "Gyrotron redundancy", 77) - fkind = ScanVariable("fkind", "Multiplier for Nth of a kind costs", 78) - eta_ecrh_injector_wall_plug = ScanVariable( - "eta_ecrh_injector_wall_plug", "ECH wall plug to injector efficiency", 79 - ) - fcoolcp = ScanVariable("fcoolcp", "Coolant fraction of TF", 80) - n_tf_coil_turns = ScanVariable("n_tf_coil_turns", "Number of turns in TF", 81) - - -class Scan: - """Perform a parameter scan using the Fortran scan module.""" - - def __init__(self, models: Model, solver: str, data: DataStructure): - """Immediately run the run_scan() method. + raise ProcessValueError("Illegal scan variable number", nwp=value) - Parameters - ---------- - models : - Physics and engineering model objects - solver : - Which solver to use, as specified in solver.py - data : - Data structure object - """ - self.models = models - self.solver = solver - self.data = data - self.solver_handler = SolverHandler(models, solver, data) - self.run_scan() - - def run_scan(self): - """Call a solver over a range of values of one of the variables. + def fname(self): + if "__" in self.name: + return self.name.replace("__", "(") + ")" + return self.name - This method calls the optimisation routine VMCON a number of times, by - performing a sweep over a range of values of a particular variable. A - number of output variable values are written to the MFILE.DAT file at - each scan point, for plotting or other post-processing purposes. - """ - if self.data.scan.isweep == 0: - # Solve single problem, rather than an array of problems (scan) - # doopt() can also run just an evaluation - start_time = time.time() - ifail = self.doopt() - write_output_files( - models=self.models, - data=self.data, - ifail=ifail, - runtime=time.time() - start_time, - ) - show_errors(constants.NOUT) - return + def set(self, data: DataStructure, sweep_val: float): + var_area = getattr(data, self.area.value) - if self.data.scan.isweep > IPNSCNS: + if self.number == 22 and var_area.i_plant_availability == 1: raise ProcessValueError( - "Illegal value of isweep", - isweep=self.data.scan.isweep, - IPNSCNS=IPNSCNS, + "Do not scan f_t_plant_available if i_plant_availability=1" ) - if self.data.scan.scan_dim == 2: - self.scan_2d() - else: - self.scan_1d() - - def doopt(self): - """Run the optimiser or solver.""" - ifail = self.solver_handler.run() - self.post_optimise(ifail) - - return ifail - - def post_optimise(self, ifail: int): - """Called after calling the optimising equation solver from Python. - - ifail : input integer : error flag - - Parameters - ---------- - ifail: int : - - """ - self.data.numerics.sqsumsq = ( - sum(r**2 for r in self.data.numerics.rcm[: self.data.numerics.neqns]) ** 0.5 - ) - - process_output.oheadr(constants.NOUT, "Numerics") - if self.solver == "fsolve": - process_output.ocmmnt( - constants.NOUT, "PROCESS has performed an fsolve (evaluation) run." - ) + if "__" in self.name: + name, index = self.name.split("__") + getattr(var_area, name)[int(index)] = sweep_val + if name == "f_nd_impurity_electrons": + var_area.f_nd_impurity_electron_array[int(index - 1)] = sweep_val else: - process_output.ocmmnt( - constants.NOUT, "PROCESS has performed a VMCON (optimisation) run." - ) - if ifail != 1: - process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) - process_output.oheadr( - constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" - ) - process_output.oblnkl(constants.IOTTY) + setattr(var_area, self.name, sweep_val) + name = self.name + + self._data_ = getattr(var_area, name) + + @DynamicClassAttribute + def data(self): + if hasattr(self, "_data_"): + return self._data_ + raise ValueError("Data not available") + + def get_val(self, mfile, scan): + return mfile.get(self.name, scan=scan) + + aspect = (1, SVE.P) + pflux_div_heat_load_max_mw = (2, SVE.D) + p_plant_electric_net_required_mw = (3, SVE.C) + hfact = (4, SVE.P) + j_tf_coil_full_area = (5, SVE.T) + pflux_fw_neutron_max_mw = (6, SVE.C) + beamfus0 = (7, SVE.P) + temp_plasma_electron_vol_avg_kev = (9, SVE.P) + boundu__14 = (10, SVE.NUM) + beta_norm_max = (11, SVE.P) + f_c_plasma_bootstrap_max = (12, SVE.CD) + boundu__10 = (13, SVE.NUM) + f_j_tf_wp_critical_max = ScanVariable(14, SVE.C) # TODO is this needed + rmajor = (16, SVE.P) + b_tf_inboard_max = (17, SVE.C) + eta_cd_norm_hcd_primary_max = (18, SVE.C) + boundl__16 = (19, SVE.NUM) + t_burn_min = (20, SVE.C) + f_t_plant_available = (22, SVE.CST) + p_fusion_total_max_mw = (24, SVE.C) + kappa = (25, SVE.P) + triang = (26, SVE.P) + tbrmin = (27, SVE.C) + b_plasma_toroidal_on_axis = (28, SVE.P) + coreradius = (29, SVE.IR) + f_alpha_energy_confinement_min = (31, SVE.C) + epsvmc = (32, SVE.NUM) + boundu__129 = (38, SVE.NUM) + boundu__131 = (39, SVE.NUM) + boundu__135 = (40, SVE.NUM) + dr_blkt_outboard = (41, SVE.B) + f_nd_impurity_electrons__9 = (42, SVE.IR) + sig_tf_case_max = (44, SVE.T) + temp_tf_superconductor_margin_min = (45, SVE.T) + boundu__152 = (46, SVE.NUM) + n_tf_wp_pancakes = (48, SVE.T) + n_tf_wp_layers = (49, SVE.T) + f_nd_impurity_electrons__13 = (50, SVE.IR) + f_p_div_lower = (51, SVE.P) + rad_fraction_sol = (52, SVE.P) + boundu__157 = (53, SVE.NUM) + b_crit_upper_nbti = (54, SVE.T) + dr_shld_inboard = (55, SVE.B) + p_cryo_plant_electric_max_mw = (56, SVE.HT) + boundl__2 = (57, SVE.NUM) + dr_fw_plasma_gap_inboard = (58, SVE.B) + dr_fw_plasma_gap_outboard = (59, SVE.B) + sig_tf_wp_max = (60, SVE.T) + copperaoh_m2_max = (61, SVE.TR) + coheof = (62, SVE.PF) + dr_cs = (63, SVE.B) + ohhghf = (64, SVE.PF) + n_cycle_min = (65, SVE.CS) + oh_steel_frac = (66, SVE.PF) + t_crack_vertical = (67, SVE.CS) + inlet_temp_liq = (68, SVE.FWBS) + outlet_temp_liq = (69, SVE.FWBS) + blpressure_liq = (70, SVE.FWBS) + n_liq_recirc = (71, SVE.FWBS) + bz_channel_conduct_liq = (72, SVE.FWBS) + pnuc_fw_ratio_dcll = (73, SVE.FWBS) + f_nuc_pow_bz_struct = (74, SVE.FWBS) + dx_fw_module = (75, SVE.FWBS) + eta_turbine = (76, SVE.HT) + startupratio = (77, SVE.CST) + fkind = (78, SVE.CST) + eta_ecrh_injector_wall_plug = (79, SVE.CD) + fcoolcp = (80, SVE.T) + n_tf_coil_turns = (81, SVE.T) - logger.critical("Solver returns with ifail /= 1. %s", ifail) - - # Error code handler for VMCON - if self.solver == "vmcon": - self.verror(ifail) - process_output.oblnkl(constants.NOUT) - process_output.oblnkl(constants.IOTTY) - else: - # Solution found - if self.solver != "fsolve": - process_output.ocmmnt( - constants.NOUT, "and found a feasible set of parameters." - ) - process_output.oheadr( - constants.IOTTY, "PROCESS found a feasible solution" - ) - else: - process_output.ocmmnt( - constants.NOUT, "and found a consistent set of parameters." - ) - process_output.oheadr( - constants.IOTTY, "PROCESS found a consistent solution" - ) - process_output.oblnkl(constants.NOUT) - process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) - if self.data.numerics.sqsumsq >= 1.0e-2: - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, - "WARNING: Constraint residues are HIGH; consider re-running", - ) - process_output.ocmmnt( - constants.NOUT, - " with lower values of EPSVMC to confirm convergence...", - ) - process_output.ocmmnt( - constants.NOUT, - " (should be able to get down to about 1.0E-8 okay)", - ) - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.IOTTY, - "WARNING: Constraint residues are HIGH; consider re-running", - ) - process_output.ocmmnt( - constants.IOTTY, - " with lower values of EPSVMC to confirm convergence...", - ) - process_output.ocmmnt( - constants.IOTTY, - " (should be able to get down to about 1.0E-8 okay)", - ) - process_output.oblnkl(constants.IOTTY) - - logger.warning( - f"High final constraint residues. {self.data.numerics.sqsumsq=}" - ) +@dataclass +class ScanRes: + iscan: int + ifail: int + solver: SolverHandler - process_output.ovarre( - constants.NOUT, - "Number of iteration variables", - "(nvar)", - self.data.numerics.nvar, - ) - process_output.ovarre( - constants.NOUT, - "Number of constraints (total)", - "(neqns+nineqns)", - self.data.numerics.neqns + self.data.numerics.nineqns, - ) - process_output.ovarre( - constants.NOUT, - "Optimisation switch", - "(ioptimz)", - self.data.numerics.ioptimz, - ) - process_output.ocmmnt( - constants.NOUT, - f" {PROCESSRunMode(self.data.numerics.ioptimz).description}", - ) - # Objective function output: none for fsolve - if self.solver != "fsolve": - process_output.ovarre( - constants.NOUT, - "Figure of merit switch", - "(minmax)", - self.data.numerics.minmax, - ) - objf_name = f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"' +class Scan: + """Perform a parameter scan + + Parameters + ---------- + models : + Physics and engineering model objects + solver : + Which solver to use, as specified in solver.py + data : + Data structure object + """ + + def __init__(self, models: Models, solver: str, data: DataStructure): + self.models = models + self.solver = solver + self.data = data - self.data.numerics.objf_name = objf_name + def _run(self, iscan, nsweep, sweep, data): + sh = SolverHandler(self.models, self.solver, data) + # TODO queue the output to avoid race condition (?) + if data.scan.nsweep is not None: + self.write_point_header(iscan) + start_time = time.time() + ifail = sh.run() + end_time = time.time() - start_time + write_output_files(models=self.models, data=data, ifail=ifail, runtime=end_time) + nums = data.numerics + nums.sqsumsq = sum(r**2 for r in nums.rcm[: nums.neqns]) ** 0.5 + + show_errors(constants.NOUT) + + logging_model_handler.clear_logs() + optimisation_output(data) + constraints.constraints_output(data, self.solver) + + return ScanRes(iscan, ifail, sh) + + def _set_v_x_label(self, iscan: list[int]): + sv = [ + self.scan_select(self.data.scan.nsweep, self.data.scan.sweep, isc) + for isc in iscan + ] + self.data.globals.vlabel = [s.fname for s in sv] + self.data.globals.xlabel = [s.data.description for s in sv] - process_output.ovarre( - constants.NOUT, - "Objective function name", - "(objf_name)", - self.data.numerics.objf_name, - ) - process_output.ovarre( - constants.NOUT, - "Normalised objective function", - "(norm_objf)", - self.data.numerics.norm_objf, - "OP ", - ) + def write_point_header(self, iscan): + self._set_v_x_label(iscan) - process_output.ovarre( - constants.NOUT, - "Square root of the sum of squares of the constraint residuals", - "(sqsumsq)", - self.data.numerics.sqsumsq, - "OP ", - ) - if self.solver != "fsolve": - process_output.ovarre( - constants.NOUT, - "VMCON convergence parameter", - "(convergence_parameter)", - self.data.globals.convergence_parameter, - "OP ", - ) - process_output.ovarre( - constants.NOUT, - "Number of optimising solver iterations", - "(nviter)", - self.data.numerics.nviter, - "OP ", - ) process_output.oblnkl(constants.NOUT) + process_output.oblnkl(constants.MFILE) - if self.solver == "fsolve": - if ifail == 1: - msg = "PROCESS has solved using fsolve." - else: - msg = "PROCESS failed to solve using fsolve." - process_output.write( - constants.NOUT, - f"{msg}\n", - ) - else: - if ifail == 1: - string1 = "PROCESS has successfully optimised" - else: - string1 = "PROCESS has failed to optimise" - - string2 = "minimise" if self.data.numerics.minmax > 0 else "maximise" - - process_output.write( - constants.NOUT, - f"{string1} the optimisation parameters to {string2} " - f"the objective function: {objf_name}\n", - ) - - written_warning = False - - # Output optimisation parameters - solution_vector_table = [] - for i in range(self.data.numerics.nvar): - self.data.numerics.xcs[i] = ( - self.data.numerics.xcm[i] * self.data.numerics.scafc[i] - ) - - name = self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1] - solution_vector_table.append([ - name, - self.data.numerics.xcs[i], - self.data.numerics.xcm[i], - ]) - - xminn = 1.01 * self.data.numerics.itv_scaled_lower_bounds[i] - xmaxx = 0.99 * self.data.numerics.itv_scaled_upper_bounds[i] - - # Write to output file if close to optimisation parameter bounds - if self.data.numerics.xcm[i] < xminn or self.data.numerics.xcm[i] > xmaxx: - if not written_warning: - written_warning = True - process_output.ocmmnt( - constants.NOUT, - ( - "Certain operating limits have been reached," - "\n as shown by the following optimisation parameters" - " that are" - "\n at or near to the edge of their prescribed range :\n" - ), - ) - - xcval = self.data.numerics.xcm[i] * self.data.numerics.scafc[i] - - if self.data.numerics.xcm[i] < xminn: - location, bound = "below", "lower" - bounds = self.data.numerics.itv_scaled_lower_bounds - else: - location, bound = "above", "upper" - bounds = self.data.numerics.itv_scaled_upper_bounds - process_output.write( - constants.NOUT, - f" {name:<30}= {xcval} is at or {location} its {bound} bound:" - f" {bounds[i] * self.data.numerics.scafc[i]}", - ) - - # Write optimisation parameters to mfile - process_output.ovarre( - constants.MFILE, - self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1], - f"(itvar{i + 1:03d})", - self.data.numerics.xcs[i], - ) - - if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]: - xnorm = 1.0 - else: - xnorm = min( - max( - ( - self.data.numerics.xcm[i] - - self.data.numerics.itv_scaled_lower_bounds[i] - ) - / ( - self.data.numerics.itv_scaled_upper_bounds[i] - - self.data.numerics.itv_scaled_lower_bounds[i] - ), - 0.0, - ), - 1.0, - ) - - process_output.ovarre( - constants.MFILE, - f"{name} (final value/initial value)", - f"(xcm{i + 1:03d})", - self.data.numerics.xcm[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} (range normalised)", - f"(nitvar{i + 1:03d})", - xnorm, - ) - process_output.ovarre( - constants.MFILE, - f"{name} (upper bound)", - f"(boundu{i + 1:03d})", - self.data.numerics.itv_scaled_upper_bounds[i] - * self.data.numerics.scafc[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} (lower bound)", - f"(boundl{i + 1:03d})", - self.data.numerics.itv_scaled_lower_bounds[i] - * self.data.numerics.scafc[i], - ) - - # Write optimisation parameter headings to output file - process_output.osubhd( - constants.NOUT, "The solution vector is comprised as follows :" - ) process_output.write( constants.NOUT, - tabulate( - solution_vector_table, - headers=["", "Final value", "Final / initial"], - numalign="left", + f"Scan point {iscan} of {np.prod(self.data.scan.isweep)} : \n".join( + f"{v} = {self.data.scan.sweep[iscan[no] - 1]}" + for no, v in enumerate(self.data.globals.vlabel) ), ) + process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan) - process_output.osubhd( - constants.NOUT, - "The following equality constraint residues should be close to zero:", - ) - - con1, con2, err, _, lab = constraints.constraint_eqns( - self.data.numerics.neqns + self.data.numerics.nineqns, -1, self.data - ) - - # Write equality constraints to mfile - equality_constraint_table = [] - for i in range(self.data.numerics.neqns): - name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1] - - equality_constraint_table.append([ - name, - "=", - f"{con2[i]} {lab[i]}", - f"{err[i]} {lab[i]}", - con1[i], - ]) - process_output.ovarre( - constants.MFILE, - f"{name:<33} normalised residue", - f"(eq_con{self.data.numerics.icc[i]:03d})", - con1[i], - ) - - process_output.ovarre( - constants.MFILE, - f"{name:<33} residual", - f"(res_eq_con{self.data.numerics.icc[i]:03d})", - err[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} constraint value", - f"(val_eq_con{self.data.numerics.icc[i]:03d})", - con2[i], - ) - - process_output.ovarre( - constants.MFILE, - f"{name} units", - f"(eq_units_con{self.data.numerics.icc[i]:03d})", - f"'{lab[i]}'", + print( + f"Starting scan point {iscan}: {self.data.globals.xlabel}, \n".join( + f"{v} = {self.data.scan.sweep[iscan[no] - 1]}" + for no, v in enumerate(self.data.globals.vlabel) ) - - # Write equality constraints to output file - process_output.write( - constants.NOUT, - tabulate( - equality_constraint_table, - headers=[ - "", - "", - "Physical constraint", - "Constraint residue", - "Normalised residue", - ], - numalign="left", - ), ) - # Write inequality constraints - if self.data.numerics.nineqns > 0: - inequality_constraint_table = [] - # Inequalities not necessarily satisfied when evaluating - process_output.osubhd( - constants.NOUT, - "Negative inequality constraint (normalised) residuals " - "indicate a constraint is satisfied.", - ) - if self.solver == "fsolve": - process_output.osubhd( - constants.NOUT, - "This MFile was produced via an evaluation, not an optimisation, " - "and so the constraints might be violated.", - ) - - for i in range( - self.data.numerics.neqns, - self.data.numerics.neqns + self.data.numerics.nineqns, - ): - name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1] - constraint = constraints.ConstraintManager.evaluate_constraint( - int(self.data.numerics.icc[i]), self.data - ) - - inequality_constraint_table.append([ - name, - f"{constraint.constraint_value} {constraint.units}", - constraint.symbol, - f"{constraint.constraint_bound} {constraint.units}", - f"{constraint.residual} {constraint.units}", - f"{constraint.normalised_residual}", - ]) - process_output.ovarre( - constants.MFILE, - f"{name} normalised residue", - f"(ineq_con{self.data.numerics.icc[i]:03d})", - -constraint.normalised_residual, - ) - process_output.ovarre( - constants.MFILE, - f"{name} physical value", - f"(ineq_value_con{self.data.numerics.icc[i]:03d})", - constraint.constraint_value, - ) - - process_output.ovarre( - constants.MFILE, - f"{name} symbol", - f"(ineq_symbol_con{self.data.numerics.icc[i]:03d})", - f"'{constraint.symbol}'", - ) - - process_output.ovarre( - constants.MFILE, - f"{name} units", - f"(ineq_units_con{self.data.numerics.icc[i]:03d})", - f"'{constraint.units}'", - ) - - process_output.ovarre( - constants.MFILE, - f"{name} physical bound", - f"(ineq_bound_con{self.data.numerics.icc[i]:03d})", - constraint.constraint_bound, - ) - - process_output.write( - constants.NOUT, - tabulate( - inequality_constraint_table, - headers=[ - "", - "Physical constraint", - "", - "Physical constraint bound", - "Constraint residue", - "Normalised residue", - ], - numalign="left", - ), - ) - - @staticmethod - def verror(ifail: int): - """Routine to print out relevant messages in the case of an - unfeasible result from a VMCON (optimisation) run - - ifail : input integer : error flag - This routine prints out relevant messages in the case of - an unfeasible result from a VMCON (optimisation) run. + def scan_select(self, nsweep, sweep, iscan): + sv = ScanVariables(nsweep) + sv.set(self.data, sweep[iscan - 1]) + return sv - Parameters - ---------- - ifail: int : + def run(self): + """Call a solver over a range of values of one of the variables. + This method calls the optimisation routine VMCON a number of times, by + performing a sweep over a range of values of a particular variable. A + number of output variable values are written to the MFILE.DAT file at + each scan point, for plotting or other post-processing purposes. """ - if ifail == -1: - process_output.ocmmnt(constants.NOUT, "User-terminated execution of VMCON.") - process_output.ocmmnt(constants.IOTTY, "User-terminated execution of VMCON.") - elif ifail == 0: - process_output.ocmmnt( - constants.NOUT, "Improper input parameters to the VMCON routine." - ) - process_output.ocmmnt(constants.NOUT, "PROCESS coding must be checked.") - - process_output.ocmmnt( - constants.IOTTY, "Improper input parameters to the VMCON routine." - ) - process_output.ocmmnt(constants.IOTTY, "PROCESS coding must be checked.") - elif ifail == 2: - process_output.ocmmnt( - constants.NOUT, - "The maximum number of calls has been reached without solution.", - ) - process_output.ocmmnt( - constants.NOUT, - "The code may be stuck in a minimum in the residual space that is " - "significantly above zero.", - ) - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, "There is either no solution possible, or the code" - ) - process_output.ocmmnt( - constants.NOUT, "is failing to escape from a deep local minimum." - ) - process_output.ocmmnt( - constants.NOUT, - "Try changing the variables in IXC, or modify their initial values.", - ) - - process_output.ocmmnt( - constants.IOTTY, - "The maximum number of calls has been reached without solution.", - ) - process_output.ocmmnt( - constants.IOTTY, - "The code may be stuck in a minimum in the residual space that is " - "significantly above zero.", - ) - process_output.oblnkl(constants.NOUT) - process_output.oblnkl(constants.IOTTY) - process_output.ocmmnt( - constants.IOTTY, "There is either no solution possible, or the code" - ) - process_output.ocmmnt( - constants.IOTTY, "is failing to escape from a deep local minimum." - ) - process_output.ocmmnt( - constants.IOTTY, - "Try changing the variables in IXC, or modify their initial values.", - ) - elif ifail == 3: - process_output.ocmmnt( - constants.NOUT, "The line search required the maximum of 10 calls." - ) - process_output.ocmmnt( - constants.NOUT, "A feasible solution may be difficult to achieve." - ) - process_output.ocmmnt( - constants.NOUT, "Try changing or adding variables to IXC." - ) - - process_output.ocmmnt( - constants.IOTTY, "The line search required the maximum of 10 calls." - ) - process_output.ocmmnt( - constants.IOTTY, "A feasible solution may be difficult to achieve." - ) - process_output.ocmmnt( - constants.IOTTY, "Try changing or adding variables to IXC." - ) - elif ifail == 4: - process_output.ocmmnt( - constants.NOUT, "An uphill search direction was found." - ) - process_output.ocmmnt( - constants.NOUT, "Try changing the equations in ICC, or" - ) - process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.") - - process_output.ocmmnt( - constants.IOTTY, "An uphill search direction was found." - ) - process_output.ocmmnt( - constants.IOTTY, "Try changing the equations in ICC, or" - ) - process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.") - elif ifail == 5: - process_output.ocmmnt( - constants.NOUT, "The quadratic programming technique was unable to" - ) - process_output.ocmmnt(constants.NOUT, "find a feasible point.") - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, "Try changing or adding variables to IXC, or modify" - ) - process_output.ocmmnt( - constants.NOUT, - "their initial values (especially if only 1 optimisation", - ) - process_output.ocmmnt(constants.NOUT, "iteration was performed).") - - process_output.ocmmnt( - constants.IOTTY, "The quadratic programming technique was unable to" - ) - process_output.ocmmnt(constants.IOTTY, "find a feasible point.") - process_output.oblnkl(constants.IOTTY) - process_output.ocmmnt( - constants.IOTTY, "Try changing or adding variables to IXC, or modify" - ) - process_output.ocmmnt( - constants.IOTTY, - "their initial values (especially if only 1 optimisation", - ) - process_output.ocmmnt(constants.IOTTY, "iteration was performed).") - elif ifail == 6: - process_output.ocmmnt( - constants.NOUT, "The quadratic programming technique was restricted" - ) - process_output.ocmmnt( - constants.NOUT, "by an artificial bound, or failed due to a singular" - ) - process_output.ocmmnt(constants.NOUT, "matrix.") - process_output.ocmmnt( - constants.NOUT, "Try changing the equations in ICC, or" - ) - process_output.ocmmnt(constants.NOUT, "adding new variables to IXC.") - - process_output.ocmmnt( - constants.IOTTY, "The quadratic programming technique was restricted" - ) - process_output.ocmmnt( - constants.IOTTY, "by an artificial bound, or failed due to a singular" - ) - process_output.ocmmnt(constants.IOTTY, "matrix.") - process_output.ocmmnt( - constants.IOTTY, "Try changing the equations in ICC, or" - ) - process_output.ocmmnt(constants.IOTTY, "adding new variables to IXC.") + # vectorise running of self._run + if self.data.scan.nsweep is not None: + for d, n, v in ( + ("Number of scan points", "(isweep)", self.data.scan.isweep), + ("Scanning variable number", "(nsweep)", self.data.scan.nsweep), + ): + process_output.ovarin(constants.MFILE, d, n, v) - def scan_1d(self): - """Run a 1-D scan.""" - # initialise dict which will contain ifail values for each scan point - scan_1d_ifail_dict = {} + # TODO copy of self.data for each vectorised run (?) + scan_res = np.vectorise(self._run)( + self.data.scan.isweep, self.data.scan.nsweep, self.data.scan.sweep, self.data + ) - for iscan in range(1, self.data.scan.isweep + 1): - self.scan_1d_write_point_header(iscan) - start_time = time.time() - ifail = self.doopt() - scan_1d_ifail_dict[iscan] = ifail - write_output_files( - models=self.models, - data=self.data, - ifail=ifail, - runtime=time.time() - start_time, - ) + if self.data.scan.nsweep is not None: + self.summary(scan_res) - show_errors(constants.NOUT) - logging_model_handler.clear_logs() + def summary(self, scan_res): + print("Scan Convergence Summary\n") + sweep_values = self.data.scan.sweep + nsweep_var = [ScanVariables(nsw) for nsw in self.data.scan.nsweep] - # outvar now contains results - self.scan_1d_write_plot(self.data.scan) - print("Scan Convergence Summary \n") - sweep_values = self.data.scan.sweep[: self.data.scan.isweep] - nsweep_var_name, _ = self.scan_select( - self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep - ) + conv_list = [] converged_count = 0 - # offsets for aligning the converged/unconverged column - max_sweep_value_length = len(str(np.max(sweep_values)).replace(".", "")) - offsets = [ - max_sweep_value_length - len(str(sweep_val).replace(".", "")) - for sweep_val in sweep_values - ] - for iscan in range(1, self.data.scan.isweep + 1): - if scan_1d_ifail_dict[iscan] == 1: + conv_str = "\u001b[3{}CONVERGED \u001b[0m" + for no, sr in enumerate(scan_res): + if sr.ifail == 1: converged_count += 1 - print( - f"Scan {iscan:02d}: {nsweep_var_name} = {sweep_values[iscan - 1]} " - + " " * offsets[iscan - 1] - + "\u001b[32mCONVERGED \u001b[0m" - ) + conv = conv_str.format("2m") else: - print( - f"Scan {iscan:02d}: {nsweep_var_name} = {sweep_values[iscan - 1]} " - + " " * offsets[iscan - 1] - + "\u001b[31mUNCONVERGED \u001b[0m" - ) - converged_percentage = converged_count / self.data.scan.isweep * 100 - print(f"\nConvergence Percentage: {converged_percentage:.2f}%") - - def scan_2d(self): - """Run a 2-D scan.""" - # Initialise intent(out) arrays - self.scan_2d_init(self.data.scan) - iscan = 1 - - # initialise array which will contain ifail values for each scan point - scan_2d_ifail_list = np.zeros( - (NOUTVARS, IPNSCNS), - dtype=np.float64, - order="F", - ) - for iscan_1 in range(1, self.data.scan.isweep + 1): - for iscan_2 in range(1, self.data.scan.isweep_2 + 1): - self.scan_2d_write_point_header(iscan, iscan_1, iscan_2) - start_time = time.time() - ifail = self.doopt() - write_output_files( - models=self.models, - data=self.data, - ifail=ifail, - runtime=time.time() - start_time, - ) - - show_errors(constants.NOUT) - logging_model_handler.clear_logs() - scan_2d_ifail_list[iscan_1][iscan_2] = ifail - iscan += 1 + conv = conv_str.format("1mUN") + conv_list.append([ + "{sr.iscan:02d}", + nsweep_var[no].fname, + sweep_values[sr.iscan], + conv, + ]) - print("Scan Convergence Summary\n") - sweep_1_values = self.data.scan.sweep[: self.data.scan.isweep] - sweep_2_values = self.data.scan.sweep_2[: self.data.scan.isweep_2] - nsweep_var_name, _ = self.scan_select( - self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep - ) - nsweep_2_var_name, _ = self.scan_select( - self.data.scan.nsweep_2, self.data.scan.sweep_2, self.data.scan.isweep_2 - ) - converged_count = 0 - scan_point = 1 - # offsets for aligning the converged/unconverged column - max_sweep1_value_length = len(str(np.max(sweep_1_values)).replace(".", "")) - max_sweep2_value_length = len(str(np.max(sweep_2_values)).replace(".", "")) - offsets = np.zeros( - (self.data.scan.isweep, self.data.scan.isweep_2), dtype=int, order="F" + print( + tabulate(conv_list, headers=["Iscan", "Sweep Var", "Sweep Val", "Converged"]) ) - for count1, sweep1 in enumerate(sweep_1_values): - for count2, sweep2 in enumerate(sweep_2_values): - offsets[count1][count2] = ( - max_sweep1_value_length - - len(str(sweep1).replace(".", "")) - + max_sweep2_value_length - - len(str(sweep2).replace(".", "")) - ) - for iscan_1 in range(1, self.data.scan.isweep + 1): - for iscan_2 in range(1, self.data.scan.isweep_2 + 1): - if scan_2d_ifail_list[iscan_1][iscan_2] == 1: - converged_count += 1 - print( - ( - f"Scan {scan_point:02d}: ({nsweep_var_name} = " - f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var_name} " - f"= {sweep_2_values[iscan_2 - 1]}) " - ) - + " " * offsets[iscan_1 - 1][iscan_2 - 1] - + "\u001b[32mCONVERGED \u001b[0m" - ) - scan_point += 1 - else: - print( - ( - f"Scan {scan_point:02d}: ({nsweep_var_name} = " - f"{sweep_1_values[iscan_1 - 1]}, {nsweep_2_var_name} = " - f"{sweep_2_values[iscan_2 - 1]}) " - ) - + " " * offsets[iscan_1 - 1][iscan_2 - 1] - + "\u001b[31mUNCONVERGED \u001b[0m" - ) - scan_point += 1 - converged_percentage = ( - converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100 - ) + converged_percentage = converged_count / np.prod(self.data.scan.isweep) * 100 print(f"\nConvergence Percentage: {converged_percentage:.2f}%") - @staticmethod - def scan_2d_init(scan_data: ScanData): - process_output.ovarre( - constants.MFILE, - "Number of first variable scan points", - "(isweep)", - scan_data.isweep, - ) - process_output.ovarre( - constants.MFILE, - "Number of second variable scan points", - "(isweep_2)", - scan_data.isweep_2, - ) - process_output.ovarre( - constants.MFILE, - "Scanning first variable number", - "(nsweep)", - scan_data.nsweep, - ) - process_output.ovarre( - constants.MFILE, - "Scanning second variable number", - "(nsweep_2)", - scan_data.nsweep_2, - ) - process_output.ovarre( - constants.MFILE, - "Scanning second variable number", - "(nsweep_2)", - scan_data.nsweep_2, - ) - process_output.ovarre( - constants.MFILE, - "Scanning second variable number", - "(nsweep_2)", - scan_data.nsweep_2, - ) - def scan_1d_write_point_header(self, iscan: int): - self.data.globals.iscan_global = iscan - self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( - self.data.scan.nsweep, self.data.scan.sweep, iscan - ) +def optimisation_output(data: DataStructure): + nums = data.numerics - process_output.oblnkl(constants.NOUT) - process_output.ostars(constants.NOUT, 110) + written_warning = False - process_output.write( - constants.NOUT, - f"***** Scan point {iscan} of {self.data.scan.isweep} : " - f"{self.data.globals.xlabel}" - f", {self.data.globals.vlabel} = {self.data.scan.sweep[iscan - 1]} " - "*****", - ) - process_output.ostars(constants.NOUT, 110) - process_output.oblnkl(constants.MFILE) - process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan) - - print( - f"Starting scan point {iscan} of {self.data.scan.isweep} : " - f"{self.data.globals.xlabel} , {self.data.globals.vlabel}" - f" = {self.data.scan.sweep[iscan - 1]}" - ) - - def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): - iscan_r = self.data.scan.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2 - - # Makes iscan available globally (read-only) - self.data.globals.iscan_global = iscan - - self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( - self.data.scan.nsweep, self.data.scan.sweep, iscan_1 - ) - self.data.globals.vlabel_2, self.data.globals.xlabel_2 = self.scan_select( - self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r - ) + # Output optimisation parameters + solution_vector_table = [] + for i in range(nums.nvar): + nums.xcs[i] = nums.xcm[i] * nums.scafc[i] - process_output.oblnkl(constants.NOUT) - process_output.ostars(constants.NOUT, 110) + name = nums.lablxc[nums.ixc[i] - 1] + solution_vector_table.append([name, nums.xcs[i], nums.xcm[i]]) - process_output.write( - constants.NOUT, - f"***** 2D Scan point {iscan} of " - f"{self.data.scan.isweep * self.data.scan.isweep_2} : " - f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]} and" - f" {self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} " - "*****", - ) - process_output.ostars(constants.NOUT, 110) - process_output.oblnkl(constants.MFILE) - process_output.ovarre(constants.MFILE, "Scan point number", "(iscan)", iscan) + xminn = 1.01 * nums.itv_scaled_lower_bounds[i] + xmaxx = 0.99 * nums.itv_scaled_upper_bounds[i] - print( - f"Starting scan point {iscan}: {self.data.globals.xlabel}, " - f"{self.data.globals.vlabel} = {self.data.scan.sweep[iscan_1 - 1]}" - f" and {self.data.globals.xlabel_2}, " - f"{self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} " - ) + # Write to output file if close to optimisation parameter bounds + if nums.xcm[i] < xminn or nums.xcm[i] > xmaxx: + if not written_warning: + written_warning = True + process_output.ocmmnt( + constants.NOUT, + ( + "Certain operating limits have been reached," + "\n as shown by the following optimisation parameters that are" + "\n at or near to the edge of their prescribed range :\n" + ), + ) - return iscan_r + xcval = nums.xcm[i] * nums.scafc[i] - @staticmethod - def scan_1d_write_plot(scan_data: ScanData): - if scan_data.first_call_1d: - process_output.ovarre( - constants.MFILE, - "Number of scan points", - "(isweep)", - scan_data.isweep, - ) - process_output.ovarre( - constants.MFILE, - "Scanning variable number", - "(nsweep)", - scan_data.nsweep, + if nums.xcm[i] < xminn: + location, bound = "below", "lower" + bounds = nums.itv_scaled_lower_bounds + else: + location, bound = "above", "upper" + bounds = nums.itv_scaled_upper_bounds + process_output.write( + constants.NOUT, + f" {name:<30}= {xcval} is at or {location} its {bound} bound:" + f" {bounds[i] * nums.scafc[i]}", + ) + + xnorm = ( + 1.0 + if nums.boundu[i] == nums.boundl[i] + else min( + max( + (nums.xcm[i] - nums.itv_scaled_lower_bounds[i]) + / ( + nums.itv_scaled_upper_bounds[i] - nums.itv_scaled_lower_bounds[i] + ), + 0.0, + ), + 1.0, ) + ) - scan_data.first_call_1d = False - - def scan_select(self, nwp, swp, iscn): - match nwp: - case 1: - self.data.physics.aspect = swp[iscn - 1] - case 2: - self.data.divertor.pflux_div_heat_load_max_mw = swp[iscn - 1] - case 3: - self.data.constraints.p_plant_electric_net_required_mw = swp[iscn - 1] - case 4: - self.data.physics.hfact = swp[iscn - 1] - case 5: - self.data.tfcoil.j_tf_coil_full_area = swp[iscn - 1] - case 6: - self.data.constraints.pflux_fw_neutron_max_mw = swp[iscn - 1] - case 7: - self.data.physics.beamfus0 = swp[iscn - 1] - case 9: - self.data.physics.temp_plasma_electron_vol_avg_kev = swp[iscn - 1] - case 10: - self.data.numerics.boundu[14] = swp[iscn - 1] - case 11: - self.data.physics.beta_norm_max = swp[iscn - 1] - case 12: - self.data.current_drive.f_c_plasma_bootstrap_max = swp[iscn - 1] - case 13: - self.data.numerics.boundu[9] = swp[iscn - 1] - case 16: - self.data.physics.rmajor = swp[iscn - 1] - case 17: - self.data.constraints.b_tf_inboard_max = swp[iscn - 1] - case 18: - self.data.constraints.eta_cd_norm_hcd_primary_max = swp[iscn - 1] - case 19: - self.data.numerics.boundl[15] = swp[iscn - 1] - case 20: - self.data.constraints.t_burn_min = swp[iscn - 1] - case 22: - if self.data.costs.i_plant_availability == 1: - raise ProcessValueError( - "Do not scan f_t_plant_available if i_plant_availability=1" - ) - self.data.costs.f_t_plant_available = swp[iscn - 1] - case 24: - self.data.constraints.p_fusion_total_max_mw = swp[iscn - 1] - case 25: - self.data.physics.kappa = swp[iscn - 1] - case 26: - self.data.physics.triang = swp[iscn - 1] - case 27: - self.data.constraints.tbrmin = swp[iscn - 1] - case 28: - self.data.physics.b_plasma_toroidal_on_axis = swp[iscn - 1] - case 29: - self.data.impurity_radiation.radius_plasma_core_norm = swp[iscn - 1] - case 31: - self.data.constraints.f_t_alpha_energy_confinement_min = swp[iscn - 1] - case 32: - self.data.numerics.epsvmc = swp[iscn - 1] - case 38: - self.data.numerics.boundu[128] = swp[iscn - 1] - case 39: - self.data.numerics.boundu[130] = swp[iscn - 1] - case 40: - self.data.numerics.boundu[134] = swp[iscn - 1] - case 41: - self.data.build.dr_blkt_outboard = swp[iscn - 1] - case 42: - self.data.impurity_radiation.f_nd_impurity_electrons[8] = swp[iscn - 1] - self.data.impurity_radiation.f_nd_impurity_electron_array[8] = ( - self.data.impurity_radiation.f_nd_impurity_electrons[8] - ) - case 44: - self.data.tfcoil.sig_tf_case_max = swp[iscn - 1] - case 45: - self.data.tfcoil.temp_tf_superconductor_margin_min = swp[iscn - 1] - case 46: - self.data.numerics.boundu[151] = swp[iscn - 1] - case 48: - self.data.tfcoil.n_tf_wp_pancakes = int(swp[iscn - 1]) - case 49: - self.data.tfcoil.n_tf_wp_layers = int(swp[iscn - 1]) - case 50: - self.data.impurity_radiation.f_nd_impurity_electrons[12] = swp[iscn - 1] - self.data.impurity_radiation.f_nd_impurity_electron_array[12] = ( - self.data.impurity_radiation.f_nd_impurity_electrons[12] - ) - case 51: - self.data.physics.f_p_div_lower = swp[iscn - 1] - case 52: - self.data.physics.rad_fraction_sol = swp[iscn - 1] - case 53: - self.data.numerics.boundu[156] = swp[iscn - 1] - case 54: - self.data.tfcoil.b_crit_upper_nbti = swp[iscn - 1] - case 55: - self.data.build.dr_shld_inboard = swp[iscn - 1] - case 56: - self.data.heat_transport.p_cryo_plant_electric_max_mw = swp[iscn - 1] - case 57: - self.data.numerics.boundl[1] = swp[iscn - 1] - case 58: - self.data.build.dr_fw_plasma_gap_inboard = swp[iscn - 1] - case 59: - self.data.build.dr_fw_plasma_gap_outboard = swp[iscn - 1] - case 60: - self.data.tfcoil.sig_tf_wp_max = swp[iscn - 1] - case 61: - self.data.rebco.copperaoh_m2_max = swp[iscn - 1] - case 62: - self.data.pf_coil.j_cs_flat_top_end = swp[iscn - 1] - case 63: - self.data.build.dr_cs = swp[iscn - 1] - case 64: - self.data.pf_coil.f_z_cs_tf_internal = swp[iscn - 1] - case 65: - self.data.cs_fatigue.n_cycle_min = swp[iscn - 1] - case 66: - self.data.pf_coil.f_a_cs_turn_steel = swp[iscn - 1] - case 67: - self.data.cs_fatigue.t_crack_vertical = swp[iscn - 1] - case 68: - self.data.fwbs.inlet_temp_liq = swp[iscn - 1] - case 69: - self.data.fwbs.outlet_temp_liq = swp[iscn - 1] - case 70: - self.data.fwbs.blpressure_liq = swp[iscn - 1] - case 71: - self.data.fwbs.n_liq_recirc = swp[iscn - 1] - case 72: - self.data.fwbs.bz_channel_conduct_liq = swp[iscn - 1] - case 73: - self.data.fwbs.pnuc_fw_ratio_dcll = swp[iscn - 1] - case 74: - self.data.fwbs.f_nuc_pow_bz_struct = swp[iscn - 1] - case 75: - self.data.fwbs.dx_fw_module = swp[iscn - 1] - case 76: - self.data.heat_transport.eta_turbine = swp[iscn - 1] - case 77: - self.data.costs.startupratio = swp[iscn - 1] - case 78: - self.data.costs.fkind = swp[iscn - 1] - case 79: - self.data.current_drive.eta_ecrh_injector_wall_plug = swp[iscn - 1] - case 80: - self.data.tfcoil.fcoolcp = swp[iscn - 1] - case 81: - self.data.tfcoil.n_tf_coil_turns = swp[iscn - 1] - case _: - raise ProcessValueError("Illegal scan variable number", nwp=nwp) + # Write optimisation parameters to mfile + for d, var, v in ( + (nums.lablxc[nums.ixc[i] - 1], f"(itvar{i + 1:03d})", nums.xcs[i]), + (f"{name} (final value/initial value)", f"(xcm{i + 1:03d})", nums.xcm[i]), + (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm), + ( + f"{name} (upper bound)", + f"(boundu{i + 1:03d})", + nums.itv_scaled_upper_bounds[i] * nums.scafc[i], + ), + ( + f"{name} (lower bound)", + f"(boundl{i + 1:03d})", + nums.itv_scaled_lower_bounds[i] * nums.scafc[i], + ), + ): + process_output.ovarre(constants.MFILE, d, var, v) - return ScanVariables(int(nwp)).value + # Write optimisation parameter headings to output file + process_output.osubhd( + constants.NOUT, "The solution vector is comprised as follows :" + ) + process_output.write( + constants.NOUT, + tabulate( + solution_vector_table, + headers=["", "Final value", "Final / initial"], + numalign="left", + ), + ) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 3c2406507e..d681c1440f 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -4,8 +4,9 @@ from typing import ClassVar, Literal import numpy as np +from tabulate import tabulate -from process.core import constants +from process.core import constants, process_output from process.core.exceptions import ProcessError, ProcessValueError from process.core.model import DataStructure from process.data_structure.build_variables import TFCSRadialConfiguration @@ -1953,3 +1954,123 @@ def constraint_eqns(m: int, ieqn: int, data: DataStructure): units.append(tmp_units) return np.array(cc), np.array(con), np.array(err), symbol, units + + +def constraints_output(data: DataStructure, solver_name: str): + nums = data.numerics + + process_output.osubhd( + constants.NOUT, + "The following equality constraint residues should be close to zero:", + ) + + con1, con2, err, _, lab = constraint_eqns(nums.neqns + nums.nineqns, -1, data) + + # Write equality constraints to mfile + equality_constraint_table = [] + for i in range(nums.neqns): + name = nums.lablcc[nums.icc[i] - 1] + + equality_constraint_table.append([ + name, + "=", + f"{con2[i]} {lab[i]}", + f"{err[i]} {lab[i]}", + con1[i], + ]) + + for d, var, v in ( + (f"{name:<33} normalised residue", f"(eq_con{nums.icc[i]:03d})", con1[i]), + (f"{name:<33} residual", f"(res_eq_con{nums.icc[i]:03d})", err[i]), + (f"{name} constraint value", f"(val_eq_con{nums.icc[i]:03d})", con2[i]), + (f"{name} units", f"(eq_units_con{nums.icc[i]:03d})", f"'{lab[i]}'"), + ): + process_output.ovarre(constants.MFILE, d, var, v) + + # Write equality constraints to output file + process_output.write( + constants.NOUT, + tabulate( + equality_constraint_table, + headers=[ + "", + "", + "Physical constraint", + "Constraint residue", + "Normalised residue", + ], + numalign="left", + ), + ) + + # Write inequality constraints + if nums.nineqns > 0: + inequality_constraint_table = [] + # Inequalities not necessarily satisfied when evaluating + process_output.osubhd( + constants.NOUT, + "Negative inequality constraint (normalised) residuals indicate a constraint is satisfied.", + ) + if solver_name == "fsolve": + process_output.osubhd( + constants.NOUT, + "This MFile was produced via an evaluation, not an optimisation, and so the constraints " + "might be violated.", + ) + + for i in range( + nums.neqns, + nums.neqns + nums.nineqns, + ): + name = nums.lablcc[nums.icc[i] - 1] + constraint = ConstraintManager.evaluate_constraint(int(nums.icc[i]), data) + + inequality_constraint_table.append([ + name, + f"{constraint.constraint_value} {constraint.units}", + constraint.symbol, + f"{constraint.constraint_bound} {constraint.units}", + f"{constraint.residual} {constraint.units}", + f"{constraint.normalised_residual}", + ]) + + for d, var, v in ( + ( + "normalised residue", + f"(ineq_con{nums.icc[i]:03d})", + -constraint.normalised_residual, + ), + ( + "physical value", + f"(ineq_value_con{nums.icc[i]:03d})", + constraint.constraint_value, + ), + ( + "symbol", + f"(ineq_symbol_con{nums.icc[i]:03d})", + f"'{constraint.symbol}'", + ), + ("units", f"(ineq_units_con{nums.icc[i]:03d})", f"'{constraint.units}'"), + ( + "physical bound", + f"(ineq_bound_con{nums.icc[i]:03d})", + constraint.constraint_bound, + ), + ): + process_output.ovarre(constants.MFILE, f"{name} {d}", var, v) + + process_output.write( + constants.NOUT, + tabulate( + inequality_constraint_table, + headers=[ + "", + "Physical constraint", + "", + "Physical constraint bound", + "Constraint residue", + "Normalised residue", + ], + numalign="left", + ), + ) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index ac79dc2789..ce6b8242bd 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -311,10 +311,7 @@ def load_iteration_variables(data): # warn of the iteration variable is also a scan variable because this will cause # the optimiser and scan to overwrite the same variable and conflict - if iteration_variable.name in { - data.globals.vlabel, - data.globals.vlabel_2, - }: + if iteration_variable.name in data.globals.vlabel: logger.critical( ( "The sweep variable is also an iteration variable and will be " diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index a38e6aecf6..a6ac114521 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -16,6 +16,7 @@ ) from scipy.optimize import fsolve +from process.core import constants, process_output from process.core.exceptions import ProcessValueError from process.core.model import DataStructure from process.core.solver.evaluators import Evaluators @@ -283,6 +284,65 @@ def _ineq_cons_satisfied( return self.info + def verror(self): + """Routine to print out relevant messages in the case of an + unfeasible result from a VMCON (optimisation) run + + This routine prints out relevant messages in the case of + an unfeasible result from a VMCON (optimisation) run. + + Parameters + ---------- + ifail: int : + + """ + strings = "\n".join( + { + -1: ("User-terminated execution of VMCON.",), + 0: ( + "Improper input parameters to the VMCON routine.", + "PROCESS coding must be checked.", + ), + 2: ( + "The maximum number of calls has been reached without solution.", + ( + "The code may be stuck in a minimum in the residual space that" + " is significantly above zero.\n" + ), + "There is either no solution possible, or the code", + "is failing to escape from a deep local minimum.", + "Try changing the variables in IXC, or modify their initial values.", + ), + 3: ( + "The line search required the maximum of 10 calls.", + "A feasible solution may be difficult to achieve.", + "Try changing or adding variables to IXC.", + ), + 4: ( + "An uphill search direction was found.", + "Try changing the equations in ICC, or", + "adding new variables to IXC.", + ), + 5: ( + "The quadratic programming technique was unable to", + "find a feasible point.\n", + "Try changing or adding variables to IXC, or modify", + "their initial values (especially if only 1 optimisation", + "iteration was performed).", + ), + 6: ( + "The quadratic programming technique was restricted", + "by an artificial bound, or failed due to a singular", + "matrix.", + "Try changing the equations in ICC, or", + "adding new variables to IXC.", + ), + }.get(self.info, "Unknown Error code") + ) + + process_output.ocmmnt(constants.NOUT, strings) + print(strings) + class VmconBounded(Vmcon): """A solver that uses VMCON but checks x is in bounds before running""" diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index b457e25ca9..cb5f962212 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -1,9 +1,16 @@ +import logging +from contextlib import contextmanager + +from process.core import constants, process_output from process.core.solver.evaluators import Evaluators from process.core.solver.iteration_variables import ( load_iteration_variables, load_scaled_bounds, ) from process.core.solver.solver import get_solver +from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode + +logger = logging.getLogger(__name__) class SolverHandler: @@ -36,14 +43,9 @@ def run(self): # Initialise iteration variables and bounds in Python: relies on Fortran # iteration variables being defined above # Trim maximum size arrays down to actually used size - n = self.data.numerics.nvar - x = self.data.numerics.xcm[:n] - bndl = self.data.numerics.itv_scaled_lower_bounds[:n] - bndu = self.data.numerics.itv_scaled_upper_bounds[:n] - - # Define total number of constraints and equality constraints - m = self.data.numerics.neqns + self.data.numerics.nineqns - meq = self.data.numerics.neqns + x = self.data.numerics.xcm[: self.data.numerics.nvar] + bndl = self.data.numerics.itv_scaled_lower_bounds[: self.data.numerics.nvar] + bndu = self.data.numerics.itv_scaled_upper_bounds[: self.data.numerics.nvar] # Evaluators() calculates the objective and constraint functions and # their gradients for a given vector x @@ -54,30 +56,18 @@ def run(self): self.solver.set_evaluators(evaluators) self.solver.set_bounds(bndl, bndu) self.solver.set_opt_params(x) - self.solver.set_constraints(m, meq) + # Define total number of constraints and equality constraints + self.solver.set_constraints( + m=self.data.numerics.neqns + self.data.numerics.nineqns, + meq=self.data.numerics.neqns, + ) ifail = self.solver.solve() # If VMCON optimisation has failed then try altering value of epsfcn if self.solver_name == "vmcon": if ifail != 1: - print("Trying again with new epsfcn") - # epsfcn is only used in evaluators.Evaluators() - # TODO epsfcn could be set in Evaluators instance now, don't need to - # set/unset in self.data.numerics module - self.data.numerics.epsfcn *= 10 # try new larger value - print("new epsfcn = ", self.data.numerics.epsfcn) - - ifail = self.solver.solve() - # First solution attempt failed (ifail != 1): supply ifail value - # to next attempt - self.data.numerics.epsfcn /= 10 # reset value - - if ifail != 1: - print("Trying again with new epsfcn") - self.data.numerics.epsfcn /= 10 # try new smaller value - print("new epsfcn = ", self.data.numerics.epsfcn) - ifail = self.solver.solve() - self.data.numerics.epsfcn *= 10 # reset value + with epsfcn_context(self.data.numerics): + self.solver.solve() # If VMCON has exited with error code 5 try another run using a multiple # of the identity matrix as input for the Hessian b(n,n) @@ -104,3 +94,135 @@ def output(self): # than required, size self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf + + nums = self.data.numerics + + process_output.oheadr(constants.NOUT, "Numerics") + process_output.ocmmnt( + constants.NOUT, + f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'} (optimisation) run.", + ) + ifail = self.solver.info + if ifail != 1: + process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail) + process_output.oheadr( + constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" + ) + print() + + logger.critical("Solver returns with ifail /= 1. %s", ifail) + + if self.solver_name == "vmcon": + self.solver.verror() + + process_output.oblnkl(constants.NOUT) + print() + else: + # Solution found + descr = "consistent" if self.solver == "fsolve" else "feasible" + process_output.ocmmnt( + constants.NOUT, f"and found a {descr} set of parameters." + ) + process_output.oheadr(constants.IOTTY, f"PROCESS found a {descr} solution") + process_output.oblnkl(constants.NOUT) + process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail) + + if nums.sqsumsq >= 1.0e-2: + string = ( + "WARNING: Constraint residues are HIGH; consider re-running\n" + " with lower values of EPSVMC to confirm convergence...\n" + " (should be able to get down to about 1.0E-8 okay)\n" + ) + process_output.ocmmnt(constants.NOUT, ("\n" + string)) + print(string) + + logger.warning(f"High final constraint residues. {nums.sqsumsq=}") + + for d, var, v in ( + ("Number of iteration variables", "(nvar)", nums.nvar), + ( + "Number of constraints (total)", + "(neqns+nineqns)", + nums.neqns + nums.nineqns, + ), + ("Optimisation switch", "(ioptimz)", nums.ioptimz), + ): + process_output.ovarin(constants.NOUT, d, var, v) + + process_output.ocmmnt( + constants.NOUT, + f" {PROCESSRunMode(nums.ioptimz).description}", + ) + + # Objective function output: none for fsolve + if self.solver_name != "fsolve": + process_output.ovarin( + constants.NOUT, + "Figure of merit switch", + "(minmax)", + nums.minmax, + ) + + nums.objf_name = f'"{FiguresOfMerit(abs(nums.minmax)).description}"' + + for d, var, v, o in ( + ("Objective function name", "(objf_name)", nums.objf_name, ""), + ("Normalised objective function", "(norm_objf)", nums.norm_objf, "OP "), + ( + "VMCON convergence parameter", + "(convergence_parameter)", + self.data.globals.convergence_parameter, + "OP ", + ), + ( + "Number of optimising solver iterations", + "(nviter)", + nums.nviter, + "OP ", + ), + ( + "Square root of the sum of squares of the constraint residuals", + "(sqsumsq)", + nums.sqsumsq, + "OP ", + ), + ): + process_output.ovarre(constants.NOUT, d, var, v, o) + + process_output.oblnkl(constants.NOUT) + + if self.solver_name == "fsolve": + process_output.write( + constants.NOUT, + "PROCESS has solved using fsolve.\n" + if ifail == 1 + else "PROCESS failed to solve using fsolve.\n", + ) + else: + process_output.write( + constants.NOUT, + ( + ( + "PROCESS has successfully optimised" + if ifail == 1 + else "PROCESS has failed to optimise" + ) + + " the optimisation parameters to" + + ("minimise" if nums.minmax > 0 else "maximise") + + f" the objective function: {nums.objf_name}\n" + ), + ) + + +@contextmanager +def epsfcn_context(numerics): + print("Trying again with new epsfcn") + # epsfcn is only used in evaluators.Evaluators() + # TODO epsfcn could be set in Evaluators instance now, don't need to + # set/unset in numerics module + numerics.epsfcn *= 10 # try new larger value + print("new epsfcn = ", numerics.epsfcn) + try: + yield + finally: + numerics.epsfcn /= 10 # reset value diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py index f3e9849c3e..100198ee33 100644 --- a/process/data_structure/global_variables.py +++ b/process/data_structure/global_variables.py @@ -1,4 +1,4 @@ -from dataclasses import dataclass +from dataclasses import dataclass, field @dataclass(slots=True) @@ -18,19 +18,11 @@ class GlobalData: output_prefix: str = "" """Output file path prefix""" - xlabel: str = "" - """Scan parameter description label""" + xlabel: list[str] = field(default_factory=lambda: [""]) + """Scan parameters description label""" - vlabel: str = "" - """Scan value name label""" - - xlabel_2: str = "" - """Scan parameter description label (2nd dimension)""" - - vlabel_2: str = "" - """Scan value name label (2nd dimension)""" - - iscan_global: int = 0 + vlabel: list[str] = field(default_factory=lambda: [""]) + """Scan values name label""" convergence_parameter: float = 0.0 """VMCON convergence parameter 'sum'""" diff --git a/process/data_structure/scan_variables.py b/process/data_structure/scan_variables.py index 04aeb2dc9f..947ce0fd04 100644 --- a/process/data_structure/scan_variables.py +++ b/process/data_structure/scan_variables.py @@ -9,129 +9,54 @@ import numpy as np +from process.core.exceptions import ProcessValueError + IPNSCNS = 1000 """Maximum number of scan points""" -IPNSCNV = 81 -"""Number of available scan variables""" - - -NOUTVARS = 84 - - @dataclass(slots=True) class ScanData: - scan_dim: int = 1 - """1-D or 2-D scan switch (1=1D, 2=2D)""" - - isweep: int = 0 + isweep: list[int] | int = 1 """Number of scan points to calculate""" - isweep_2: int = 0 - """Number of 2D scan points to calculate""" - - nsweep: int = 1 - """Switch denoting quantity to scan: -
  • 68 `inlet_temp_liq' : Inlet temperature of blanket liquid metal coolant/breeder (K) -
  • 69 `outlet_temp_liq' : Outlet temperature of blanket liquid metal coolant/breeder (K) -
  • 70 `blpressure_liq' : Blanket liquid metal breeder/coolant pressure (Pa) -
  • 71 `n_liq_recirc' : Selected number of liquid metal breeder recirculations per day -
  • 72 `bz_channel_conduct_liq' : Conductance of liquid metal breeder duct walls (A V-1 m-1) -
  • 73 `pnuc_fw_ratio_dcll' : Ratio of FW nuclear power as fraction of total (FW+BB) -
  • 74 `f_nuc_pow_bz_struct' : Fraction of BZ power cooled by primary coolant for dual-coolant blanket -
  • 75 dx_fw_module : pitch of first wall cooling channels (m) -
  • 76 eta_turbine : Thermal conversion eff. -
  • 77 startupratio : Gyrotron redundancy -
  • 78 fkind : Multiplier for Nth of a kind costs -
  • 79 eta_ecrh_injector_wall_plug : ECH wall plug to injector efficiency - """ - - nsweep_2: int = 3 - """nsweep_2 /3/ : switch denoting quantity to scan for 2D scan:""" - - sweep: list[float] = field( - default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64) - ) - """sweep(IPNSCNS) /../: actual values to use in scan""" + nsweep: list[int] | int | None = None + """Switch denoting quantity to scan - sweep_2: list[float] = field( - default_factory=lambda: np.zeros(IPNSCNS, dtype=np.float64) - ) - """sweep_2(IPNSCNS) /../: actual values to use in 2D scan""" - - # Vars in subroutines scan_1d and scan_2d requiring re-initialising before - # each new run - - first_call_1d: bool = True + see `process.core.scan.ScanVariables` for available options + """ - first_call_2d: bool = True + sweep: np.ndarray = field(default_factory=lambda: np.zeros(1, dtype=np.float64)) + """Actual values to use in scan""" + + def __post_init__(self): + if isinstance(self.isweep, int): + # avoid old 0 default + self.isweep = [self.isweep or 1] + + if len(self.isweep) > 2 or len(self.sweep.shape) > 2: + raise NotImplementedError("N-D Scans not currently supported") + + if max(self.isweep) > IPNSCNS: + raise ProcessValueError( + "Illegal value of isweep", + isweep=self.isweep, + IPNSCNS=IPNSCNS, + ) + if self.nsweep != len(self.isweep): + raise ValueError( + "Number of sweep variables not equal to scan point dimensions" + ) + if self.sweep.shape != self.isweep: + if self.isweep != 1: + self.isweep = list(self.sweep.shape) + else: + print("Unset sweep values set to zero") + # TODO append to size instead of resetting + self.sweep = np.zeros(self.isweep, dtype=np.float64) + + self.nsweep = np.asarray(self.nsweep, dtype=int) + self.isweep = np.asarray(self.isweep, dtype=int) CREATE_DICTS_FROM_DATACLASS = ScanData diff --git a/process/main.py b/process/main.py index 7d2715c2c9..bcd83f34d2 100644 --- a/process/main.py +++ b/process/main.py @@ -55,6 +55,7 @@ from process.core.model import DataStructure, Model from process.core.process_output import OutputFileManager, oheadr from process.core.scan import Scan +from process.data_structure.blanket_variables import BlktModelTypes from process.data_structure.cost_variables import CostModels from process.data_structure.numerics import PROCESSRunMode from process.models.availability import Availability @@ -110,7 +111,7 @@ from process.models.stellarator.neoclassics import Neoclassics from process.models.stellarator.stellarator import Stellarator from process.models.structure import Structure -from process.models.tfcoil.base import TFCoil +from process.models.tfcoil.base import TFCoil, TFConductorModel from process.models.tfcoil.resistive import ( AluminiumTFCoil, CopperTFCoil, @@ -120,6 +121,7 @@ CICCSuperconductingTFCoil, CROCOSuperconductingTFCoil, SuperconductingTFCoil, + SuperconductingTFTurnType, ) from process.models.vacuum import Vacuum, VacuumVessel from process.models.water_use import WaterUse @@ -434,11 +436,8 @@ def initialise(self): def run_scan(self): """Create scan object if required.""" - # TODO Move this solver logic up to init? - # ioptimz == 1: optimisation if self.data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: pass - # ioptimz == -2: evaluation elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION: # No optimisation: # solve equality (consistency) constraints only using fsolve (HYBRD) @@ -449,6 +448,7 @@ def run_scan(self): "select either 1 (optimise) or -2 (no optimisation)." ) self.scan = Scan(self.models, self.solver, self.data) + self.scan.run() @staticmethod def show_errors(): @@ -802,6 +802,145 @@ def setup_data_structure(self): for model in self.models: model.data = self.data + def write(self, data, _outfile): + """Write the results to the main output file (OUT.DAT). + + Write the program results to a file, in a tidy format. + + Parameters + ---------- + self : process.main.Models + physics and engineering model objects + _outfile : int + Fortran output unit identifier + + """ + # ensure we are capturing warnings that occur in the 'output' stage as these are warnings + # that occur at our solution point. So we clear existing warnings + logging_model_handler.start_capturing() + logging_model_handler.clear_logs() + + # Call stellarator output routine instead if relevant + if data.stellarator.istell != 0: + self.stellarator.output() + return + + # Call IFE output routine instead if relevant + if data.ife.ife != 0: + self.ife.output() + return + + # Costs model + # Cost switch values + # No. | model + # ---- | ------ + # 0 | 1990 costs model + # 1 | 2015 Kovari model + # 2 | Custom model + self.costs.output() + + # Availability model + self.availability.output() + + # Physics model + self.physics.output() + + # Detailed physics, currently only done at final point as values are not used + # by any other functions + self.physics_detailed.output() + + # TODO what is this? Not in caller.py? + self.current_drive.output() + + # Pulsed reactor model + self.pulse.output() + + self.divertor.output() + + # Machine Build Model + self.build.output() + + # Cryostat build + self.cryostat.output() + + # Toroidal field coil copper model + if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER: + self.copper_tf_coil.output() + + # Toroidal field coil superconductor model + if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: + tf_turn_type = SuperconductingTFTurnType( + data.superconducting_tfcoil.i_tf_turn_type + ) + if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT: + self.cicc_sctfcoil.output() + elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR: + self.croco_sctfcoil.output() + else: + raise ValueError( + "Unsupported superconducting TF turn type: " + f"{data.superconducting_tfcoil.i_tf_turn_type}" + ) + + # Toroidal field coil aluminium model + if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM: + self.aluminium_tf_coil.output() + + # Tight aspect ratio machine model + if ( + data.physics.itart == 1 + and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING + ): + self.tfcoil.output() + + # Poloidal field coil model + self.pfcoil.output() + + # Structure Model + self.structure.output() + + # Blanket model + # Blanket switch values + # No. | model + # ---- | ------ + # 1 | CCFE HCPB model + # 2 | KIT HCPB model + # 3 | CCFE HCPB model with Tritium Breeding Ratio calculation + # 4 | KIT HCLL model + # 5 | DCLL model + + self.shield.output() + self.vacuum_vessel.output() + + # First wall geometry + self.fw.output() + + if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB: + # CCFE HCPB model + self.ccfe_hcpb.output() + + elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL: + # DCLL model + self.dcll.output() + + # FISPACT and LOCA model (not used)- removed + + # Power model + self.power.output() + + # Vacuum model + self.vacuum.output() + + # Buildings model + self.buildings.output() + + # Water usage in secondary cooling system + self.water_use.output() + + # stop capturing warnings so that Outfile does not end up with + # a lot of non-model logs + logging_model_handler.stop_capturing() + # setup handlers for writing to terminal (on warnings+) # or writing to the log file (on info+)