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 384e2e3323..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: @@ -161,15 +160,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 +201,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 +229,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( @@ -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 79cbd5e575..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 @@ -21,7 +20,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 @@ -45,16 +43,16 @@ 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() + process_output.OutputFileManager.open_files(data.globals.output_prefix) # Input any desired new initial values inputs = parse_input_file(data) # Set active constraints - set_active_constraints() + set_active_constraints(data) # set the device type (icase) set_device_type(data) @@ -65,7 +63,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 +101,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 +137,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,49 +150,45 @@ 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) 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_structure.global_variables.maxcal}", + 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) @@ -215,12 +209,12 @@ def run_summary(): 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 ) @@ -232,9 +226,6 @@ def init_all_module_vars(): than a 'run-once' executable. """ logging_model_handler.clear_logs() - data_structure.numerics.init_numerics() - data_structure.global_variables.init_global_variables() - init_scan_variables() constants.init_constants() @@ -246,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", @@ -270,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( @@ -310,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)" @@ -355,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" ) @@ -407,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 @@ -471,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" ) @@ -635,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" ) @@ -665,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 ): @@ -688,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" @@ -737,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)" @@ -748,7 +699,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 @@ -761,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 ( @@ -781,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 ( @@ -942,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 * ( @@ -953,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 @@ -969,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, @@ -1098,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 @@ -1112,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" @@ -1127,33 +1062,28 @@ 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_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..664f331d55 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -10,15 +10,16 @@ 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 if TYPE_CHECKING: from collections.abc import Callable @@ -32,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 @@ -68,11 +73,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`. @@ -103,21 +109,21 @@ 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]), - "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)), - "minmax": InputVariable(data_structure.numerics, int), + "runtitle": InputVariable("globals", str), + "verbose": InputVariable("globals", int, choices=[0, 1]), + "run_tests": InputVariable("globals", int, choices=[0, 1]), + "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("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)), @@ -581,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)), @@ -1049,7 +1053,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 +1095,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, @@ -1121,30 +1125,30 @@ 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) ), } 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 != "": @@ -1255,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 a289cc13b3..4c6904d389 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__) @@ -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) @@ -261,8 +263,9 @@ 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.data = data self.echo() self.prepare_wdir() @@ -275,10 +278,10 @@ def setup(self): self.generator = default_rng(seed=self.u_seed) - data_structure.global_variables.output_prefix = str( + self.data.globals.output_prefix = str( self.wdir / self.outfile.strip("MFILE.DAT") ) - data_structure.global_variables.fileprefix = str(self.wdir / self.infile) + 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 dff9ccc1f1..3ddc7bd3ac 100644 --- a/process/core/model.py +++ b/process/core/model.py @@ -14,10 +14,12 @@ 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 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 @@ -26,6 +28,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 @@ -75,6 +78,9 @@ class DataStructure: rebco: RebcoData = initialise_later tfcoil: TFData = initialise_later 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 ff8956bd8b..3faff8c0c2 100644 --- a/process/core/process_output.py +++ b/process/core/process_output.py @@ -4,38 +4,37 @@ import numpy as np from process.core import constants -from process.data_structure import global_variables, 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): @@ -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 30c157676b..b2da1471e0 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -15,12 +15,8 @@ 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 ( - global_variables, - numerics, - scan_variables, -) 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 @@ -217,7 +213,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() @@ -231,14 +227,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() @@ -260,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": @@ -304,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, @@ -333,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 ", ) @@ -378,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": @@ -386,14 +396,14 @@ 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( constants.NOUT, "Number of optimising solver iterations", "(nviter)", - numerics.nviter, + self.data.numerics.nviter, "OP ", ) process_output.oblnkl(constants.NOUT) @@ -413,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, @@ -424,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( @@ -446,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, ), @@ -487,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, @@ -499,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 @@ -527,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, @@ -545,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]}'", ) @@ -586,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( @@ -600,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([ @@ -617,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, ) @@ -826,7 +850,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() @@ -842,13 +866,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 @@ -857,7 +881,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( @@ -871,23 +895,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() @@ -906,13 +930,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 @@ -920,7 +944,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): @@ -931,8 +955,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( @@ -949,53 +973,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): - global_variables.iscan_global = iscan - global_variables.vlabel, global_variables.xlabel = self.scan_select( - scan_variables.nsweep, scan_variables.sweep, iscan + self.data.globals.iscan_global = iscan + self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( + self.data.scan.nsweep, self.data.scan.sweep, iscan ) process_output.oblnkl(constants.NOUT) @@ -1003,8 +1027,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 {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) @@ -1012,22 +1036,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"{global_variables.xlabel} , {global_variables.vlabel}" - f" = {scan_variables.sweep[iscan - 1]}" + f"Starting scan point {iscan} of {self.data.scan.isweep} : " + f"{self.data.globals.xlabel} , {self.data.globals.vlabel}" + f" = {self.data.scan.sweep[iscan - 1]}" ) def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): - iscan_r = 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) - global_variables.iscan_global = iscan + self.data.globals.iscan_global = iscan - global_variables.vlabel, global_variables.xlabel = self.scan_select( - scan_variables.nsweep, scan_variables.sweep, iscan_1 + self.data.globals.vlabel, self.data.globals.xlabel = self.scan_select( + self.data.scan.nsweep, self.data.scan.sweep, iscan_1 ) - global_variables.vlabel_2, global_variables.xlabel_2 = self.scan_select( - scan_variables.nsweep_2, scan_variables.sweep_2, iscan_r + self.data.globals.vlabel_2, self.data.globals.xlabel_2 = self.scan_select( + self.data.scan.nsweep_2, self.data.scan.sweep_2, iscan_r ) process_output.oblnkl(constants.NOUT) @@ -1035,9 +1059,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"{global_variables.vlabel} = {scan_variables.sweep[iscan_1 - 1]} and" - f" {global_variables.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,31 +1069,31 @@ 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} = {self.data.scan.sweep[iscan_1 - 1]}" + f" and {self.data.globals.xlabel_2}, " + f"{self.data.globals.vlabel_2} = {self.data.scan.sweep_2[iscan_r - 1]} " ) 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: @@ -1090,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: @@ -1104,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: @@ -1128,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: @@ -1147,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: @@ -1162,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: @@ -1170,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 d892bb6e61..9138b6a3c9 100644 --- a/process/core/solver/evaluators.py +++ b/process/core/solver/evaluators.py @@ -5,8 +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,16 +62,16 @@ 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 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 = }") @@ -128,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 4b858c7a33..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,17 +292,17 @@ 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, ) # 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( ( @@ -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 a285de39c7..2c3d8eced3 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -17,8 +17,8 @@ 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 logger = logging.getLogger(__name__) @@ -26,12 +26,15 @@ class _Solver(ABC): """Base class for different solver implementations.""" - def __init__(self): + 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 def set_evaluators(self, evaluators: Evaluators): """Set objective and constraint functions and their gradient evaluators. @@ -179,11 +182,11 @@ 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 - global_variables.convergence_parameter = convergence_param + self.data.numerics.nviter = i + 1 + self.convergence_parameter = convergence_param print( f"{i + 1} | Convergence Parameter: {convergence_param:.3E}", end="\r", @@ -224,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( @@ -232,13 +237,13 @@ 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, 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: @@ -256,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}") @@ -340,7 +347,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,11 +363,13 @@ def get_solver(solver_name: str = "vmcon") -> _Solver: solver: _Solver if solver_name == "vmcon": - solver = Vmcon() + solver = Vmcon(data=data, maxcal=data.globals.maxcal) elif solver_name == "vmcon_bounded": - solver = VmconBounded() + 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 5094c5a9e2..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,26 +31,26 @@ 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 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) @@ -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/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/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)