From 633b8145c4a24b8f6299dc72a2de57dcc11cb0b5 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 1/2] 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)