diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md index f6cbd1a625..49e01f6ac8 100644 --- a/documentation/source/development/add-vars.md +++ b/documentation/source/development/add-vars.md @@ -105,36 +105,19 @@ objective_function(): After following the instruction to add an input variable, you can make the variable a scan variable by following these steps: -1. Increment the parameter `IPNSCNV` defined in `scan_variables.py` in the data_structure directory, to accommodate the new scanning variable. The incremented value will identify your scan variable. +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. -2. Add a short description of the new scanning variable in the `nsweep` comment in `scan_variables.py`, alongside its identification number. - -3. 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. - -4. 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. - - -`nsweep` comment example: -```fortran - - integer :: nsweep = 1 - !! nsweep /1/ : switch denoting quantity to scan: -``` +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. `SCAN_VARIABLES` case example: ```python - class ScanVariables(Enum): - 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), + class ScanVariables(ScanVariable, Enum): + aspect = (1, SVE.P), + pflux_div_heat_load_max_mw = (2, SVE.D) ... - Bc2_0K: ScanVariable("Bc2(0K)", "GL_NbTi Bc2(0K)", 54), - dr_shld_inboard : ScanVariable("dr_shld_inboard", "Inboard neutronic shield", 55), + b_crit_upper_nbti = (54, SVE.T) + dr_shld_inboard = (55, SVE.B), ``` --------------- diff --git a/process/core/caller.py b/process/core/caller.py index 24ae04f9de..4d9dce666f 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -8,13 +8,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 644078ef83..506ecddf8b 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1099,24 +1099,15 @@ def __post_init__(self): "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 9820a6b967..6251f5c08e 100644 --- a/process/core/io/plot/cli.py +++ b/process/core/io/plot/cli.py @@ -238,10 +238,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 09d7ea0d62..0c906bb734 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,904 +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 varible 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 defintion 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: "fiooic", - 16: "rmajor", - 17: "b_tf_inboard_peak_symmetric", # b_tf_inboard_max the maximum T field upper limit is the scan variable - 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_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", - 57: "b_plasma_toroidal_on_axis", # Genuinly b_plasma_toroidal_on_axis lower bound - 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 imput 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( - f"Warning : Non-convergent scan point : {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} : {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 " - f"{meta[scan_var_name].latex if scan_var_name in meta else {scan_var_name}}", - 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 ''}" - 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): - 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 - 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.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( - 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." - ) + 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}" - ) - 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]}" - - # 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}" - ) + 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 1368eaae7f..e69de29bb2 100644 --- a/process/core/output.py +++ b/process/core/output.py @@ -1,146 +0,0 @@ -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 b2da1471e0..11dbd30104 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,1235 +16,367 @@ 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) - fiooic = ScanVariable("fiooic", "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_alpha_energy_confinement_min = ScanVariable( - "f_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. - - 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() + raise ProcessValueError("Illegal scan variable number", nwp=value) - 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, - ) - - 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." - ) - else: - process_output.ocmmnt( - constants.NOUT, "PROCESS has performed a VMCON (optimisation) run." + "Do not scan f_t_plant_available if i_plant_availability=1" ) - if ifail != 1: - process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail) - process_output.oheadr( - constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" - ) - process_output.oblnkl(constants.IOTTY) - - 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) + 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: - # 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.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail) + 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) + 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) - 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.ovarin( - constants.NOUT, - "Number of iteration variables", - "(nvar)", - self.data.numerics.nvar, - ) - process_output.ovarin( - constants.NOUT, - "Number of constraints (total)", - "(neqns+nineqns)", - self.data.numerics.neqns + self.data.numerics.nineqns, - ) - process_output.ovarin( - 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.ovarin( - 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.ovarst( - 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.ovarin( - 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} 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" + 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} = {sweep_1_values[iscan_1 - 1]}, {nsweep_2_var_name} = {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} = {sweep_1_values[iscan_1 - 1]}, {nsweep_2_var_name} = {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.ovarin( - constants.MFILE, - "Number of first variable scan points", - "(isweep)", - scan_data.isweep, - ) - process_output.ovarin( - constants.MFILE, - "Number of second variable scan points", - "(isweep_2)", - scan_data.isweep_2, - ) - process_output.ovarin( - constants.MFILE, - "Scanning first variable number", - "(nsweep)", - scan_data.nsweep, - ) - process_output.ovarin( - constants.MFILE, - "Scanning second variable number", - "(nsweep_2)", - scan_data.nsweep_2, - ) - process_output.ovarin( - constants.MFILE, - "Scanning second variable number", - "(nsweep_2)", - scan_data.nsweep_2, - ) - process_output.ovarin( - 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} : {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.ovarin(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 + # Output optimisation parameters + solution_vector_table = [] + for i in range(nums.nvar): + nums.xcs[i] = nums.xcm[i] * nums.scafc[i] - # 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 - ) - - 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 {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.ovarin(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.ovarin( - constants.MFILE, - "Number of scan points", - "(isweep)", - scan_data.isweep, - ) - process_output.ovarin( - 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_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 1585ef1ee1..6b9d9200a4 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 @@ -1831,3 +1832,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 06a1e7ac34..af5e9c4531 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -300,10 +300,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: warn( ( "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 2d12f59148..c28bd5d77c 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 9e8703ffbe..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 balnket -
  • 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 1fd681b9a0..92d723ad36 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.numerics import PROCESSRunMode from process.models.availability import Availability from process.models.blankets.blanket_library import BlanketLibrary @@ -109,7 +110,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, @@ -119,6 +120,7 @@ CICCSuperconductingTFCoil, CROCOSuperconductingTFCoil, SuperconductingTFCoil, + SuperconductingTFTurnType, ) from process.models.vacuum import Vacuum, VacuumVessel from process.models.water_use import WaterUse @@ -432,11 +434,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) self.solver = "fsolve" @@ -446,6 +445,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(): @@ -783,6 +783,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+)