From 1a187385da5539d555a75de25054248bb8ab3392 Mon Sep 17 00:00:00 2001 From: Clair Mould <86794332+clmould@users.noreply.github.com> Date: Thu, 4 Jun 2026 13:58:45 +0100 Subject: [PATCH 1/3] convert global vars to data class --- process/core/caller.py | 14 +-- process/core/init.py | 23 ++-- process/core/input.py | 10 +- process/core/io/vary_run/config.py | 10 +- process/core/model.py | 2 + process/core/process_output.py | 18 +-- process/core/scan.py | 31 +++--- process/core/solver/evaluators.py | 3 +- process/core/solver/iteration_variables.py | 4 +- process/core/solver/solver.py | 17 +-- process/core/solver/solver_handler.py | 2 +- process/data_structure/global_variables.py | 109 ++++++++----------- process/main.py | 25 +++-- process/models/costs/costs_2015.py | 3 +- process/models/stellarator/stellarator.py | 7 +- process/models/tfcoil/base.py | 3 +- tests/integration/test_vmcon.py | 10 +- tests/unit/core/test_input.py | 22 ++-- tests/unit/core/test_mfile.py | 9 +- tests/unit/models/blankets/test_ccfe_hcpb.py | 5 +- tests/unit/models/tfcoil/test_sctfcoil.py | 3 +- tests/unit/test_main.py | 12 +- 22 files changed, 161 insertions(+), 181 deletions(-) diff --git a/process/core/caller.py b/process/core/caller.py index 384e2e3323..963d8d4af5 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -161,15 +161,13 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): for _ in range(10): # Divert OUT.DAT and MFILE.DAT output to scratch files for # idempotence checking - OutputFileManager.open_idempotence_files() + OutputFileManager.open_idempotence_files(self.data.globals.output_prefix) self._call_models_once(xc) # Write mfile finalise(self.models, self.data, ifail) # Extract data from intermediate idempotence-checking mfile - mfile_path = ( - data_structure.global_variables.output_prefix - ) + "IDEM_MFILE.DAT" + mfile_path = (self.data.globals.output_prefix) + "IDEM_MFILE.DAT" mfile = MFile(mfile_path) # Create mfile dict of float values: only compare floats mfile_data = { @@ -204,7 +202,9 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): logger.debug("Mfiles idempotent, returning") # Divert OUT.DAT and MFILE.DAT output back to original files # now idempotence checking complete - OutputFileManager.close_idempotence_files() + OutputFileManager.close_idempotence_files( + self.data.globals.output_prefix + ) # Write final output file and mfile finalise(self.models, self.data, ifail) return @@ -230,12 +230,12 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): ) # Close idempotence files, write final output file and mfile - OutputFileManager.close_idempotence_files() + OutputFileManager.close_idempotence_files(self.data.globals.output_prefix) except Exception: # If exception in model evaluations delete intermediate idempotence # files to clean up - OutputFileManager.close_idempotence_files() + OutputFileManager.close_idempotence_files(self.data.globals.output_prefix) raise else: finalise( diff --git a/process/core/init.py b/process/core/init.py index 79cbd5e575..497654e1fb 100644 --- a/process/core/init.py +++ b/process/core/init.py @@ -48,7 +48,7 @@ def init_process(data: DataStructure): iteration_variables.initialise_iteration_variables() # Creating and open the files MFile and OUTFile - process_output.OutputFileManager.open_files() + process_output.OutputFileManager.open_files(data.globals.output_prefix) # Input any desired new initial values inputs = parse_input_file(data) @@ -65,7 +65,7 @@ def init_process(data: DataStructure): # Check input data for errors/ambiguities check_process(inputs, data) - run_summary() + run_summary(data) def get_git_summary() -> tuple[str, str]: @@ -103,7 +103,7 @@ def get_git_summary() -> tuple[str, str]: return git_branch, git_tag -def run_summary(): +def run_summary(data: DataStructure): """Write a summary of the PROCESS run to the output file and MFile""" # Outfile and terminal # for outfile in [constants.NOUT, constants.IOTTY]: @@ -139,12 +139,12 @@ def run_summary(): process_output.ocmmnt(outfile, f"Computer : {machine}") process_output.ocmmnt(outfile, f"Directory : {Path.cwd()}") - fileprefix = data_structure.global_variables.fileprefix + fileprefix = data.globals.fileprefix process_output.ocmmnt( outfile, f"Input : {fileprefix}", ) - runtitle = data_structure.global_variables.runtitle + runtitle = data.globals.runtitle process_output.ocmmnt( outfile, f"Run title : {runtitle}", @@ -152,7 +152,7 @@ def run_summary(): process_output.ocmmnt( outfile, - f"Run type : Reactor concept design: {data_structure.global_variables.icase}, (c) UK Atomic Energy Authority", + f"Run type : Reactor concept design: {data.globals.icase}, (c) UK Atomic Energy Authority", ) process_output.oblnkl(outfile) @@ -177,7 +177,7 @@ def run_summary(): if data_structure.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: process_output.ocmmnt( outfile, - f"Max iterations : {data_structure.global_variables.maxcal}", + f"Max iterations : {data.globals.maxcal}", ) if data_structure.numerics.minmax > 0: @@ -233,7 +233,6 @@ def init_all_module_vars(): """ logging_model_handler.clear_logs() data_structure.numerics.init_numerics() - data_structure.global_variables.init_global_variables() init_scan_variables() constants.init_constants() @@ -561,7 +560,7 @@ def check_process(inputs, data): # noqa: ARG001 # Tight aspect ratio options (ST) if data.physics.itart == 1: - data_structure.global_variables.icase = "Tight aspect ratio tokamak model" + data.globals.icase = "Tight aspect ratio tokamak model" # Disabled Forcing that no inboard breeding blanket is used # Disabled i_blkt_inboard = 0 @@ -748,7 +747,7 @@ def check_process(inputs, data): # noqa: ARG001 # Pulsed power plant model if data.pulse.i_pulsed_plant == 1: - data_structure.global_variables.icase = "Pulsed tokamak model" + data.globals.icase = "Pulsed tokamak model" else: data.buildings.esbldgm3 = 0.0 @@ -1154,6 +1153,6 @@ def set_active_constraints(): def set_device_type(data): if data.ife.ife == 1: - data_structure.global_variables.icase = "Inertial Fusion model" + data.globals.icase = "Inertial Fusion model" elif data.stellarator.istell != 0: - data_structure.global_variables.icase = "Stellarator model" + data.globals.icase = "Stellarator model" diff --git a/process/core/input.py b/process/core/input.py index ead47bdc8e..60a7ae6afe 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -103,15 +103,15 @@ def __post_init__(self): INPUT_VARIABLES = { - "runtitle": InputVariable(data_structure.global_variables, str), - "verbose": InputVariable(data_structure.global_variables, int, choices=[0, 1]), - "run_tests": InputVariable(data_structure.global_variables, int, choices=[0, 1]), + "runtitle": InputVariable("globals", str), + "verbose": InputVariable("globals", int, choices=[0, 1]), + "run_tests": InputVariable("globals", int, choices=[0, 1]), "ioptimz": InputVariable(data_structure.numerics, int, choices=[1, -2]), "epsvmc": InputVariable(data_structure.numerics, float, range=(0.0, 1.0)), "boundl": InputVariable(data_structure.numerics, float, array=True), "boundu": InputVariable(data_structure.numerics, float, array=True), "epsfcn": InputVariable(data_structure.numerics, float, range=(0.0, 1.0)), - "maxcal": InputVariable(data_structure.global_variables, int, range=(0, 10000)), + "maxcal": InputVariable("globals", int, range=(0, 10000)), "minmax": InputVariable(data_structure.numerics, int), "neqns": InputVariable( data_structure.numerics, int, range=(0, ConstraintManager.num_constraints()) @@ -1144,7 +1144,7 @@ def __post_init__(self): def parse_input_file(data_structure_obj: DataStructure): - input_file = data_structure.global_variables.fileprefix + input_file = data_structure_obj.globals.fileprefix input_file_path = Path("IN.DAT") if input_file != "": diff --git a/process/core/io/vary_run/config.py b/process/core/io/vary_run/config.py index a289cc13b3..fef0bab0d5 100644 --- a/process/core/io/vary_run/config.py +++ b/process/core/io/vary_run/config.py @@ -13,7 +13,6 @@ import click from numpy.random import default_rng -from process import data_structure from process.core.io.in_dat import InDat from process.core.io.mfile import MFile from process.core.io.vary_run.tools import ( @@ -26,6 +25,7 @@ process_warnings, set_variable_in_indat, ) +from process.core.model import DataStructure logger = logging.getLogger(__name__) @@ -261,7 +261,7 @@ def error_status2readme(self, mfile): def modify_in_dat(self): """Modifies the original IN.DAT file""" - def setup(self): + def setup(self, data: DataStructure): """Sets up the program for running""" self.echo() @@ -275,10 +275,8 @@ def setup(self): self.generator = default_rng(seed=self.u_seed) - data_structure.global_variables.output_prefix = str( - self.wdir / self.outfile.strip("MFILE.DAT") - ) - data_structure.global_variables.fileprefix = str(self.wdir / self.infile) + data.globals.output_prefix = str(self.wdir / self.outfile.strip("MFILE.DAT")) + data.globals.fileprefix = str(self.wdir / self.infile) @staticmethod def run_process(input_path: Path, solver: str = "vmcon"): diff --git a/process/core/model.py b/process/core/model.py index dff9ccc1f1..b8acd843d5 100644 --- a/process/core/model.py +++ b/process/core/model.py @@ -14,6 +14,7 @@ from process.data_structure.divertor_variables import DivertorData from process.data_structure.first_wall_variables import FirstWallData from process.data_structure.fwbs_variables import FWBSData +from process.data_structure.global_variables import GlobalData from process.data_structure.heat_transport_variables import HeatTransportData from process.data_structure.ife_variables import IFEData from process.data_structure.impurity_radiation_variables import ImpurityRadiationData @@ -75,6 +76,7 @@ class DataStructure: rebco: RebcoData = initialise_later tfcoil: TFData = initialise_later superconducting_tfcoil: SuperconductingTFData = initialise_later + globals: GlobalData = initialise_later def __post_init__(self): for f in fields(self): diff --git a/process/core/process_output.py b/process/core/process_output.py index ff8956bd8b..fb8340c617 100644 --- a/process/core/process_output.py +++ b/process/core/process_output.py @@ -4,38 +4,38 @@ import numpy as np from process.core import constants -from process.data_structure import global_variables, numerics +from process.data_structure import numerics class OutputFileManager: @classmethod - def open_files(cls, *, mode="w"): + def open_files(cls, output_prefix: str, *, mode="w"): cls._outfile = open( # noqa: SIM115 - Path(global_variables.output_prefix + "OUT.DAT"), mode + Path(output_prefix + "OUT.DAT"), mode ) cls._mfile = open( # noqa: SIM115 - Path(global_variables.output_prefix + "MFILE.DAT"), mode + Path(output_prefix + "MFILE.DAT"), mode ) @classmethod - def open_idempotence_files(cls): + def open_idempotence_files(cls, output_prefix: str): cls._outfile.close() cls._mfile.close() cls._outfile = open( # noqa: SIM115 - Path(global_variables.output_prefix + "IDEM_OUT.DAT"), "w" + Path(output_prefix + "IDEM_OUT.DAT"), "w" ) cls._mfile = open( # noqa: SIM115 - Path(global_variables.output_prefix + "IDEM_MFILE.DAT"), "w" + Path(output_prefix + "IDEM_MFILE.DAT"), "w" ) @classmethod - def close_idempotence_files(cls): + def close_idempotence_files(cls, output_prefix: str): Path(cls._outfile.name).unlink() Path(cls._mfile.name).unlink() cls._outfile.close() cls._mfile.close() - cls.open_files(mode="a") + cls.open_files(mode="a", output_prefix=output_prefix) @classmethod def finish(cls): diff --git a/process/core/scan.py b/process/core/scan.py index 30c157676b..d4f7d3b3ad 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -16,7 +16,6 @@ from process.core.solver import constraints from process.core.solver.solver_handler import SolverHandler from process.data_structure import ( - global_variables, numerics, scan_variables, ) @@ -386,7 +385,7 @@ def post_optimise(self, ifail: int): constants.NOUT, "VMCON convergence parameter", "(convergence_parameter)", - global_variables.convergence_parameter, + self.data.globals.convergence_parameter, "OP ", ) process_output.ovarin( @@ -993,8 +992,8 @@ def scan_2d_init(): ) def scan_1d_write_point_header(self, iscan: int): - global_variables.iscan_global = iscan - global_variables.vlabel, global_variables.xlabel = self.scan_select( + self.data.globals.iscan_global = iscan + self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( scan_variables.nsweep, scan_variables.sweep, iscan ) @@ -1003,8 +1002,8 @@ def scan_1d_write_point_header(self, iscan: int): process_output.write( constants.NOUT, - f"***** Scan point {iscan} of {scan_variables.isweep} : {global_variables.xlabel}" - f", {global_variables.vlabel} = {scan_variables.sweep[iscan - 1]} " + f"***** Scan point {iscan} of {scan_variables.isweep} : {self.data.globals.xlabel}" + f", {self.data.globals.vlabel} = {scan_variables.sweep[iscan - 1]} " "*****", ) process_output.ostars(constants.NOUT, 110) @@ -1013,7 +1012,7 @@ def scan_1d_write_point_header(self, iscan: int): print( f"Starting scan point {iscan} of {scan_variables.isweep} : " - f"{global_variables.xlabel} , {global_variables.vlabel}" + f"{self.data.globals.xlabel} , {self.data.globals.vlabel}" f" = {scan_variables.sweep[iscan - 1]}" ) @@ -1021,12 +1020,12 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): iscan_r = scan_variables.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2 # Makes iscan available globally (read-only) - global_variables.iscan_global = iscan + self.data.globals.iscan_global = iscan - global_variables.vlabel, global_variables.xlabel = self.scan_select( + self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( scan_variables.nsweep, scan_variables.sweep, iscan_1 ) - global_variables.vlabel_2, global_variables.xlabel_2 = self.scan_select( + self.data.globals.vlabel_2, self.data.globals.xlabel_2 = self.scan_select( scan_variables.nsweep_2, scan_variables.sweep_2, iscan_r ) @@ -1036,8 +1035,8 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): process_output.write( constants.NOUT, f"***** 2D Scan point {iscan} of {scan_variables.isweep * scan_variables.isweep_2} : " - f"{global_variables.vlabel} = {scan_variables.sweep[iscan_1 - 1]} and" - f" {global_variables.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} " + f"{self.data.globals.vlabel} = {scan_variables.sweep[iscan_1 - 1]} and" + f" {self.data.globals.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} " "*****", ) process_output.ostars(constants.NOUT, 110) @@ -1045,10 +1044,10 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan) print( - f"Starting scan point {iscan}: {global_variables.xlabel}, " - f"{global_variables.vlabel} = {scan_variables.sweep[iscan_1 - 1]}" - f" and {global_variables.xlabel_2}, " - f"{global_variables.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} " + f"Starting scan point {iscan}: {self.data.globals.xlabel}, " + f"{self.data.globals.vlabel} = {scan_variables.sweep[iscan_1 - 1]}" + f" and {self.data.globals.xlabel_2}, " + f"{self.data.globals.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} " ) return iscan_r diff --git a/process/core/solver/evaluators.py b/process/core/solver/evaluators.py index d892bb6e61..8677e0edd6 100644 --- a/process/core/solver/evaluators.py +++ b/process/core/solver/evaluators.py @@ -5,7 +5,6 @@ from process.core.caller import Caller from process.core.model import DataStructure -from process.data_structure import global_variables as gv from process.data_structure import numerics logger = logging.getLogger(__name__) @@ -64,7 +63,7 @@ def fcnvmc1(self, _n, m, xv, ifail): objf, conf = self.caller.call_models(xv, m) # Verbose diagnostics - if gv.verbose == 1: + if self.data.globals.verbose == 1: summ = 0.0 for i in range(m): summ += conf[i] ** 2 diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 4b858c7a33..26ef603195 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -302,8 +302,8 @@ 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_structure.global_variables.vlabel, - data_structure.global_variables.vlabel_2, + data.globals.vlabel, + data.globals.vlabel_2, }: warn( ( diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index a285de39c7..167ff4b2cb 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -17,8 +17,9 @@ from scipy.optimize import fsolve from process.core.exceptions import ProcessValueError +from process.core.model import DataStructure from process.core.solver.evaluators import Evaluators -from process.data_structure import global_variables, numerics +from process.data_structure import numerics logger = logging.getLogger(__name__) @@ -26,12 +27,14 @@ class _Solver(ABC): """Base class for different solver implementations.""" - def __init__(self): + def __init__(self, *, maxcal: int | None = None): """Initialise a solver.""" # Exit code for the solver self.ifail = 0 self.tolerance = numerics.epsvmc self.b: float | None = None + self.convergence_parameter: float | None = None + self.maxcal = maxcal def set_evaluators(self, evaluators: Evaluators): """Set objective and constraint functions and their gradient evaluators. @@ -183,7 +186,7 @@ def solve(self) -> int: def _solver_callback(i: int, _result, _x, convergence_param: float): numerics.nviter = i + 1 - global_variables.convergence_parameter = convergence_param + self.convergence_parameter = convergence_param print( f"{i + 1} | Convergence Parameter: {convergence_param:.3E}", end="\r", @@ -232,7 +235,7 @@ def _ineq_cons_satisfied( np.array(self.x_0), np.array(self.bndl), np.array(self.bndu), - max_iter=global_variables.maxcal, + max_iter=self.maxcal, epsilon=self.tolerance, qsp_options={"solver": cvxpy.CLARABEL}, initial_B=bb, @@ -340,7 +343,7 @@ def solve(self) -> int: return self.info -def get_solver(solver_name: str = "vmcon") -> _Solver: +def get_solver(data: DataStructure, solver_name: str = "vmcon") -> _Solver: """Return a solver instance. Parameters @@ -356,9 +359,9 @@ def get_solver(solver_name: str = "vmcon") -> _Solver: solver: _Solver if solver_name == "vmcon": - solver = Vmcon() + solver = Vmcon(maxcal=data.globals.maxcal) elif solver_name == "vmcon_bounded": - solver = VmconBounded() + solver = VmconBounded(maxcal=data.globals.maxcal) elif solver_name == "fsolve": solver = FSolve() else: diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index 5094c5a9e2..03e76a687a 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -51,7 +51,7 @@ def run(self): evaluators = Evaluators(self.models, self.data, x) # Configure solver for problem - self.solver = get_solver(self.solver_name) + self.solver = get_solver(self.data, self.solver_name) self.solver.set_evaluators(evaluators) self.solver.set_bounds(bndl, bndu) self.solver.set_opt_params(x) diff --git a/process/data_structure/global_variables.py b/process/data_structure/global_variables.py index b599b85961..80bfc7370f 100644 --- a/process/data_structure/global_variables.py +++ b/process/data_structure/global_variables.py @@ -1,66 +1,43 @@ -icase: str = None -"""Power plant type""" - -runtitle: str -"""A short descriptive title for the run""" - -run_tests: int = None - -verbose: int = None - -maxcal: int = None -"""Maximum number of solver iterations""" - -fileprefix: str = None -"""Path to input file""" - -output_prefix: str = None -"""Output file path prefix""" - -xlabel: str = None -"""Scan parameter description label""" - -vlabel: str = None -"""Scan value name label""" - -xlabel_2: str = None -"""Scan parameter description label (2nd dimension)""" - -vlabel_2: str = None -"""Scan value name label (2nd dimension)""" - -iscan_global: int = None - -convergence_parameter: int = None -"""VMCON convergence parameter 'sum'""" - - -def init_global_variables(): - global \ - icase, \ - runtitle, \ - run_tests, \ - verbose, \ - maxcal, \ - fileprefix, \ - output_prefix, \ - xlabel, \ - vlabel, \ - xlabel_2, \ - vlabel_2, \ - iscan_global, \ - convergence_parameter - - icase = "Steady-state tokamak model" - runtitle = "Run Title (change this line using input variable 'runtitle')" - run_tests = 0 - verbose = 0 - maxcal = 200 - fileprefix = "" - output_prefix = "" - xlabel = "" - vlabel = "" - xlabel_2 = "" - vlabel_2 = "" - iscan_global = 0 - convergence_parameter = 0.0 +from dataclasses import dataclass + + +@dataclass(slots=True) +class GlobalData: + icase: str = "Steady-state tokamak model" + """Power plant type""" + + runtitle: str = "Run Title (change this line using input variable 'runtitle')" + """A short descriptive title for the run""" + + run_tests: int = 0 + + verbose: int = 0 + + maxcal: int = 200 + """Maximum number of solver iterations""" + + fileprefix: str = "" + """Path to input file""" + + output_prefix: str = "" + """Output file path prefix""" + + xlabel: str = "" + """Scan parameter 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 + + convergence_parameter: float = 0.0 + """VMCON convergence parameter 'sum'""" + + +CREATE_DICTS_FROM_DATACLASS = GlobalData diff --git a/process/main.py b/process/main.py index 6318866185..ef0c6b8528 100644 --- a/process/main.py +++ b/process/main.py @@ -267,7 +267,12 @@ class VaryRun: README.txt - contains comments from config file """ - def __init__(self, config_file: str, solver: str = "vmcon"): + def __init__( + self, + config_file: str, + solver: str = "vmcon", + data_structure: DataStructure | None = None, + ): """Initialise and perform a VaryRun. Parameters @@ -280,7 +285,7 @@ def __init__(self, config_file: str, solver: str = "vmcon"): # Store the absolute path to the config file immediately: various # dir changes happen in old run_process code self.config = RunProcessConfig.from_file(Path(config_file).resolve(), solver) - self.data = DataStructure() + self.data = data_structure or DataStructure() @property def mfile_path(self): @@ -295,7 +300,7 @@ def run(self): if input file doesn't exist """ init.init_all_module_vars() - self.config.setup() + self.config.setup(self.data) setup_loggers(Path(self.config.wdir) / "process.log") @@ -328,11 +333,11 @@ def __init__( which solver to use, as specified in solver.py """ self.input_file = Path(input_file) + self.data = data_structure or DataStructure() self.validate_input(update_obsolete) self.init_module_vars() self.set_filenames(filepath_out) - self.data = data_structure or DataStructure() self.initialise() self.models = Models(self.data) self.solver = solver @@ -368,7 +373,7 @@ def set_filenames(self, filepath_out): .replace("MFILE.DAT", "") ).strip() self.set_input() - data_structure.global_variables.output_prefix = ( + self.data.globals.output_prefix = ( f"{Path(self.filepath).as_posix().strip()}/" if not self.filename_prefix else Path(self.filepath, self.filename_prefix).as_posix().strip() @@ -397,22 +402,18 @@ def set_input(self): ) # Set the input file in the Fortran - data_structure.global_variables.fileprefix = self.input_path.resolve() + self.data.globals.fileprefix = self.input_path.resolve() def set_output(self): """Set the output file name. Set Path object on the Process object, and set the prefix in the Fortran. """ - self.output_path = Path( - data_structure.global_variables.output_prefix + "OUT.DAT" - ) + self.output_path = Path(self.data.globals.output_prefix + "OUT.DAT") def set_mfile(self): """Set the mfile filename.""" - self.mfile_path = Path( - data_structure.global_variables.output_prefix + "MFILE.DAT" - ) + self.mfile_path = Path(self.data.globals.output_prefix + "MFILE.DAT") def initialise(self): """Run the init module to call all initialisation routines.""" diff --git a/process/models/costs/costs_2015.py b/process/models/costs/costs_2015.py index 579270d23d..e62a0e28ee 100644 --- a/process/models/costs/costs_2015.py +++ b/process/models/costs/costs_2015.py @@ -5,7 +5,6 @@ from process.core import constants from process.core import process_output as po from process.core.model import Model -from process.data_structure import global_variables logger = logging.getLogger(__name__) @@ -147,7 +146,7 @@ def calc_fwbs_costs(self): tail_li6 = feed_li6 * 0.75e0 # Built-in test - if global_variables.run_tests == 1: + if self.data.globals.run_tests == 1: product_li6 = 0.3 feed_to_product_mass_ratio = (product_li6 - tail_li6) / (feed_li6 - tail_li6) tail_to_product_mass_ratio = (product_li6 - feed_li6) / (feed_li6 - tail_li6) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index a15a6285d7..7478983af9 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -13,10 +13,7 @@ from process.core.coolprop_interface import FluidProperties from process.core.exceptions import ProcessValueError from process.core.model import Model -from process.data_structure import ( - global_variables, - numerics, -) +from process.data_structure import numerics from process.models.engineering.pumping import CoolantType from process.models.physics.physics import Physics, rether from process.models.power import PumpingPowerModelTypes @@ -208,7 +205,7 @@ def st_new_config(self): """ load_stellarator_config( self.data.stellarator.istell, - Path(f"{global_variables.output_prefix}stella_conf.json"), + Path(f"{self.data.globals.output_prefix}stella_conf.json"), self.data, ) diff --git a/process/models/tfcoil/base.py b/process/models/tfcoil/base.py index c3ace698bd..dc7b0c8878 100644 --- a/process/models/tfcoil/base.py +++ b/process/models/tfcoil/base.py @@ -16,7 +16,6 @@ from process.core import process_output as po from process.core.exceptions import ProcessValueError from process.core.model import DataStructure, Model -from process.data_structure import global_variables from process.data_structure.physics_variables import DivertorNumberModels logger = logging.getLogger(__name__) @@ -3609,7 +3608,7 @@ def table_format_arrays(a, mult=1, delim="\t\t"): for k, v in sig_file_data.items() } - sig_filename = global_variables.output_prefix + "SIG_TF.json" + sig_filename = self.data.globals.output_prefix + "SIG_TF.json" with open(sig_filename, "w") as f: json.dump(sig_file_data, f) diff --git a/tests/integration/test_vmcon.py b/tests/integration/test_vmcon.py index 1b179dc338..112644e3e0 100644 --- a/tests/integration/test_vmcon.py +++ b/tests/integration/test_vmcon.py @@ -13,6 +13,7 @@ import pytest from process.core.init import init_all_module_vars +from process.core.model import DataStructure from process.core.solver.evaluators import Evaluators from process.core.solver.solver import get_solver @@ -26,6 +27,11 @@ def reinit(): init_all_module_vars() +@pytest.fixture +def data_structure_obj(): + return DataStructure() + + class Case: """A Vmcon test case.""" @@ -610,7 +616,7 @@ def case(request): return case_fn() -def test_vmcon(case, solver_name): +def test_vmcon(case, solver_name, data_structure_obj): """Integration test for Vmcon. :param case: a Vmcon scenario and its expected result @@ -626,7 +632,7 @@ def test_vmcon(case, solver_name): logger.debug(f"x[{i}] = {case.solver_args.x[i]}") # Configure solver for problem - solver = get_solver(solver_name) + solver = get_solver(data_structure_obj, solver_name) solver.set_evaluators(case.evaluator) solver.set_opt_params(case.solver_args.x) solver.set_bounds( diff --git a/tests/unit/core/test_input.py b/tests/unit/core/test_input.py index ffcba21edb..d11a15415b 100644 --- a/tests/unit/core/test_input.py +++ b/tests/unit/core/test_input.py @@ -66,7 +66,7 @@ def test_parse_real(epsvmc, expected, tmp_path, data_structure_obj): Program to get the expected value for 0.008 provided at https://github.com/ukaea/PROCESS/pull/3067 """ - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, f"epsvmc = {epsvmc}" ) init.init_process(data_structure_obj) @@ -91,7 +91,7 @@ def test_exact_parsing(value, tmp_path, data_structure_obj): These tests failed using the old input parser and serve to show that the Python parser generally produces more accurate floats and accumulates less error. """ - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, f"epsvmc = {value}" ) init.init_process(data_structure_obj) @@ -100,20 +100,20 @@ def test_exact_parsing(value, tmp_path, data_structure_obj): def test_parse_input(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("runtitle = my run title\nioptimz = -2\nepsvmc = 0.6\nboundl(1) = 0.5"), ) init.init_process(data_structure_obj) - assert data_structure.global_variables.runtitle == "my run title" + assert data_structure_obj.globals.runtitle == "my run title" assert data_structure.numerics.ioptimz == PROCESSRunMode.EVALUATION assert pytest.approx(data_structure.numerics.epsvmc) == 0.6 assert pytest.approx(data_structure.numerics.boundl[0]) == 0.5 def test_input_choices(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("ioptimz = -1") ) @@ -125,7 +125,7 @@ def test_input_choices(tmp_path, data_structure_obj): ("input_file_value"), [(-0.01,), (1.1,)], ids=("violate lower", "violate upper") ) def test_input_range(tmp_path, input_file_value, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, (f"epsvmc = {input_file_value}") ) @@ -137,7 +137,7 @@ def test_input_range(tmp_path, input_file_value, data_structure_obj): def test_input_array_when_not(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("epsvmc(1) = 0.5") ) @@ -146,7 +146,7 @@ def test_input_array_when_not(tmp_path, data_structure_obj): def test_input_not_array_when_is(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("boundl = 0.5") ) @@ -155,7 +155,7 @@ def test_input_not_array_when_is(tmp_path, data_structure_obj): def test_input_float_when_int(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("ioptimz = 0.5") ) @@ -164,7 +164,7 @@ def test_input_float_when_int(tmp_path, data_structure_obj): def test_input_array(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("boundl = 0.1, 0.2, 1.0, 0.0, 1.0e2") ) @@ -175,7 +175,7 @@ def test_input_array(tmp_path, data_structure_obj): def test_input_on_new_data_structure(tmp_path, data_structure_obj): - data_structure.global_variables.fileprefix = _create_input_file( + data_structure_obj.globals.fileprefix = _create_input_file( tmp_path, ("windspeed = 1.22") ) diff --git a/tests/unit/core/test_mfile.py b/tests/unit/core/test_mfile.py index 8426e187a8..ed5961a7ac 100644 --- a/tests/unit/core/test_mfile.py +++ b/tests/unit/core/test_mfile.py @@ -9,12 +9,17 @@ from process.core.model import DataStructure -def test_get_mfile_initial_ixc_values(input_file, tmp_path): +@pytest.fixture +def data_structure_obj(): + return DataStructure() + + +def test_get_mfile_initial_ixc_values(input_file, tmp_path, data_structure_obj): tmp_input_file = tmp_path / "IN.DAT" shutil.copy(input_file, tmp_input_file) iteration_variable_names, iteration_variable_values = get_mfile_initial_ixc_values( - Path(tmp_input_file), DataStructure() + Path(tmp_input_file), data_structure_obj ) assert iteration_variable_names[0] == "b_plasma_toroidal_on_axis" diff --git a/tests/unit/models/blankets/test_ccfe_hcpb.py b/tests/unit/models/blankets/test_ccfe_hcpb.py index 24e34d32be..dd06cf00e2 100644 --- a/tests/unit/models/blankets/test_ccfe_hcpb.py +++ b/tests/unit/models/blankets/test_ccfe_hcpb.py @@ -2,7 +2,6 @@ import pytest -from process.data_structure import global_variables from process.data_structure.build_variables import InboardBlanketConfiguration from process.models.engineering.pumping import CoolantType @@ -339,7 +338,9 @@ def test_nuclear_heating_magnets(nuclearheatingmagnetsparam, monkeypatch, ccfe_h ccfe_hcpb.data.tfcoil, "whttflgs", nuclearheatingmagnetsparam.whttflgs ) - monkeypatch.setattr(global_variables, "verbose", nuclearheatingmagnetsparam.verbose) + monkeypatch.setattr( + ccfe_hcpb.data.globals, "verbose", nuclearheatingmagnetsparam.verbose + ) monkeypatch.setattr( ccfe_hcpb.data.ccfe_hcpb, diff --git a/tests/unit/models/tfcoil/test_sctfcoil.py b/tests/unit/models/tfcoil/test_sctfcoil.py index 75aedf4398..1f4ed0339c 100644 --- a/tests/unit/models/tfcoil/test_sctfcoil.py +++ b/tests/unit/models/tfcoil/test_sctfcoil.py @@ -3,7 +3,6 @@ import numpy as np import pytest -from process.data_structure import global_variables from process.models.tfcoil import superconducting as sctf @@ -429,7 +428,7 @@ def test_supercon(superconparam, monkeypatch, cicc_sctfcoil): superconparam.f_b_tf_inboard_peak_ripple_symmetric, ) - monkeypatch.setattr(global_variables, "run_tests", superconparam.run_tests) + monkeypatch.setattr(cicc_sctfcoil.data.globals, "run_tests", superconparam.run_tests) tf_limits = cicc_sctfcoil.tf_cable_in_conduit_superconductor_properties( i_tf_superconductor=superconparam.i_tf_superconductor, diff --git a/tests/unit/test_main.py b/tests/unit/test_main.py index a4bab469a4..9529eefaed 100644 --- a/tests/unit/test_main.py +++ b/tests/unit/test_main.py @@ -5,7 +5,6 @@ import pytest -from process import data_structure from process.core.model import DataStructure from process.main import SingleRun, VaryRun @@ -71,7 +70,7 @@ def test_set_input(single_run, monkeypatch, input_file): # Mocks set up, can now run set_input() single_run.set_input() # Check path has been set - assert data_structure.global_variables.fileprefix == expected + assert single_run.data.globals.fileprefix == expected def test_set_output(single_run, monkeypatch): @@ -85,14 +84,11 @@ def test_set_output(single_run, monkeypatch): # Expected output prefix expected = "output_prefix" # Mock self.filename_prefix on single_run with the value of expected - monkeypatch.setattr(data_structure.global_variables, "output_prefix", expected) + monkeypatch.setattr(single_run.data.globals, "output_prefix", expected) - # Mocking undo trys to set the value as none - # monkeypatch.setattr(data_structure.global_variables, "output_prefix", None) - # Run the method, and extract the value single_run.set_output() - assert Path(data_structure.global_variables.output_prefix).name == expected + assert Path(single_run.data.globals.output_prefix).name == expected def test_initialise(single_run, monkeypatch): @@ -117,7 +113,7 @@ def test_set_mfile(single_run, monkeypatch): prefix = "test" expected = Path(prefix + "MFILE.DAT") # Mock filename_prefix and run - monkeypatch.setattr(data_structure.global_variables, "output_prefix", prefix) + monkeypatch.setattr(single_run.data.globals, "output_prefix", prefix) single_run.set_mfile() assert single_run.mfile_path.name == expected.name From 86cb81f6ea84d2a0794406246c2ad50cbf6d91e2 Mon Sep 17 00:00:00 2001 From: Clair Mould <86794332+clmould@users.noreply.github.com> Date: Tue, 9 Jun 2026 13:22:55 +0100 Subject: [PATCH 2/3] Convert scan variables to dataclass --- process/core/init.py | 2 - process/core/input.py | 32 +-- process/core/model.py | 2 + process/core/scan.py | 102 +++++---- process/data_structure/scan_variables.py | 267 +++++++++++------------ 5 files changed, 191 insertions(+), 214 deletions(-) diff --git a/process/core/init.py b/process/core/init.py index 497654e1fb..7724f6f681 100644 --- a/process/core/init.py +++ b/process/core/init.py @@ -21,7 +21,6 @@ from process.data_structure.impurity_radiation_variables import N_IMPURITIES from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode from process.data_structure.physics_variables import DivertorNumberModels -from process.data_structure.scan_variables import init_scan_variables from process.models.pfcoil import PFLocationTypes from process.models.physics.profiles import DensityProfilePedestalType from process.models.stellarator.initialization import st_init @@ -233,7 +232,6 @@ def init_all_module_vars(): """ logging_model_handler.clear_logs() data_structure.numerics.init_numerics() - init_scan_variables() constants.init_constants() diff --git a/process/core/input.py b/process/core/input.py index 60a7ae6afe..97b524cfda 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -19,6 +19,7 @@ from process.data_structure.impurity_radiation_variables import N_IMPURITIES from process.data_structure.pfcoil_variables import N_PF_GROUPS_MAX from process.data_structure.physics_variables import N_CONFINEMENT_SCALINGS +from process.data_structure.scan_variables import IPNSCNS, IPNSCNV if TYPE_CHECKING: from collections.abc import Callable @@ -68,11 +69,12 @@ class InputVariable: been cast to the specified `type`. """ additional_actions: ( - Callable[[str, ValidInputTypes, int | None, InputVariable], None] | None + Callable[[str, ValidInputTypes, int | None, InputVariable, DataStructure], None] + | None ) = None - """A function that takes the input variable: name, value, array index, and config (this dataclass) - as input and performs some additional action in addition to the default actions prescribed by the variables - config. May raise a ProcessValidationError. + """A function that takes the input variable: name, value, array index, config (this dataclass), and + the data structure object as input and performs some additional action in addition to the default + actions prescribed by the variables config. May raise a ProcessValidationError. NOTE: The value passed in has already been cleaned in the default way and has been cast to the specified `type`. @@ -1049,7 +1051,7 @@ def __post_init__(self): "output_costs": InputVariable("costs", int, choices=[0, 1]), "i_p_coolant_pumping": InputVariable("fwbs", int, range=(0, 3)), "reinke_mode": InputVariable("reinke", int, choices=[0, 1]), - "scan_dim": InputVariable(data_structure.scan_variables, int, range=(1, 2)), + "scan_dim": InputVariable("scan", int, range=(1, 2)), "i_thermal_electric_conversion": InputVariable("fwbs", int, range=(0, 4)), "secondary_cycle_liq": InputVariable("fwbs", int, range=(2, 4)), "supercond_cost_model": InputVariable("costs", int, choices=[0, 1]), @@ -1091,27 +1093,27 @@ def __post_init__(self): "v2matf": InputVariable("ife", float, array=True), "v3matf": InputVariable("ife", float, array=True), "isweep": InputVariable( - data_structure.scan_variables, + "scan", int, - choices=range(data_structure.scan_variables.IPNSCNS + 1), + choices=range(IPNSCNS + 1), ), "nsweep": InputVariable( - data_structure.scan_variables, + "scan", int, - choices=range(1, data_structure.scan_variables.IPNSCNV + 1), + choices=range(1, IPNSCNV + 1), ), "isweep_2": InputVariable( - data_structure.scan_variables, + "scan", int, - choices=range(data_structure.scan_variables.IPNSCNS + 1), + choices=range(IPNSCNS + 1), ), "nsweep_2": InputVariable( - data_structure.scan_variables, + "scan", int, - choices=range(1, data_structure.scan_variables.IPNSCNV + 1), + choices=range(1, IPNSCNV + 1), ), - "sweep": InputVariable(data_structure.scan_variables, float, array=True), - "sweep_2": InputVariable(data_structure.scan_variables, float, array=True), + "sweep": InputVariable("scan", float, array=True), + "sweep_2": InputVariable("scan", float, array=True), "impvardiv": InputVariable( "reinke", int, diff --git a/process/core/model.py b/process/core/model.py index b8acd843d5..e4fa92fb05 100644 --- a/process/core/model.py +++ b/process/core/model.py @@ -27,6 +27,7 @@ from process.data_structure.pulse_variables import PulseData from process.data_structure.rebco_variables import RebcoData from process.data_structure.reinke_variables import ReinkeData +from process.data_structure.scan_variables import ScanData from process.data_structure.stellarator_configuration import StellaratorConfigData from process.data_structure.stellarator_variables import StellaratorData from process.data_structure.structure_variables import StructureData @@ -77,6 +78,7 @@ class DataStructure: tfcoil: TFData = initialise_later superconducting_tfcoil: SuperconductingTFData = initialise_later globals: GlobalData = initialise_later + scan: ScanData = initialise_later def __post_init__(self): for f in fields(self): diff --git a/process/core/scan.py b/process/core/scan.py index d4f7d3b3ad..51a6c11693 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -15,11 +15,9 @@ 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 import ( - numerics, - scan_variables, -) +from process.data_structure import numerics 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 @@ -216,7 +214,7 @@ def run_scan(self): number of output variable values are written to the MFILE.DAT file at each scan point, for plotting or other post-processing purposes. """ - if scan_variables.isweep == 0: + 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() @@ -230,14 +228,14 @@ def run_scan(self): show_errors(constants.NOUT) return - if scan_variables.isweep > scan_variables.IPNSCNS: + if self.data.scan.isweep > IPNSCNS: raise ProcessValueError( "Illegal value of isweep", - isweep=scan_variables.isweep, - IPNSCNS=scan_variables.IPNSCNS, + isweep=self.data.scan.isweep, + IPNSCNS=IPNSCNS, ) - if scan_variables.scan_dim == 2: + if self.data.scan.scan_dim == 2: self.scan_2d() else: self.scan_1d() @@ -825,7 +823,7 @@ def scan_1d(self): # initialise dict which will contain ifail values for each scan point scan_1d_ifail_dict = {} - for iscan in range(1, scan_variables.isweep + 1): + for iscan in range(1, self.data.scan.isweep + 1): self.scan_1d_write_point_header(iscan) start_time = time.time() ifail = self.doopt() @@ -841,13 +839,13 @@ def scan_1d(self): logging_model_handler.clear_logs() # outvar now contains results - self.scan_1d_write_plot() + self.scan_1d_write_plot(self.data.scan) print( " ****************************************** Scan Convergence Summary ****************************************** \n" ) - sweep_values = scan_variables.sweep[: scan_variables.isweep] + sweep_values = self.data.scan.sweep[: self.data.scan.isweep] nsweep_var_name, _ = self.scan_select( - scan_variables.nsweep, scan_variables.sweep, scan_variables.isweep + self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep ) converged_count = 0 # offsets for aligning the converged/unconverged column @@ -856,7 +854,7 @@ def scan_1d(self): max_sweep_value_length - len(str(sweep_val).replace(".", "")) for sweep_val in sweep_values ] - for iscan in range(1, scan_variables.isweep + 1): + for iscan in range(1, self.data.scan.isweep + 1): if scan_1d_ifail_dict[iscan] == 1: converged_count += 1 print( @@ -870,23 +868,23 @@ def scan_1d(self): + " " * offsets[iscan - 1] + "\u001b[31mUNCONVERGED \u001b[0m" ) - converged_percentage = converged_count / scan_variables.isweep * 100 + 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.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( - (scan_variables.NOUTVARS, scan_variables.IPNSCNS), + (NOUTVARS, IPNSCNS), dtype=np.float64, order="F", ) - for iscan_1 in range(1, scan_variables.isweep + 1): - for iscan_2 in range(1, scan_variables.isweep_2 + 1): + 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() @@ -905,13 +903,13 @@ def scan_2d(self): print( " ****************************************** Scan Convergence Summary ****************************************** \n" ) - sweep_1_values = scan_variables.sweep[: scan_variables.isweep] - sweep_2_values = scan_variables.sweep_2[: scan_variables.isweep_2] + 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( - scan_variables.nsweep, scan_variables.sweep, scan_variables.isweep + self.data.scan.nsweep, self.data.scan.sweep, self.data.scan.isweep ) nsweep_2_var_name, _ = self.scan_select( - scan_variables.nsweep_2, scan_variables.sweep_2, scan_variables.isweep_2 + self.data.scan.nsweep_2, self.data.scan.sweep_2, self.data.scan.isweep_2 ) converged_count = 0 scan_point = 1 @@ -919,7 +917,7 @@ def scan_2d(self): 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( - (scan_variables.isweep, scan_variables.isweep_2), dtype=int, order="F" + (self.data.scan.isweep, self.data.scan.isweep_2), dtype=int, order="F" ) for count1, sweep1 in enumerate(sweep_1_values): for count2, sweep2 in enumerate(sweep_2_values): @@ -930,8 +928,8 @@ def scan_2d(self): - len(str(sweep2).replace(".", "")) ) - for iscan_1 in range(1, scan_variables.isweep + 1): - for iscan_2 in range(1, scan_variables.isweep_2 + 1): + 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( @@ -948,53 +946,53 @@ def scan_2d(self): ) scan_point += 1 converged_percentage = ( - converged_count / (scan_variables.isweep * scan_variables.isweep_2) * 100 + converged_count / (self.data.scan.isweep * self.data.scan.isweep_2) * 100 ) print(f"\nConvergence Percentage: {converged_percentage:.2f}%") @staticmethod - def scan_2d_init(): + def scan_2d_init(scan_data: ScanData): process_output.ovarin( constants.MFILE, "Number of first variable scan points", "(isweep)", - scan_variables.isweep, + scan_data.isweep, ) process_output.ovarin( constants.MFILE, "Number of second variable scan points", "(isweep_2)", - scan_variables.isweep_2, + scan_data.isweep_2, ) process_output.ovarin( constants.MFILE, "Scanning first variable number", "(nsweep)", - scan_variables.nsweep, + scan_data.nsweep, ) process_output.ovarin( constants.MFILE, "Scanning second variable number", "(nsweep_2)", - scan_variables.nsweep_2, + scan_data.nsweep_2, ) process_output.ovarin( constants.MFILE, "Scanning second variable number", "(nsweep_2)", - scan_variables.nsweep_2, + scan_data.nsweep_2, ) process_output.ovarin( constants.MFILE, "Scanning second variable number", "(nsweep_2)", - scan_variables.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( - scan_variables.nsweep, scan_variables.sweep, iscan + self.data.scan.nsweep, self.data.scan.sweep, iscan ) process_output.oblnkl(constants.NOUT) @@ -1002,8 +1000,8 @@ def scan_1d_write_point_header(self, iscan: int): process_output.write( constants.NOUT, - f"***** Scan point {iscan} of {scan_variables.isweep} : {self.data.globals.xlabel}" - f", {self.data.globals.vlabel} = {scan_variables.sweep[iscan - 1]} " + 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) @@ -1011,22 +1009,22 @@ def scan_1d_write_point_header(self, iscan: int): process_output.ovarin(constants.MFILE, "Scan point number", "(iscan)", iscan) print( - f"Starting scan point {iscan} of {scan_variables.isweep} : " + f"Starting scan point {iscan} of {self.data.scan.isweep} : " f"{self.data.globals.xlabel} , {self.data.globals.vlabel}" - f" = {scan_variables.sweep[iscan - 1]}" + f" = {self.data.scan.sweep[iscan - 1]}" ) def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): - iscan_r = scan_variables.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2 + iscan_r = self.data.scan.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2 # Makes iscan available globally (read-only) self.data.globals.iscan_global = iscan self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( - scan_variables.nsweep, scan_variables.sweep, iscan_1 + self.data.scan.nsweep, self.data.scan.sweep, iscan_1 ) self.data.globals.vlabel_2, self.data.globals.xlabel_2 = self.scan_select( - scan_variables.nsweep_2, scan_variables.sweep_2, iscan_r + self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r ) process_output.oblnkl(constants.NOUT) @@ -1034,9 +1032,9 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): process_output.write( constants.NOUT, - f"***** 2D Scan point {iscan} of {scan_variables.isweep * scan_variables.isweep_2} : " - f"{self.data.globals.vlabel} = {scan_variables.sweep[iscan_1 - 1]} and" - f" {self.data.globals.vlabel_2} = {scan_variables.sweep_2[iscan_r - 1]} " + 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) @@ -1045,30 +1043,30 @@ def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): print( f"Starting scan point {iscan}: {self.data.globals.xlabel}, " - f"{self.data.globals.vlabel} = {scan_variables.sweep[iscan_1 - 1]}" + 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} = {scan_variables.sweep_2[iscan_r - 1]} " + f"{self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} " ) return iscan_r @staticmethod - def scan_1d_write_plot(): - if scan_variables.first_call_1d: + 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_variables.isweep, + scan_data.isweep, ) process_output.ovarin( constants.MFILE, "Scanning variable number", "(nsweep)", - scan_variables.nsweep, + scan_data.nsweep, ) - scan_variables.first_call_1d = False + scan_data.first_call_1d = False def scan_select(self, nwp, swp, iscn): match nwp: diff --git a/process/data_structure/scan_variables.py b/process/data_structure/scan_variables.py index f4cdec5d5a..9e8703ffbe 100644 --- a/process/data_structure/scan_variables.py +++ b/process/data_structure/scan_variables.py @@ -5,156 +5,133 @@ over a range of values of a particular scanning variable. """ +from dataclasses import dataclass, field + import numpy as np -IPNSCNS: int = 1000 +IPNSCNS = 1000 """Maximum number of scan points""" -IPNSCNV: int = 81 +IPNSCNV = 81 """Number of available scan variables""" -NOUTVARS: int = 84 - - -WIDTH: int = 110 - - -scan_dim: int = None -"""1-D or 2-D scan switch (1=1D, 2=2D)""" - - -isweep: int = None -"""Number of scan points to calculate""" - - -isweep_2: int = None -"""Number of 2D scan points to calculate""" - - -nsweep: int = None -"""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 = None -"""nsweep_2 /3/ : switch denoting quantity to scan for 2D scan:""" - - -sweep: list[float] = None -"""sweep(IPNSCNS) /../: actual values to use in scan""" - - -sweep_2: list[float] = None -"""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 = None - -first_call_2d: bool = None - - -def init_scan_variables(): - """Initialise the scan module""" - global \ - scan_dim, \ - isweep, \ - isweep_2, \ - nsweep, \ - nsweep_2, \ - sweep, \ - sweep_2, \ - first_call_1d, \ - first_call_2d - - scan_dim = 1 - isweep = 0 - isweep_2 = 0 - nsweep = 1 - nsweep_2 = 3 - sweep = np.zeros(1000, dtype=np.float64) - sweep_2 = np.zeros(1000, dtype=np.float64) - first_call_1d = True - first_call_2d = True +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 + """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""" + + 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 + + first_call_2d: bool = True + + +CREATE_DICTS_FROM_DATACLASS = ScanData From f0b9e2aa479ad7984ebf944ea479cad3eb61acdd Mon Sep 17 00:00:00 2001 From: Clair Mould <86794332+clmould@users.noreply.github.com> Date: Thu, 11 Jun 2026 14:37:35 +0100 Subject: [PATCH 3/3] convert numerics vars to dataclass --- documentation/source/development/add-vars.md | 4 +- process/core/caller.py | 9 +- process/core/final.py | 26 +- process/core/init.py | 213 ++-- process/core/input.py | 47 +- process/core/io/mfile/base.py | 5 +- process/core/io/plot/solutions.py | 2 - process/core/io/vary_run/config.py | 11 +- process/core/io/vary_run/tools.py | 14 +- process/core/model.py | 2 + process/core/process_output.py | 11 +- process/core/scan.py | 153 ++- process/core/solver/constraints.py | 3 +- process/core/solver/evaluators.py | 9 +- process/core/solver/iteration_variables.py | 47 +- process/core/solver/solver.py | 30 +- process/core/solver/solver_handler.py | 39 +- process/data_structure/numerics.py | 1110 ++++++++---------- process/main.py | 11 +- process/models/build.py | 5 +- process/models/physics/l_h_transition.py | 5 +- process/models/physics/physics.py | 9 +- process/models/power.py | 5 +- process/models/pulse.py | 3 +- process/models/stellarator/initialization.py | 3 +- process/models/stellarator/stellarator.py | 7 +- process/models/tfcoil/superconducting.py | 3 +- tests/unit/core/test_input.py | 13 +- tests/unit/models/test_power.py | 8 +- tests/unit/models/test_pulse.py | 6 +- 30 files changed, 851 insertions(+), 962 deletions(-) diff --git a/documentation/source/development/add-vars.md b/documentation/source/development/add-vars.md index 81ff1b49a1..f6cbd1a625 100644 --- a/documentation/source/development/add-vars.md +++ b/documentation/source/development/add-vars.md @@ -55,7 +55,7 @@ Code example in the `input.py` file: To add a `PROCESS` iteration variable please follow the steps below, in addition to the instructions for adding an input variable: -1. The parameter `ipnvars` in module `numerics` of `numerics.f90` will normally be greater than the actual number of iteration variables, and does not need to be changed. +1. The parameter `IPNVARS` in module `numerics` of `numerics.f90` will normally be greater than the actual number of iteration variables, and does not need to be changed. 2. Append a new iteration number key to the end of the `ITERATION_VARIABLES` dictionary in `iteration_variables.py`. The associated variable is the corresponding key value. 3. Set the variable origin file and then the associated lower and upper bounds 4. Update the `lablxc` description in `numerics.f90`. @@ -80,7 +80,7 @@ ITERATION_VARIABLES = { New figures of merit are added to `PROCESS` in the following way: -1. Increment the parameter `ipnfoms` in module `numerics` in source file `numerics.py` to accommodate the new figure of merit. +1. Increment the parameter `IPNFOMS` in module `numerics` in source file `numerics.py` to accommodate the new figure of merit. 2. Assign the new integer value and description string of the new figure of merit to the `FiguresOfMerit` enumerator in `numerics.py`. diff --git a/process/core/caller.py b/process/core/caller.py index 963d8d4af5..041242c69b 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -7,7 +7,6 @@ import numpy as np from tabulate import tabulate -from process import data_structure from process.core import constants from process.core.final import finalise from process.core.io.mfile import MFile @@ -99,7 +98,7 @@ def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: for _ in range(10): self._call_models_once(xc) # Evaluate objective function and constraints - objf = objective_function(data_structure.numerics.minmax, self.data) + objf = objective_function(self.data.numerics.minmax, self.data) conf, _, _, _, _ = constraints.constraint_eqns(m, -1, self.data) if objf_prev is None and conf_prev is None: @@ -261,7 +260,7 @@ def _call_models_once(self, xc: np.ndarray): nvars = len(xc) # Increment the call counter - data_structure.numerics.ncalls += 1 + self.data.numerics.ncalls += 1 # Convert variables set_scaled_iteration_variable(xc, nvars, self.data) @@ -409,8 +408,8 @@ def write_output_files( ifail : int solver return code """ - n = data_structure.numerics.nvar - x = data_structure.numerics.xcm[:n] + n = data.numerics.nvar + x = data.numerics.xcm[:n] # Call models, ensuring output mfiles are fully idempotent caller = Caller(models, data) if runtime is not None: diff --git a/process/core/final.py b/process/core/final.py index a7955cf6c5..d5da7593bc 100644 --- a/process/core/final.py +++ b/process/core/final.py @@ -7,7 +7,6 @@ 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 import numerics from process.data_structure.numerics import PROCESSRunMode @@ -33,7 +32,7 @@ def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None): po.oheadr(constants.NOUT, "Final UNFEASIBLE Point") # Output relevant to no optimisation - if numerics.ioptimz == PROCESSRunMode.EVALUATION: + if data.numerics.ioptimz == PROCESSRunMode.EVALUATION: output_evaluation(data) # Print non-idempotence warning to OUT.DAT only @@ -58,18 +57,21 @@ def output_evaluation(data): po.oblnkl(constants.NOUT) # Evaluate objective function - norm_objf = objective_function(numerics.minmax, data) + 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( - numerics.neqns + numerics.nineqns, -1, data + data.numerics.neqns + data.numerics.nineqns, -1, data ) labels = [ - numerics.lablcc[j] - for j in [i - 1 for i in numerics.icc[: numerics.neqns + numerics.nineqns]] + 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)] @@ -84,8 +86,8 @@ def output_evaluation(data): po.write(constants.NOUT, tabulate(table_data, headers="keys")) - for i in range(numerics.neqns): - constraint_id = numerics.icc[i] + for i in range(data.numerics.neqns): + constraint_id = data.numerics.icc[i] po.ovarre( constants.MFILE, f"{labels[i]} normalised residue", @@ -93,11 +95,11 @@ def output_evaluation(data): residual_error[i], ) - for i in range(numerics.nineqns): - constraint_id = numerics.icc[numerics.neqns + i] + for i in range(data.numerics.nineqns): + constraint_id = data.numerics.icc[data.numerics.neqns + i] po.ovarre( constants.MFILE, - f"{labels[numerics.neqns + i]}", + f"{labels[data.numerics.neqns + i]}", f"(ineq_con{constraint_id:03d})", - residual_error[numerics.neqns + i], + residual_error[data.numerics.neqns + i], ) diff --git a/process/core/init.py b/process/core/init.py index 7724f6f681..c4bd99ea7e 100644 --- a/process/core/init.py +++ b/process/core/init.py @@ -9,7 +9,6 @@ from warnings import warn import process -from process import data_structure from process.core import constants, process_output from process.core.exceptions import ProcessValidationError from process.core.input import parse_input_file @@ -44,7 +43,7 @@ def init_process(data: DataStructure): the input file, and checks the run parameters for consistency. """ # Initialise the program variables - iteration_variables.initialise_iteration_variables() + iteration_variables.initialise_iteration_variables(data) # Creating and open the files MFile and OUTFile process_output.OutputFileManager.open_files(data.globals.output_prefix) @@ -53,7 +52,7 @@ def init_process(data: DataStructure): inputs = parse_input_file(data) # Set active constraints - set_active_constraints() + set_active_constraints(data) # set the device type (icase) set_device_type(data) @@ -158,42 +157,38 @@ def run_summary(data: DataStructure): process_output.ostars(outfile, 110) process_output.oblnkl(outfile) - process_output.ocmmnt( - outfile, f"Equality constraints : {data_structure.numerics.neqns}" - ) + process_output.ocmmnt(outfile, f"Equality constraints : {data.numerics.neqns}") process_output.ocmmnt( outfile, - f"Inequality constraints : {data_structure.numerics.nineqns}", + f"Inequality constraints : {data.numerics.nineqns}", ) process_output.ocmmnt( outfile, - f"Total constraints : {data_structure.numerics.nineqns + data_structure.numerics.neqns}", - ) - process_output.ocmmnt( - outfile, f"Iteration variables : {data_structure.numerics.nvar}" + f"Total constraints : {data.numerics.nineqns + data.numerics.neqns}", ) + process_output.ocmmnt(outfile, f"Iteration variables : {data.numerics.nvar}") # If optimising, write objective function and convergence parameter - if data_structure.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: + if data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: process_output.ocmmnt( outfile, f"Max iterations : {data.globals.maxcal}", ) - if data_structure.numerics.minmax > 0: + if data.numerics.minmax > 0: minmax_string = " -- minimise " minmax_sign = "+" else: minmax_string = " -- maximise " minmax_sign = "-" - fom_string = FiguresOfMerit(abs(data_structure.numerics.minmax)).description + fom_string = FiguresOfMerit(abs(data.numerics.minmax)).description process_output.ocmmnt( outfile, - f"Figure of merit : {minmax_sign}{abs(data_structure.numerics.minmax)}{minmax_string}{fom_string}", + f"Figure of merit : {minmax_sign}{abs(data.numerics.minmax)}{minmax_string}{fom_string}", ) process_output.ocmmnt( outfile, - f"Convergence parameter : {data_structure.numerics.epsvmc}", + f"Convergence parameter : {data.numerics.epsvmc}", ) process_output.oblnkl(outfile) @@ -214,12 +209,12 @@ def run_summary(data: DataStructure): process_output.ovarst(mfile, "Input filename", "(fileprefix)", f'"{fileprefix}"') process_output.ovarin( - mfile, "Optimisation switch", "(ioptimz)", data_structure.numerics.ioptimz + mfile, "Optimisation switch", "(ioptimz)", data.numerics.ioptimz ) # If optimising, write figure of merit switch - if data_structure.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: + if data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: process_output.ovarin( - mfile, "Figure of merit switch", "(minmax)", data_structure.numerics.minmax + mfile, "Figure of merit switch", "(minmax)", data.numerics.minmax ) @@ -231,7 +226,6 @@ def init_all_module_vars(): than a 'run-once' executable. """ logging_model_handler.clear_logs() - data_structure.numerics.init_numerics() constants.init_constants() @@ -243,23 +237,23 @@ def check_process(inputs, data): # noqa: ARG001 and ensures other dependent variables are given suitable values. """ # Check that there are sufficient iteration variables - if data_structure.numerics.nvar < data_structure.numerics.neqns: + if data.numerics.nvar < data.numerics.neqns: raise ProcessValidationError( "Insufficient iteration variables to solve the problem! NVAR < NEQNS", - nvar=data_structure.numerics.nvar, - neqns=data_structure.numerics.neqns, + nvar=data.numerics.nvar, + neqns=data.numerics.neqns, ) # Check that sufficient elements of ixc and icc have been specified - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 0).any(): + if (data.numerics.ixc[: data.numerics.nvar] == 0).any(): raise ProcessValidationError( "The number of iteration variables specified is smaller than the number stated in ixc", - nvar=data_structure.numerics.nvar, + nvar=data.numerics.nvar, ) # Check that dr_tf_wp_with_insulation (ixc = 140) and dr_tf_inboard (ixc = 13) are not being used simultaneously as iteration variables - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 13).any() and ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 140 + if (data.numerics.ixc[: data.numerics.nvar] == 13).any() and ( + data.numerics.ixc[: data.numerics.nvar] == 140 ).any(): raise ProcessValidationError( "Iteration variables 13 and 140 cannot be used simultaneously", @@ -267,38 +261,31 @@ def check_process(inputs, data): # noqa: ARG001 # Can't use c_tf_turn as interation var, constraint or input if i_tf_turns_integer == 1 if ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 60 + data.numerics.ixc[: data.numerics.nvar] == 60 ).any() and data.tfcoil.i_tf_turns_integer == 1: raise ProcessValidationError( "Iteration variable 60 (TF current per turn, c_tf_turn) cannot be used with the TF coil integer turn model (i_tf_turns_integer == 1) as it is a calculated output instead for this model. However, the maximum current per turn can be constrained with constraint 77." ) # Can't have icc 77 and ixc 60 at the same time - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 60).any() and ( - data_structure.numerics.icc[: data_structure.numerics.nvar] == 77 + if (data.numerics.ixc[: data.numerics.nvar] == 60).any() and ( + data.numerics.icc[: data.numerics.nvar] == 77 ).any(): raise ProcessValidationError( "Cannot use iteration variable 60 (TF coil current per turn, c_tf_turn) and constraint 77 (maximum TF current per turn) simultaneously." ) - if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 0 - ).any(): + if (data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 0).any(): raise ProcessValidationError( "The number of constraints specified is smaller than the number stated in neqns+nineqns", - neqns=data_structure.numerics.neqns, - nineqns=data_structure.numerics.nineqns, + neqns=data.numerics.neqns, + nineqns=data.numerics.nineqns, ) # Deprecate constraints for depcrecated_constraint in [3, 4, 10, 74, 42]: if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == depcrecated_constraint ).any(): raise ProcessValidationError( @@ -307,10 +294,7 @@ def check_process(inputs, data): # noqa: ARG001 # MDK Report error if constraint 63 is used with old vacuum model if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 63 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 63 ).any() and data.vacuum.i_vacuum_pumping != "simple": raise ProcessValidationError( "Constraint 63 is requested without the correct vacuum model (simple)" @@ -352,7 +336,7 @@ def check_process(inputs, data): # noqa: ARG001 # Stop the run if j_tf_coil_full_area is used as an optimisation variable # As the current density is now calculated from b_plasma_toroidal_on_axis without constraint 10 - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 12).any(): + if (data.numerics.ixc[: data.numerics.nvar] == 12).any(): raise ProcessValidationError( "The 1/R toroidal B field dependency constraint is being depreciated" ) @@ -404,20 +388,17 @@ def check_process(inputs, data): # noqa: ARG001 ) if ( - data_structure.numerics.ioptimz == PROCESSRunMode.OPTIMISATION - and (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 4).any() - and data_structure.numerics.boundl[3] - < data.physics.temp_plasma_pedestal_kev * 1.001 + data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION + and (data.numerics.ixc[: data.numerics.nvar] == 4).any() + and data.numerics.boundl[3] < data.physics.temp_plasma_pedestal_kev * 1.001 ): warn( "Lower limit of volume averaged electron temperature (temp_plasma_electron_vol_avg_kev) has been raised to ensure temp_plasma_electron_vol_avg_kev > temp_plasma_pedestal_kev", stacklevel=2, ) - data_structure.numerics.boundl[3] = ( - data.physics.temp_plasma_pedestal_kev * 1.001 - ) - data_structure.numerics.boundu[3] = max( - data_structure.numerics.boundu[3], data_structure.numerics.boundl[3] + data.numerics.boundl[3] = data.physics.temp_plasma_pedestal_kev * 1.001 + data.numerics.boundu[3] = max( + data.numerics.boundu[3], data.numerics.boundl[3] ) # Density checks @@ -468,22 +449,17 @@ def check_process(inputs, data): # noqa: ARG001 # Issue #862 : Variable nd_plasma_electron_on_axis/nd_plasma_pedestal_electron ratio without constraint eq 81 (nd_plasma_electron_on_axis>nd_plasma_pedestal_electron) # -> Potential hollowed density profile if ( - data_structure.numerics.ioptimz == PROCESSRunMode.OPTIMISATION + data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION and not ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 81 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 81 ).any() ): - if ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 145 - ).any(): + if (data.numerics.ixc[: data.numerics.nvar] == 145).any(): warn( "nd_plasma_pedestal_electron set with f_nd_plasma_pedestal_greenwald without constraint eq 81 (nd_plasma_pedestal_electron 273.15 K" ) @@ -632,8 +592,8 @@ def check_process(inputs, data): # noqa: ARG001 # Check if conductor upper limit is properly set to 50 K or below if ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 20 - ).any() and data_structure.numerics.boundu[19] > 50.0: + data.numerics.ixc[: data.numerics.nvar] == 20 + ).any() and data.numerics.boundu[19] > 50.0: raise ProcessValidationError( "Too large CP conductor temperature (temp_cp_average). Upper limit for cryo-al < 50 K" ) @@ -662,10 +622,7 @@ def check_process(inputs, data): # noqa: ARG001 if ( data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM and ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 85 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 85 ).any() and data.physics.itart == 1 ): @@ -685,7 +642,7 @@ def check_process(inputs, data): # noqa: ARG001 # Checking the CP TF top radius if ( abs(data.build.r_cp_top) > 0 - or (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 174).any() + or (data.numerics.ixc[: data.numerics.nvar] == 174).any() ) and data.build.i_r_cp_top != 1: raise ProcessValidationError( "To set the TF CP top value, you must use i_r_cp_top = 1" @@ -734,10 +691,7 @@ def check_process(inputs, data): # noqa: ARG001 # Constraint 10 is dedicated to ST designs with demountable joints if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 10 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 10 ).any(): raise ProcessValidationError( "Constraint equation 10 (CP lifetime) to used with ST desing (itart=1)" @@ -758,13 +712,9 @@ def check_process(inputs, data): # noqa: ARG001 if ( ( not ( - (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 16).any() - or ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 29 - ).any() - or ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 42 - ).any() + (data.numerics.ixc[: data.numerics.nvar] == 16).any() + or (data.numerics.ixc[: data.numerics.nvar] == 29).any() + or (data.numerics.ixc[: data.numerics.nvar] == 42).any() ) ) # No dr_bore,dr_cs_tf_gap, dr_cs iteration and ( @@ -778,16 +728,10 @@ def check_process(inputs, data): # noqa: ARG001 ) # dr_bore + dr_cs_tf_gap + dr_cs = 0 and ( ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 31 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 31 ).any() or ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 32 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 32 ).any() ) # Stress constraints (31 or 32) is used and ( @@ -939,7 +883,7 @@ def check_process(inputs, data): # noqa: ARG001 # Check if the WP/conductor radial thickness (dr_tf_wp_with_insulation) is large enough # To contains the insulation, cooling and the support structure # Rem : Only verified if the WP thickness is used - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 140).any(): + if (data.numerics.ixc[: data.numerics.nvar] == 140).any(): # Minimal WP thickness if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: dr_tf_wp_min = 2.0 * ( @@ -950,8 +894,8 @@ def check_process(inputs, data): # noqa: ARG001 ) # Steel conduit thickness (can be an iteration variable) - if (data_structure.numerics.ixc[: data_structure.numerics.nvar] == 58).any(): - dr_tf_wp_min += 2.0 * data_structure.numerics.boundl[57] + if (data.numerics.ixc[: data.numerics.nvar] == 58).any(): + dr_tf_wp_min += 2.0 * data.numerics.boundl[57] else: dr_tf_wp_min += 2.0 * data.tfcoil.dx_tf_turn_steel @@ -966,7 +910,7 @@ def check_process(inputs, data): # noqa: ARG001 + 4.0 * data.tfcoil.radius_cp_coolant_channel ) - if data_structure.numerics.boundl[139] < dr_tf_wp_min: + if data.numerics.boundl[139] < dr_tf_wp_min: raise ProcessValidationError( "The TF coil WP thickness (dr_tf_wp_with_insulation) must be at least", dr_tf_wp_min=dr_tf_wp_min, @@ -1095,10 +1039,7 @@ def check_process(inputs, data): # noqa: ARG001 # Cannot use temperature margin constraint with REBCO TF coils if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 36 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 36 ).any() and ( SuperconductorModel(data.tfcoil.i_tf_sc_mat).sc_type == SuperconductorMaterial.REBCO @@ -1109,10 +1050,7 @@ def check_process(inputs, data): # noqa: ARG001 # Cannot use temperature margin constraint with REBCO CS coils if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 60 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 60 ).any() and data.pf_coil.i_cs_superconductor == 8: raise ProcessValidationError( "turn off CS temperature margin constraint icc = 60 when using REBCO" @@ -1124,32 +1062,27 @@ def check_process(inputs, data): # noqa: ARG001 # Cannot use TF coil strain limit if i_str_wp is off: if ( - data_structure.numerics.icc[ - : data_structure.numerics.neqns + data_structure.numerics.nineqns - ] - == 88 + data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] == 88 ).any() and data.tfcoil.i_str_wp == 0: raise ProcessValidationError("Can't use constraint 88 if i_strain_tf == 0") -def set_active_constraints(): +def set_active_constraints(data: DataStructure): """Set constraints provided in the input file as 'active'""" num_constraints = 0 for i in range(ConstraintManager.num_constraints()): - if data_structure.numerics.icc[i] != 0: - data_structure.numerics.active_constraints[ - data_structure.numerics.icc[i] - 1 - ] = True + if data.numerics.icc[i] != 0: + data.numerics.active_constraints[data.numerics.icc[i] - 1] = True num_constraints += 1 - if data_structure.numerics.neqns < 0: + if data.numerics.neqns < 0: # The value of neqns has not been set in the input file. Default = 0. - data_structure.numerics.neqns = num_constraints - data_structure.numerics.nineqns + data.numerics.neqns = num_constraints - data.numerics.nineqns else: - data_structure.numerics.nineqns = num_constraints - data_structure.numerics.neqns + data.numerics.nineqns = num_constraints - data.numerics.neqns -def set_device_type(data): +def set_device_type(data: DataStructure): if data.ife.ife == 1: data.globals.icase = "Inertial Fusion model" elif data.stellarator.istell != 0: diff --git a/process/core/input.py b/process/core/input.py index 97b524cfda..664f331d55 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -10,13 +10,13 @@ from warnings import warn import process -from process import data_structure from process.core.exceptions import ( ProcessValidationError, ProcessValueError, ) from process.core.solver.constraints import ConstraintManager from process.data_structure.impurity_radiation_variables import N_IMPURITIES +from process.data_structure.numerics import IPEQNS, IPNVARS from process.data_structure.pfcoil_variables import N_PF_GROUPS_MAX from process.data_structure.physics_variables import N_CONFINEMENT_SCALINGS from process.data_structure.scan_variables import IPNSCNS, IPNSCNV @@ -33,14 +33,18 @@ """Valid input variable types alias""" -def _ixc_additional_actions(_name, value: int, _array_index, _config): - data_structure.numerics.ixc[data_structure.numerics.nvar] = value - data_structure.numerics.nvar += 1 +def _ixc_additional_actions( + _name, value: int, _array_index, _config, data: DataStructure +): + data.numerics.ixc[data.numerics.nvar] = value + data.numerics.nvar += 1 -def _icc_additional_actions(_name, value: int, _array_index, _config): - data_structure.numerics.icc[data_structure.numerics.n_constraints] = value - data_structure.numerics.n_constraints += 1 +def _icc_additional_actions( + _name, value: int, _array_index, _config, data: DataStructure +): + data.numerics.icc[data.numerics.n_constraints] = value + data.numerics.n_constraints += 1 @dataclass @@ -108,18 +112,18 @@ def __post_init__(self): "runtitle": InputVariable("globals", str), "verbose": InputVariable("globals", int, choices=[0, 1]), "run_tests": InputVariable("globals", int, choices=[0, 1]), - "ioptimz": InputVariable(data_structure.numerics, int, choices=[1, -2]), - "epsvmc": InputVariable(data_structure.numerics, float, range=(0.0, 1.0)), - "boundl": InputVariable(data_structure.numerics, float, array=True), - "boundu": InputVariable(data_structure.numerics, float, array=True), - "epsfcn": InputVariable(data_structure.numerics, float, range=(0.0, 1.0)), + "ioptimz": InputVariable("numerics", int, choices=[1, -2]), + "epsvmc": InputVariable("numerics", float, range=(0.0, 1.0)), + "boundl": InputVariable("numerics", float, array=True), + "boundu": InputVariable("numerics", float, array=True), + "epsfcn": InputVariable("numerics", float, range=(0.0, 1.0)), "maxcal": InputVariable("globals", int, range=(0, 10000)), - "minmax": InputVariable(data_structure.numerics, int), + "minmax": InputVariable("numerics", int), "neqns": InputVariable( - data_structure.numerics, int, range=(0, ConstraintManager.num_constraints()) + "numerics", int, range=(0, ConstraintManager.num_constraints()) ), "nineqns": InputVariable( - data_structure.numerics, int, range=(0, ConstraintManager.num_constraints()) + "numerics", int, range=(0, ConstraintManager.num_constraints()) ), "alphaj": InputVariable("physics", float, range=(0.0, 10.0)), "alphan": InputVariable("physics", float, range=(0.0, 10.0)), @@ -583,9 +587,7 @@ def __post_init__(self): ), "nd_plasma_pedestal_electron": InputVariable("physics", float, range=(0.0, 1e21)), "nd_plasma_separatrix_electron": InputVariable("physics", float, range=(0.0, 1e21)), - "i_nd_plasma_pedestal_separatrix": InputVariable( - data_structure.physics_variables, int, choices=[0, 1] - ), + "i_nd_plasma_pedestal_separatrix": InputVariable("physics", int, choices=[0, 1]), "nflutfmax": InputVariable("constraints", float, range=(0.0, 1e24)), "j_tf_coil_full_area": InputVariable("tfcoil", float, range=(10000.0, 1000000000.0)), "f_a_cs_turn_steel": InputVariable("pf_coil", float, range=(0.001, 0.999)), @@ -1123,24 +1125,24 @@ def __post_init__(self): "ixc": InputVariable( None, int, - range=(1, data_structure.numerics.ipnvars), + range=(1, IPNVARS), additional_actions=_ixc_additional_actions, set_variable=False, ), "icc": InputVariable( None, int, - range=(1, data_structure.numerics.ipeqns), + range=(1, IPEQNS), additional_actions=_icc_additional_actions, set_variable=False, ), "force_vmcon_inequality_satisfication": InputVariable( - data_structure.numerics, + "numerics", int, choices=(0, 1), ), "force_vmcon_inequality_tolerance": InputVariable( - data_structure.numerics, float, range=(0.0, 1e10) + "numerics", float, range=(0.0, 1e10) ), } @@ -1257,6 +1259,7 @@ def parse_input_file(data_structure_obj: DataStructure): clean_variable_value, array_index_clean, variable_config, + data_structure_obj, ) # add the variable to a dictionary indexed by the variable name (in the input file) diff --git a/process/core/io/mfile/base.py b/process/core/io/mfile/base.py index f43f92306d..59532ed469 100644 --- a/process/core/io/mfile/base.py +++ b/process/core/io/mfile/base.py @@ -13,7 +13,6 @@ import numpy as np -from process import data_structure from process.core.model import DataStructure from process.core.solver import iteration_variables @@ -603,8 +602,8 @@ def get_mfile_initial_ixc_values(file_path: Path, data: DataStructure): iteration_variable_names = [] iteration_variable_values = [] - for i in range(data_structure.numerics.nvar): - ivar = data_structure.numerics.ixc[i].item() + for i in range(data.numerics.nvar): + ivar = data.numerics.ixc[i].item() itv = iteration_variables.ITERATION_VARIABLES[ivar] module = getattr(data, itv.module) if isinstance(itv.module, str) else itv.module diff --git a/process/core/io/plot/solutions.py b/process/core/io/plot/solutions.py index f4fce3a19d..d3760e75ea 100644 --- a/process/core/io/plot/solutions.py +++ b/process/core/io/plot/solutions.py @@ -21,7 +21,6 @@ import seaborn as sns from process.core.io.mfile import MFile -from process.data_structure import numerics from process.data_structure.numerics import FiguresOfMerit # Variables of interest in mfiles and subsequent dataframes @@ -445,7 +444,6 @@ def _plot_solutions( ): objf_list = norm_objf_df[NORM_OBJF_NAME].unique() else: - numerics.init_numerics() objf_list = { FiguresOfMerit(abs(int(minmax))).description for minmax in diffs_df["minmax"] } diff --git a/process/core/io/vary_run/config.py b/process/core/io/vary_run/config.py index fef0bab0d5..4c6904d389 100644 --- a/process/core/io/vary_run/config.py +++ b/process/core/io/vary_run/config.py @@ -146,7 +146,9 @@ def __iter__(self): def __next__(self): _neqns, itervars = get_neqns_itervars(in_dat=self.infile, wdir=self.wdir) - lbs, ubs = get_variable_range(itervars, self.factor, self.infile, self.wdir) + lbs, ubs = get_variable_range( + itervars, self.factor, self.infile, self.data, self.wdir + ) self.run_process(self.wdir / self.infile, self.solver) check_input_error(mfile=self.outfile, wdir=self.wdir) @@ -263,6 +265,7 @@ def modify_in_dat(self): def setup(self, data: DataStructure): """Sets up the program for running""" + self.data = data self.echo() self.prepare_wdir() @@ -275,8 +278,10 @@ def setup(self, data: DataStructure): self.generator = default_rng(seed=self.u_seed) - data.globals.output_prefix = str(self.wdir / self.outfile.strip("MFILE.DAT")) - data.globals.fileprefix = str(self.wdir / self.infile) + self.data.globals.output_prefix = str( + self.wdir / self.outfile.strip("MFILE.DAT") + ) + self.data.globals.fileprefix = str(self.wdir / self.infile) @staticmethod def run_process(input_path: Path, solver: str = "vmcon"): diff --git a/process/core/io/vary_run/tools.py b/process/core/io/vary_run/tools.py index 641367f6c3..9d871213a7 100644 --- a/process/core/io/vary_run/tools.py +++ b/process/core/io/vary_run/tools.py @@ -10,7 +10,7 @@ from process.core.io.data_structure_dicts import get_dicts from process.core.io.in_dat import InDat from process.core.io.mfile import MFile -from process.data_structure import numerics +from process.core.model import DataStructure logger = logging.getLogger(__name__) @@ -55,7 +55,7 @@ def update_ixc_bounds(indat, wdir="."): dicts["DICT_IXC_BOUNDS"][name]["ub"] = float(value["u"]) -def get_variable_range(itervars, factor, indat, wdir="."): +def get_variable_range(itervars, factor, indat, data: DataStructure, wdir="."): """Returns the lower and upper bounds of the variable range for each iteration variable. @@ -70,6 +70,10 @@ def get_variable_range(itervars, factor, indat, wdir="."): setting them to value * factor and value / factor respectively while taking their process bounds into account. + indat : str + input file name + data : DataStructure + data structure object wdir : (Default value = ".") """ @@ -80,12 +84,12 @@ def get_variable_range(itervars, factor, indat, wdir="."): lbs = [] ubs = [] - iteration_variables = numerics.lablxc + iteration_variables = data.numerics.lablxc for varname in itervars: iteration_variable_index = iteration_variables.index(varname) - lb = numerics.boundl[iteration_variable_index] - ub = numerics.boundu[iteration_variable_index] + lb = data.numerics.boundl[iteration_variable_index] + ub = data.numerics.boundu[iteration_variable_index] # for f-values we set the same range as in process if varname[0] == "f" and (varname not in dicts["NON_F_VALUES"]): lbs += [lb] diff --git a/process/core/model.py b/process/core/model.py index e4fa92fb05..3ddc7bd3ac 100644 --- a/process/core/model.py +++ b/process/core/model.py @@ -19,6 +19,7 @@ from process.data_structure.ife_variables import IFEData from process.data_structure.impurity_radiation_variables import ImpurityRadiationData from process.data_structure.neoclassics_variables import NeoclassicsData +from process.data_structure.numerics import NumericsData from process.data_structure.pf_power_variables import PFPowerData from process.data_structure.pfcoil_variables import PFCoilData from process.data_structure.physics_variables import PhysicsData @@ -79,6 +80,7 @@ class DataStructure: superconducting_tfcoil: SuperconductingTFData = initialise_later globals: GlobalData = initialise_later scan: ScanData = initialise_later + numerics: NumericsData = initialise_later def __post_init__(self): for f in fields(self): diff --git a/process/core/process_output.py b/process/core/process_output.py index fb8340c617..3faff8c0c2 100644 --- a/process/core/process_output.py +++ b/process/core/process_output.py @@ -4,7 +4,6 @@ import numpy as np from process.core import constants -from process.data_structure import numerics class OutputFileManager: @@ -183,10 +182,12 @@ def ovarre(file, descr: str, varnam: str, value, output_flag: str = ""): format_value = f"{value:.17e}" if isinstance(value, float) else f"{value: >12}" - if varnam.strip("()") in numerics.name_xc: - # MDK add ITV label if it is an iteration variable - # The ITV flag overwrites the output_flag - output_flag = "ITV" + # TODO need to find a way to identify iteration variables at a higher level + # in the data structure + # if varnam.strip("()") in numerics.name_xc: + # # MDK add ITV label if it is an iteration variable + # # The ITV flag overwrites the output_flag + # output_flag = "ITV" line = f"{description}{replacement_character} {varname}{replacement_character} {format_value} {output_flag}" write(file, line) diff --git a/process/core/scan.py b/process/core/scan.py index 51a6c11693..b2da1471e0 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -15,7 +15,6 @@ 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 import numerics from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData @@ -257,7 +256,9 @@ def post_optimise(self, ifail: int): ifail: int : """ - numerics.sqsumsq = sum(r**2 for r in numerics.rcm[: numerics.neqns]) ** 0.5 + 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": @@ -301,7 +302,7 @@ def post_optimise(self, ifail: int): process_output.oblnkl(constants.NOUT) process_output.ovarin(constants.NOUT, "Error flag", "(ifail)", ifail) - if numerics.sqsumsq >= 1.0e-2: + if self.data.numerics.sqsumsq >= 1.0e-2: process_output.oblnkl(constants.NOUT) process_output.ocmmnt( constants.NOUT, @@ -330,44 +331,56 @@ def post_optimise(self, ifail: int): ) process_output.oblnkl(constants.IOTTY) - logger.warning(f"High final constraint residues. {numerics.sqsumsq=}") + logger.warning( + f"High final constraint residues. {self.data.numerics.sqsumsq=}" + ) process_output.ovarin( - constants.NOUT, "Number of iteration variables", "(nvar)", numerics.nvar + constants.NOUT, + "Number of iteration variables", + "(nvar)", + self.data.numerics.nvar, ) process_output.ovarin( constants.NOUT, "Number of constraints (total)", "(neqns+nineqns)", - numerics.neqns + numerics.nineqns, + self.data.numerics.neqns + self.data.numerics.nineqns, ) process_output.ovarin( - constants.NOUT, "Optimisation switch", "(ioptimz)", numerics.ioptimz + constants.NOUT, + "Optimisation switch", + "(ioptimz)", + self.data.numerics.ioptimz, ) process_output.ocmmnt( - constants.NOUT, f" {PROCESSRunMode(numerics.ioptimz).description}" + 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)", numerics.minmax + constants.NOUT, + "Figure of merit switch", + "(minmax)", + self.data.numerics.minmax, ) - objf_name = f'"{FiguresOfMerit(abs(numerics.minmax)).description}"' + objf_name = f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"' - numerics.objf_name = objf_name + self.data.numerics.objf_name = objf_name process_output.ovarst( constants.NOUT, "Objective function name", "(objf_name)", - numerics.objf_name, + self.data.numerics.objf_name, ) process_output.ovarre( constants.NOUT, "Normalised objective function", "(norm_objf)", - numerics.norm_objf, + self.data.numerics.norm_objf, "OP ", ) @@ -375,7 +388,7 @@ def post_optimise(self, ifail: int): constants.NOUT, "Square root of the sum of squares of the constraint residuals", "(sqsumsq)", - numerics.sqsumsq, + self.data.numerics.sqsumsq, "OP ", ) if self.solver != "fsolve": @@ -390,7 +403,7 @@ def post_optimise(self, ifail: int): constants.NOUT, "Number of optimising solver iterations", "(nviter)", - numerics.nviter, + self.data.numerics.nviter, "OP ", ) process_output.oblnkl(constants.NOUT) @@ -410,7 +423,7 @@ def post_optimise(self, ifail: int): else: string1 = "PROCESS has failed to optimise" - string2 = "minimise" if numerics.minmax > 0 else "maximise" + string2 = "minimise" if self.data.numerics.minmax > 0 else "maximise" process_output.write( constants.NOUT, @@ -421,17 +434,23 @@ def post_optimise(self, ifail: int): # Output optimisation parameters solution_vector_table = [] - for i in range(numerics.nvar): - numerics.xcs[i] = numerics.xcm[i] * numerics.scafc[i] + for i in range(self.data.numerics.nvar): + self.data.numerics.xcs[i] = ( + self.data.numerics.xcm[i] * self.data.numerics.scafc[i] + ) - name = numerics.lablxc[numerics.ixc[i] - 1] - solution_vector_table.append([name, numerics.xcs[i], numerics.xcm[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 * numerics.itv_scaled_lower_bounds[i] - xmaxx = 0.99 * numerics.itv_scaled_upper_bounds[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 numerics.xcm[i] < xminn or numerics.xcm[i] > xmaxx: + if self.data.numerics.xcm[i] < xminn or self.data.numerics.xcm[i] > xmaxx: if not written_warning: written_warning = True process_output.ocmmnt( @@ -443,37 +462,40 @@ def post_optimise(self, ifail: int): ), ) - xcval = numerics.xcm[i] * numerics.scafc[i] + xcval = self.data.numerics.xcm[i] * self.data.numerics.scafc[i] - if numerics.xcm[i] < xminn: + if self.data.numerics.xcm[i] < xminn: location, bound = "below", "lower" - bounds = numerics.itv_scaled_lower_bounds + bounds = self.data.numerics.itv_scaled_lower_bounds else: location, bound = "above", "upper" - bounds = numerics.itv_scaled_upper_bounds + 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] * numerics.scafc[i]}", + f" {bounds[i] * self.data.numerics.scafc[i]}", ) # Write optimisation parameters to mfile process_output.ovarre( constants.MFILE, - numerics.lablxc[numerics.ixc[i] - 1], + self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1], f"(itvar{i + 1:03d})", - numerics.xcs[i], + self.data.numerics.xcs[i], ) - if numerics.boundu[i] == numerics.boundl[i]: + if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]: xnorm = 1.0 else: xnorm = min( max( - (numerics.xcm[i] - numerics.itv_scaled_lower_bounds[i]) + ( + self.data.numerics.xcm[i] + - self.data.numerics.itv_scaled_lower_bounds[i] + ) / ( - numerics.itv_scaled_upper_bounds[i] - - numerics.itv_scaled_lower_bounds[i] + self.data.numerics.itv_scaled_upper_bounds[i] + - self.data.numerics.itv_scaled_lower_bounds[i] ), 0.0, ), @@ -484,7 +506,7 @@ def post_optimise(self, ifail: int): constants.MFILE, f"{name} (final value/initial value)", f"(xcm{i + 1:03d})", - numerics.xcm[i], + self.data.numerics.xcm[i], ) process_output.ovarre( constants.MFILE, @@ -496,13 +518,15 @@ def post_optimise(self, ifail: int): constants.MFILE, f"{name} (upper bound)", f"(boundu{i + 1:03d})", - numerics.itv_scaled_upper_bounds[i] * numerics.scafc[i], + 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})", - numerics.itv_scaled_lower_bounds[i] * numerics.scafc[i], + self.data.numerics.itv_scaled_lower_bounds[i] + * self.data.numerics.scafc[i], ) # Write optimisation parameter headings to output file @@ -524,13 +548,13 @@ def post_optimise(self, ifail: int): ) con1, con2, err, _, lab = constraints.constraint_eqns( - numerics.neqns + numerics.nineqns, -1, self.data + self.data.numerics.neqns + self.data.numerics.nineqns, -1, self.data ) # Write equality constraints to mfile equality_constraint_table = [] - for i in range(numerics.neqns): - name = numerics.lablcc[numerics.icc[i] - 1] + for i in range(self.data.numerics.neqns): + name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1] equality_constraint_table.append([ name, @@ -542,27 +566,27 @@ def post_optimise(self, ifail: int): process_output.ovarre( constants.MFILE, f"{name:<33} normalised residue", - f"(eq_con{numerics.icc[i]:03d})", + f"(eq_con{self.data.numerics.icc[i]:03d})", con1[i], ) process_output.ovarre( constants.MFILE, f"{name:<33} residual", - f"(res_eq_con{numerics.icc[i]:03d})", + 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{numerics.icc[i]:03d})", + f"(val_eq_con{self.data.numerics.icc[i]:03d})", con2[i], ) process_output.ovarre( constants.MFILE, f"{name} units", - f"(eq_units_con{numerics.icc[i]:03d})", + f"(eq_units_con{self.data.numerics.icc[i]:03d})", f"'{lab[i]}'", ) @@ -583,7 +607,7 @@ def post_optimise(self, ifail: int): ) # Write inequality constraints - if numerics.nineqns > 0: + if self.data.numerics.nineqns > 0: inequality_constraint_table = [] # Inequalities not necessarily satisfied when evaluating process_output.osubhd( @@ -597,10 +621,13 @@ def post_optimise(self, ifail: int): "might be violated.", ) - for i in range(numerics.neqns, numerics.neqns + numerics.nineqns): - name = numerics.lablcc[numerics.icc[i] - 1] + 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(numerics.icc[i]), self.data + int(self.data.numerics.icc[i]), self.data ) inequality_constraint_table.append([ @@ -614,34 +641,34 @@ def post_optimise(self, ifail: int): process_output.ovarre( constants.MFILE, f"{name} normalised residue", - f"(ineq_con{numerics.icc[i]:03d})", + 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{numerics.icc[i]:03d})", + 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{numerics.icc[i]:03d})", + 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{numerics.icc[i]:03d})", + 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{numerics.icc[i]:03d})", + f"(ineq_bound_con{self.data.numerics.icc[i]:03d})", constraint.constraint_bound, ) @@ -1087,13 +1114,13 @@ def scan_select(self, nwp, swp, iscn): case 9: self.data.physics.temp_plasma_electron_vol_avg_kev = swp[iscn - 1] case 10: - numerics.boundu[14] = swp[iscn - 1] + 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: - numerics.boundu[9] = swp[iscn - 1] + self.data.numerics.boundu[9] = swp[iscn - 1] case 16: self.data.physics.rmajor = swp[iscn - 1] case 17: @@ -1101,7 +1128,7 @@ def scan_select(self, nwp, swp, iscn): case 18: self.data.constraints.eta_cd_norm_hcd_primary_max = swp[iscn - 1] case 19: - numerics.boundl[15] = swp[iscn - 1] + self.data.numerics.boundl[15] = swp[iscn - 1] case 20: self.data.constraints.t_burn_min = swp[iscn - 1] case 22: @@ -1125,13 +1152,13 @@ def scan_select(self, nwp, swp, iscn): case 31: self.data.constraints.f_alpha_energy_confinement_min = swp[iscn - 1] case 32: - numerics.epsvmc = swp[iscn - 1] + self.data.numerics.epsvmc = swp[iscn - 1] case 38: - numerics.boundu[128] = swp[iscn - 1] + self.data.numerics.boundu[128] = swp[iscn - 1] case 39: - numerics.boundu[130] = swp[iscn - 1] + self.data.numerics.boundu[130] = swp[iscn - 1] case 40: - numerics.boundu[134] = swp[iscn - 1] + self.data.numerics.boundu[134] = swp[iscn - 1] case 41: self.data.build.dr_blkt_outboard = swp[iscn - 1] case 42: @@ -1144,7 +1171,7 @@ def scan_select(self, nwp, swp, iscn): case 45: self.data.tfcoil.temp_tf_superconductor_margin_min = swp[iscn - 1] case 46: - numerics.boundu[151] = swp[iscn - 1] + self.data.numerics.boundu[151] = swp[iscn - 1] case 48: self.data.tfcoil.n_tf_wp_pancakes = int(swp[iscn - 1]) case 49: @@ -1159,7 +1186,7 @@ def scan_select(self, nwp, swp, iscn): case 52: self.data.physics.rad_fraction_sol = swp[iscn - 1] case 53: - numerics.boundu[156] = swp[iscn - 1] + self.data.numerics.boundu[156] = swp[iscn - 1] case 54: self.data.tfcoil.b_crit_upper_nbti = swp[iscn - 1] case 55: @@ -1167,7 +1194,7 @@ def scan_select(self, nwp, swp, iscn): case 56: self.data.heat_transport.p_cryo_plant_electric_max_mw = swp[iscn - 1] case 57: - numerics.boundl[1] = swp[iscn - 1] + self.data.numerics.boundl[1] = swp[iscn - 1] case 58: self.data.build.dr_fw_plasma_gap_inboard = swp[iscn - 1] case 59: diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 5b294cb0ba..4f33c2a930 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -5,7 +5,6 @@ import numpy as np -from process import data_structure from process.core import constants from process.core.exceptions import ProcessError, ProcessValueError from process.core.model import DataStructure @@ -1802,7 +1801,7 @@ def constraint_eqns(m: int, ieqn: int, data: DataStructure): cc, con, err, symbol, units = [], [], [], [], [] for i in range(i1, i2): - constraint_id = data_structure.numerics.icc[i] + constraint_id = data.numerics.icc[i] result = ConstraintManager.evaluate_constraint(constraint_id, data) tmp_cc, tmp_con, tmp_err = ( diff --git a/process/core/solver/evaluators.py b/process/core/solver/evaluators.py index 8677e0edd6..9138b6a3c9 100644 --- a/process/core/solver/evaluators.py +++ b/process/core/solver/evaluators.py @@ -5,7 +5,6 @@ from process.core.caller import Caller from process.core.model import DataStructure -from process.data_structure import numerics logger = logging.getLogger(__name__) @@ -70,9 +69,9 @@ def fcnvmc1(self, _n, m, xv, ifail): sqsumconfsq = math.sqrt(summ) logger.debug("Key evaluator values:") - logger.debug(f"{numerics.nviter = }") + logger.debug(f"{self.data.numerics.nviter = }") logger.debug(f"{(1 - (ifail % 7)) - 1 = }") - logger.debug(f"{(numerics.nviter % 2) - 1 = }") + logger.debug(f"{(self.data.numerics.nviter % 2) - 1 = }") logger.debug(f"{self.data.physics.temp_plasma_electron_vol_avg_kev = }") logger.debug(f"{self.data.costs.coe = }") logger.debug(f"{self.data.physics.rmajor = }") @@ -127,8 +126,8 @@ def fcnvmc2(self, n, m, xv, lcnorm): xfor[j] = xv[j] xbac[j] = xv[j] if i == j: - xfor[i] = xv[j] * (1.0 + numerics.epsfcn) - xbac[i] = xv[j] * (1.0 - numerics.epsfcn) + xfor[i] = xv[j] * (1.0 + self.data.numerics.epsfcn) + xbac[i] = xv[j] * (1.0 - self.data.numerics.epsfcn) # Evaluate at (x+dx) ffor, cfor = self.caller.call_models(xfor, m) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 26ef603195..06a1e7ac34 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -5,7 +5,6 @@ import numpy as np -from process import data_structure from process.core.exceptions import ProcessValueError from process.core.model import DataStructure @@ -259,8 +258,8 @@ def check_iteration_variable(iteration_variable_value, name: str = ""): def load_iteration_variables(data): """Loads the physics and engineering variables into the optimisation variable array.""" - for i in range(data_structure.numerics.nvar): - variable_index = data_structure.numerics.ixc[i] + for i in range(data.numerics.nvar): + variable_index = data.numerics.ixc[i] iteration_variable = ITERATION_VARIABLES[variable_index] # use ... as the default return value because None might be a valid return from Fortran? @@ -293,9 +292,9 @@ def load_iteration_variables(data): iteration_variable.array_index ] - data_structure.numerics.xcm[i] = iteration_variable_value - data_structure.numerics.name_xc[i] = ( - data_structure.numerics.name_xc[i], + data.numerics.xcm[i] = iteration_variable_value + data.numerics.name_xc[i] = ( + data.numerics.name_xc[i], iteration_variable.name, ) @@ -320,12 +319,10 @@ def load_iteration_variables(data): name=f"{variable_index} ({iteration_variable.name})", ) - data_structure.numerics.scale[i] = 1.0 / iteration_variable_value - data_structure.numerics.scafc[i] = 1.0 / data_structure.numerics.scale[i] + data.numerics.scale[i] = 1.0 / iteration_variable_value + data.numerics.scafc[i] = 1.0 / data.numerics.scale[i] - data_structure.numerics.xcm[i] = ( - iteration_variable_value * data_structure.numerics.scale[i] - ) + data.numerics.xcm[i] = iteration_variable_value * data.numerics.scale[i] def set_scaled_iteration_variable(xc, nn: int, data: DataStructure): @@ -344,10 +341,10 @@ def set_scaled_iteration_variable(xc, nn: int, data: DataStructure): # there is less error handling here than in load_iteration_variables # because many errors will be caught in load_iteration_variables which is # run first. This verifies the variables exist and the module target is correct. - variable_index = data_structure.numerics.ixc[i] + variable_index = data.numerics.ixc[i] iteration_variable = ITERATION_VARIABLES[variable_index] - ratio = xc[i] / data_structure.numerics.scale[i] + ratio = xc[i] / data.numerics.scale[i] module = ( getattr(data, iteration_variable.module) @@ -380,24 +377,22 @@ def set_scaled_iteration_variable(xc, nn: int, data: DataStructure): ) -def load_scaled_bounds(): +def load_scaled_bounds(data: DataStructure): """Sets the scaled bounds of the iteration variables.""" - for i in range(data_structure.numerics.nvar): - variable_index = data_structure.numerics.ixc[i] - 1 - data_structure.numerics.itv_scaled_lower_bounds[i] = ( - data_structure.numerics.boundl[variable_index] - * data_structure.numerics.scale[i] + for i in range(data.numerics.nvar): + variable_index = data.numerics.ixc[i] - 1 + data.numerics.itv_scaled_lower_bounds[i] = ( + data.numerics.boundl[variable_index] * data.numerics.scale[i] ) - data_structure.numerics.itv_scaled_upper_bounds[i] = ( - data_structure.numerics.boundu[variable_index] - * data_structure.numerics.scale[i] + data.numerics.itv_scaled_upper_bounds[i] = ( + data.numerics.boundu[variable_index] * data.numerics.scale[i] ) -def initialise_iteration_variables(): +def initialise_iteration_variables(data: DataStructure): """Initialise the iteration variables (label and default bounds)""" for itv_index, itv in ITERATION_VARIABLES.items(): - data_structure.numerics.lablxc[itv_index - 1] = itv.name + data.numerics.lablxc[itv_index - 1] = itv.name - data_structure.numerics.boundl[itv_index - 1] = itv.lower_bound - data_structure.numerics.boundu[itv_index - 1] = itv.upper_bound + data.numerics.boundl[itv_index - 1] = itv.lower_bound + data.numerics.boundu[itv_index - 1] = itv.upper_bound diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index 167ff4b2cb..2c3d8eced3 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -19,7 +19,6 @@ from process.core.exceptions import ProcessValueError from process.core.model import DataStructure from process.core.solver.evaluators import Evaluators -from process.data_structure import numerics logger = logging.getLogger(__name__) @@ -27,11 +26,12 @@ class _Solver(ABC): """Base class for different solver implementations.""" - def __init__(self, *, maxcal: int | None = None): + def __init__(self, *, data: DataStructure, maxcal: int | None = None): """Initialise a solver.""" # Exit code for the solver self.ifail = 0 - self.tolerance = numerics.epsvmc + self.data = data + self.tolerance = self.data.numerics.epsvmc self.b: float | None = None self.convergence_parameter: float | None = None self.maxcal = maxcal @@ -182,10 +182,10 @@ def solve(self) -> int: bb = None if self.b is not None: - bb = np.identity(numerics.nvar) * self.b + bb = np.identity(self.data.numerics.nvar) * self.b def _solver_callback(i: int, _result, _x, convergence_param: float): - numerics.nviter = i + 1 + self.data.numerics.nviter = i + 1 self.convergence_parameter = convergence_param print( f"{i + 1} | Convergence Parameter: {convergence_param:.3E}", @@ -227,7 +227,9 @@ def _ineq_cons_satisfied( # negative constraint value = violated # Check all ineqs are satisfied to within the tolerance # E.g. the relative violations are no more than v=0-tolerance - return bool(np.all(result.ie >= -numerics.force_vmcon_inequality_tolerance)) + return bool( + np.all(result.ie >= -self.data.numerics.force_vmcon_inequality_tolerance) + ) try: x, _, _, res = solve( @@ -241,7 +243,7 @@ def _ineq_cons_satisfied( initial_B=bb, callback=_solver_callback, additional_convergence=_ineq_cons_satisfied - if numerics.force_vmcon_inequality_satisfication + if self.data.numerics.force_vmcon_inequality_satisfication else None, ) except VMCONConvergenceException as e: @@ -259,8 +261,10 @@ def _ineq_cons_satisfied( except ValueError: itervar_name_list = "" - for count, iter_var in enumerate(numerics.ixc[: numerics.nvar]): - itervar_name = numerics.lablxc[iter_var - 1] + for count, iter_var in enumerate( + self.data.numerics.ixc[: self.data.numerics.nvar] + ): + itervar_name = self.data.numerics.lablxc[iter_var - 1] itervar_name_list += f"{count}: {itervar_name} \n" logger.warning(f"Active iteration variables are : \n{itervar_name_list}") @@ -359,11 +363,13 @@ def get_solver(data: DataStructure, solver_name: str = "vmcon") -> _Solver: solver: _Solver if solver_name == "vmcon": - solver = Vmcon(maxcal=data.globals.maxcal) + solver = Vmcon(data=data, maxcal=data.globals.maxcal) elif solver_name == "vmcon_bounded": - solver = VmconBounded(maxcal=data.globals.maxcal) + solver = VmconBounded(data=data, maxcal=data.globals.maxcal) elif solver_name == "fsolve": - solver = FSolve() + solver = FSolve( + data=data, + ) else: try: solver = load_external_solver(solver_name) diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index 03e76a687a..b457e25ca9 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -4,7 +4,6 @@ load_scaled_bounds, ) from process.core.solver.solver import get_solver -from process.data_structure import numerics class SolverHandler: @@ -32,19 +31,19 @@ def run(self): """Run solver and retry if it fails in certain ways.""" # Initialise iteration variables and bounds in Fortran load_iteration_variables(self.data) - load_scaled_bounds() + load_scaled_bounds(self.data) # 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 = numerics.nvar - x = numerics.xcm[:n] - bndl = numerics.itv_scaled_lower_bounds[:n] - bndu = numerics.itv_scaled_upper_bounds[:n] + 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 = numerics.neqns + numerics.nineqns - meq = numerics.neqns + m = self.data.numerics.neqns + self.data.numerics.nineqns + meq = self.data.numerics.neqns # Evaluators() calculates the objective and constraint functions and # their gradients for a given vector x @@ -64,26 +63,26 @@ def run(self): 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) + # 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 - numerics.epsfcn /= 10 # reset value + self.data.numerics.epsfcn /= 10 # reset value if ifail != 1: print("Trying again with new epsfcn") - numerics.epsfcn /= 10 # try new smaller value - print("new epsfcn = ", numerics.epsfcn) + self.data.numerics.epsfcn /= 10 # try new smaller value + print("new epsfcn = ", self.data.numerics.epsfcn) ifail = self.solver.solve() - numerics.epsfcn *= 10 # reset value + self.data.numerics.epsfcn *= 10 # reset value # 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) # Only do this if VMCON has not iterated (nviter=1) - if ifail == 5 and numerics.nviter < 2: + if ifail == 5 and self.data.numerics.nviter < 2: print( "VMCON error code = 5. Rerunning VMCON with a new initial " "estimate of the second derivative matrix." @@ -96,12 +95,12 @@ def run(self): return ifail def output(self): - """Store results back in Fortran numerics module. + """Store results back in Fortran self.data.numerics module. Objective function value, solution vector and constraints vector. """ - numerics.norm_objf = self.solver.objf + self.data.numerics.norm_objf = self.solver.objf # Slicing required due to Fortran arrays being maximum possible, rather # than required, size - numerics.xcm[: self.solver.x.shape[0]] = self.solver.x - numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf + self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x + self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 8499302c80..427467d30e 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -1,3 +1,4 @@ +from dataclasses import dataclass, field from enum import IntEnum from types import DynamicClassAttribute @@ -99,603 +100,526 @@ def description(self): return self._description_ -ipnvars: int = 177 +IPNVARS = 177 """total number of variables available for iteration""" -ipeqns: int = 92 +IPEQNS = 92 """number of constraint equations available""" -ipnfoms: int = 19 +IPNFOMS = 19 """number of available figures of merit""" -ipvlam: int = ipeqns + 2 * ipnvars + 1 -iptnt: int = (ipeqns * (3 * ipeqns + 13)) / 2 -ipvp1: int = ipnvars + 1 - -ioptimz: int = None -"""Code operation switch: -* -2 for evaluation mode (i.e. no optimisation) -* 1 for optimisation mode (e.g. via VMCON) -""" - -minmax: int = None -""" -Switch for figure-of-merit (see `FiguresOfMerit` for descriptions) -negative => maximise, positive => minimise -""" - -n_constraints: int = None -"""Total number of constraints (neqns + nineqns)""" - - -ncalls: int = None -"""number of function calls during solution""" - -neqns: int = None -"""number of equality constraints to be satisfied""" - -nfev1: int = None -"""number of calls to FCNHYB (HYBRD function caller) made""" - -nfev2: int = None -"""number of calls to FCNVMC1 (VMCON function caller) made""" - -nineqns: int = None -"""number of inequality constraints VMCON must satisfy -(leave at zero for now) -""" - -nvar: int = None -"""number of iteration variables to use""" - -nviter: int = None -"""number of optimisation iterations performed""" - -icc: list[int] = None - -active_constraints: list[bool] = True -"""Logical array showing which constraints are active""" - -# TODO Do not change the comments for lablcc: they are used to create the -# Python-Fortran dictionaries. This must be improved on. - -lablcc: list[str] = None -"""Labels describing constraint equations (corresponding itvs)