From 1615b5e9f407473dd5a7fc5b584a93c5bfd1a3c0 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:28:19 +0100 Subject: [PATCH 01/10] rework solve slightly --- process/core/caller.py | 100 ++++++++++++++++++- process/core/final.py | 105 -------------------- process/core/output.py | 153 ----------------------------- process/core/solver/constraints.py | 123 ++++++++++++++++++++++- process/core/solver/solver.py | 60 +++++++++++ process/main.py | 144 ++++++++++++++++++++++++++- 6 files changed, 423 insertions(+), 262 deletions(-) delete mode 100644 process/core/final.py delete mode 100644 process/core/output.py diff --git a/process/core/caller.py b/process/core/caller.py index e618e51262..1634777e53 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -9,13 +9,14 @@ from tabulate import tabulate from process.core import constants -from process.core.final import finalise +from process.core import process_output as po from process.core.io.mfile import MFile from process.core.process_output import OutputFileManager, ovarre from process.core.solver import constraints from process.core.solver.iteration_variables import set_scaled_iteration_variable from process.core.solver.objectives import objective_function from process.data_structure.blanket_variables import BlktModelTypes +from process.data_structure.numerics import PROCESSRunMode from process.models.tfcoil.base import TFConductorModel from process.models.tfcoil.superconducting import SuperconductingTFTurnType @@ -396,6 +397,103 @@ def _call_models_once(self, xc: np.ndarray): # FISPACT and LOCA model (not used)- removed +def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None): + """Routine to print out the final point in the scan. + + Writes to OUT.DAT and MFILE.DAT. + + Parameters + ---------- + models : process.main.Models + physics and engineering model objects + data: DataStructure + data structure object to provide data to evaluate the constraints + ifail : int + error flag + non_idempotent_msg : None | str, optional + warning about non-idempotent variables, defaults to None + """ + if ifail == 1: + po.oheadr(constants.NOUT, "Final Feasible Point") + else: + po.oheadr(constants.NOUT, "Final UNFEASIBLE Point") + + # Output relevant to no optimisation + if data.numerics.ioptimz == PROCESSRunMode.EVALUATION: + output_evaluation(data) + + # Print non-idempotence warning to OUT.DAT only + if non_idempotent_msg: + po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES") + po.ocmmnt(constants.NOUT, non_idempotent_msg) + + # Write output to OUT.DAT and MFILE.DAT + models.write(data, constants.NOUT) + + +def output_evaluation(data): + """Write output for an evaluation run of PROCESS + + Parameters + ---------- + data: DataStructure + data structure object to provide data to evaluate the constraints + """ + po.oheadr(constants.NOUT, "Numerics") + po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.") + po.oblnkl(constants.NOUT) + + # Evaluate objective function + norm_objf = objective_function(data.numerics.minmax, data) + po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf) + + # Print the residuals of the constraint equations + + residual_error, value, residual, symbols, units = constraints.constraint_eqns( + data.numerics.neqns + data.numerics.nineqns, -1, data + ) + + labels = [ + data.numerics.lablcc[j - 1] + for j in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] + ] + + def _fmt(a, units): + return [f"{c} {u}" for c, u in zip(a, units, strict=False)] + + po.write( + constants.NOUT, + tabulate( + { + "Constraint Name": labels, + "Constraint Type": symbols, + "Physical constraint": _fmt(value, units), + "Constraint residual": _fmt(residual, units), + "Normalised residual": residual_error, + }, + headers="keys", + ), + ) + + for i in range(data.numerics.neqns): + constraint_id = data.numerics.icc[i] + po.ovarre( + constants.MFILE, + f"{labels[i]} normalised residue", + f"(eq_con{constraint_id:03d})", + residual_error[i], + ) + + for i in range(data.numerics.nineqns): + constraint_id = data.numerics.icc[data.numerics.neqns + i] + po.ovarre( + constants.MFILE, + f"{labels[data.numerics.neqns + i]}", + f"(ineq_con{constraint_id:03d})", + residual_error[data.numerics.neqns + i], + ) + + def write_output_files( models: Models, data: DataStructure, ifail: int, *, runtime: float | None = None ): diff --git a/process/core/final.py b/process/core/final.py deleted file mode 100644 index 55a42e8fc2..0000000000 --- a/process/core/final.py +++ /dev/null @@ -1,105 +0,0 @@ -"""Final output at the end of a scan.""" - -from tabulate import tabulate - -from process.core import constants -from process.core import output as op -from process.core import process_output as po -from process.core.solver import constraints -from process.core.solver.objectives import objective_function -from process.data_structure.numerics import PROCESSRunMode, SolverOutputCondition - - -def finalise(models, data, ifail: int, non_idempotent_msg: str | None = None): - """Routine to print out the final point in the scan. - - Writes to OUT.DAT and MFILE.DAT. - - Parameters - ---------- - models : process.main.Models - physics and engineering model objects - data: DataStructure - data structure object to provide data to evaluate the constraints - ifail : int - error flag - non_idempotent_msg : None | str, optional - warning about non-idempotent variables, defaults to None - """ - if ifail == SolverOutputCondition.CONVERGED: - po.oheadr(constants.NOUT, "Final Feasible Point") - else: - po.oheadr(constants.NOUT, "Final UNFEASIBLE Point") - - # Output relevant to no optimisation - if data.numerics.ioptimz == PROCESSRunMode.EVALUATION: - output_evaluation(data) - - # Print non-idempotence warning to OUT.DAT only - if non_idempotent_msg: - po.oheadr(constants.NOUT, "NON-IDEMPOTENT VARIABLES") - po.ocmmnt(constants.NOUT, non_idempotent_msg) - - # Write output to OUT.DAT and MFILE.DAT - op.write(models, data, constants.NOUT) - - -def output_evaluation(data): - """Write output for an evaluation run of PROCESS - - Parameters - ---------- - data: DataStructure - data structure object to provide data to evaluate the constraints - """ - po.oheadr(constants.NOUT, "Numerics") - po.ocmmnt(constants.NOUT, "PROCESS has performed an evaluation run.") - po.oblnkl(constants.NOUT) - - # Evaluate objective function - norm_objf = objective_function(data.numerics.minmax, data) - po.ovarre(constants.MFILE, "Normalised objective function", "(norm_objf)", norm_objf) - - # Print the residuals of the constraint equations - - residual_error, value, residual, symbols, units = constraints.constraint_eqns( - data.numerics.neqns + data.numerics.nineqns, -1, data - ) - - labels = [ - data.numerics.lablcc[j] - for j in [ - i - 1 - for i in data.numerics.icc[: data.numerics.neqns + data.numerics.nineqns] - ] - ] - physical_constraint = [f"{c} {u}" for c, u in zip(value, units, strict=False)] - physical_residual = [f"{c} {u}" for c, u in zip(residual, units, strict=False)] - - table_data = { - "Constraint Name": labels, - "Constraint Type": symbols, - "Physical constraint": physical_constraint, - "Constraint residual": physical_residual, - "Normalised residual": residual_error, - } - - po.write(constants.NOUT, tabulate(table_data, headers="keys")) - - for i in range(data.numerics.neqns): - constraint_id = data.numerics.icc[i] - po.ovarre( - constants.MFILE, - f"{labels[i]} normalised residue", - f"(eq_con{constraint_id:03d})", - residual_error[i], - ) - - for i in range(data.numerics.nineqns): - constraint_id = data.numerics.icc[data.numerics.neqns + i] - po.ovarre( - constants.MFILE, - f"{labels[data.numerics.neqns + i]}", - f"(ineq_con{constraint_id:03d})", - residual_error[data.numerics.neqns + i], - ) diff --git a/process/core/output.py b/process/core/output.py deleted file mode 100644 index 770c1f486d..0000000000 --- a/process/core/output.py +++ /dev/null @@ -1,153 +0,0 @@ -"""Module containing routine to write the results to the main output file (OUT.DAT)""" - -from process.core.log import logging_model_handler -from process.data_structure.blanket_variables import BlktModelTypes -from process.models.tfcoil.base import TFConductorModel -from process.models.tfcoil.superconducting import ( - SuperconductingTFTurnType, -) - - -def write(models, data, _outfile): - """Write the results to the main output file (OUT.DAT). - - Write the program results to a file, in a tidy format. - - Parameters - ---------- - models : process.main.Models - physics and engineering model objects - _outfile : int - Fortran output unit identifier - - Raises - ------ - ValueError - If unsupported superconducting TF turn type is used - """ - # ensure we are capturing warnings that occur in the 'output' stage - # as these are warnings that occur at our solution point. - # So we clear existing warnings - logging_model_handler.start_capturing() - logging_model_handler.clear_logs() - - # Call stellarator output routine instead if relevant - if data.stellarator.istell != 0: - models.stellarator.output() - return - - # Call IFE output routine instead if relevant - if data.ife.ife != 0: - models.ife.output() - return - - # Costs model - # Cost switch values - # No. | model - # ---- | ------ - # 0 | 1990 costs model - # 1 | 2015 Kovari model - # 2 | Custom model - models.costs.output() - - # Availability model - models.availability.output() - - # Physics model - models.physics.output() - - # Detailed physics, currently only done at final point as values are not used - # by any other functions - models.physics_detailed.output() - - # TODO what is this? Not in caller.py? - models.current_drive.output() - - # Pulsed reactor model - models.pulse.output() - - models.divertor.output() - - # Machine Build Model - models.build.output() - - # Cryostat build - models.cryostat.output() - - # Toroidal field coil copper model - if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER: - models.copper_tf_coil.output() - - # Toroidal field coil superconductor model - if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: - tf_turn_type = SuperconductingTFTurnType( - data.superconducting_tfcoil.i_tf_turn_type - ) - if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT: - models.cicc_sctfcoil.output() - elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR: - models.croco_sctfcoil.output() - else: - raise ValueError( - "Unsupported superconducting TF turn type: " - f"{data.superconducting_tfcoil.i_tf_turn_type}" - ) - - # Toroidal field coil aluminium model - if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM: - models.aluminium_tf_coil.output() - - # Tight aspect ratio machine model - if ( - data.physics.itart == 1 - and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING - ): - models.tfcoil.output() - - # Poloidal field coil model - models.pfcoil.output() - - # Structure Model - models.structure.output() - - # Blanket model - # Blanket switch values - # No. | model - # ---- | ------ - # 1 | CCFE HCPB model - # 2 | KIT HCPB model - # 3 | CCFE HCPB model with Tritium Breeding Ratio calculation - # 4 | KIT HCLL model - # 5 | DCLL model - - models.shield.output() - models.vacuum_vessel.output() - - # First wall geometry - models.fw.output() - - if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB: - # CCFE HCPB model - models.ccfe_hcpb.output() - - elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL: - # DCLL model - models.dcll.output() - - # FISPACT and LOCA model (not used)- removed - - # Power model - models.power.output() - - # Vacuum model - models.vacuum.output() - - # Buildings model - models.buildings.output() - - # Water usage in secondary cooling system - models.water_use.output() - - # stop capturing warnings so that Outfile does not end up with - # a lot of non-model logs - logging_model_handler.stop_capturing() diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index e8def90ef6..dbef355604 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -6,8 +6,9 @@ from typing import ClassVar, Literal import numpy as np +from tabulate import tabulate -from process.core import constants +from process.core import constants, process_output from process.core.exceptions import ProcessError, ProcessValueError from process.core.model import DataStructure from process.data_structure.build_variables import TFCSRadialConfiguration @@ -2010,3 +2011,123 @@ def constraint_eqns(m: int, ieqn: int, data: DataStructure): units.append(tmp_units) return np.array(cc), np.array(con), np.array(err), symbol, units + + +def constraints_output(data: DataStructure, solver_name: str): + nums = data.numerics + + process_output.osubhd( + constants.NOUT, + "The following equality constraint residues should be close to zero:", + ) + + con1, con2, err, _, lab = constraint_eqns(nums.neqns + nums.nineqns, -1, data) + + # Write equality constraints to mfile + equality_constraint_table = [] + for i in range(nums.neqns): + name = nums.lablcc[nums.icc[i] - 1] + + equality_constraint_table.append([ + name, + "=", + f"{con2[i]} {lab[i]}", + f"{err[i]} {lab[i]}", + con1[i], + ]) + + for d, var, v in ( + (f"{name:<33} normalised residue", f"(eq_con{nums.icc[i]:03d})", con1[i]), + (f"{name:<33} residual", f"(res_eq_con{nums.icc[i]:03d})", err[i]), + (f"{name} constraint value", f"(val_eq_con{nums.icc[i]:03d})", con2[i]), + (f"{name} units", f"(eq_units_con{nums.icc[i]:03d})", f"'{lab[i]}'"), + ): + process_output.ovarre(constants.MFILE, d, var, v) + + # Write equality constraints to output file + process_output.write( + constants.NOUT, + tabulate( + equality_constraint_table, + headers=[ + "", + "", + "Physical constraint", + "Constraint residue", + "Normalised residue", + ], + numalign="left", + ), + ) + + # Write inequality constraints + if nums.nineqns > 0: + inequality_constraint_table = [] + # Inequalities not necessarily satisfied when evaluating + process_output.osubhd( + constants.NOUT, + "Negative inequality constraint (normalised) residuals indicate a constraint is satisfied.", + ) + if solver_name == "fsolve": + process_output.osubhd( + constants.NOUT, + "This MFile was produced via an evaluation, not an optimisation, and so the constraints " + "might be violated.", + ) + + for i in range( + nums.neqns, + nums.neqns + nums.nineqns, + ): + name = nums.lablcc[nums.icc[i] - 1] + constraint = ConstraintManager.evaluate_constraint(int(nums.icc[i]), data) + + inequality_constraint_table.append([ + name, + f"{constraint.constraint_value} {constraint.units}", + constraint.symbol, + f"{constraint.constraint_bound} {constraint.units}", + f"{constraint.residual} {constraint.units}", + f"{constraint.normalised_residual}", + ]) + + for d, var, v in ( + ( + "normalised residue", + f"(ineq_con{nums.icc[i]:03d})", + -constraint.normalised_residual, + ), + ( + "physical value", + f"(ineq_value_con{nums.icc[i]:03d})", + constraint.constraint_value, + ), + ( + "symbol", + f"(ineq_symbol_con{nums.icc[i]:03d})", + f"'{constraint.symbol}'", + ), + ("units", f"(ineq_units_con{nums.icc[i]:03d})", f"'{constraint.units}'"), + ( + "physical bound", + f"(ineq_bound_con{nums.icc[i]:03d})", + constraint.constraint_bound, + ), + ): + process_output.ovarre(constants.MFILE, f"{name} {d}", var, v) + + process_output.write( + constants.NOUT, + tabulate( + inequality_constraint_table, + headers=[ + "", + "Physical constraint", + "", + "Physical constraint bound", + "Constraint residue", + "Normalised residue", + ], + numalign="left", + ), + ) diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index d96d14a7c4..dd3b0bd6f4 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -16,6 +16,7 @@ ) from scipy.optimize import fsolve +from process.core import constants, process_output from process.core.exceptions import ProcessValueError from process.core.model import DataStructure from process.core.solver.evaluators import Evaluators @@ -294,6 +295,65 @@ def _ineq_cons_satisfied( return self.info + def verror(self): + """Routine to print out relevant messages in the case of an + unfeasible result from a VMCON (optimisation) run + + This routine prints out relevant messages in the case of + an unfeasible result from a VMCON (optimisation) run. + + Parameters + ---------- + ifail: int : + + """ + strings = "\n".join( + { + -1: ("User-terminated execution of VMCON.",), + 0: ( + "Improper input parameters to the VMCON routine.", + "PROCESS coding must be checked.", + ), + 2: ( + "The maximum number of calls has been reached without solution.", + ( + "The code may be stuck in a minimum in the residual space that" + " is significantly above zero.\n" + ), + "There is either no solution possible, or the code", + "is failing to escape from a deep local minimum.", + "Try changing the variables in IXC, or modify their initial values.", + ), + 3: ( + "The line search required the maximum of 10 calls.", + "A feasible solution may be difficult to achieve.", + "Try changing or adding variables to IXC.", + ), + 4: ( + "An uphill search direction was found.", + "Try changing the equations in ICC, or", + "adding new variables to IXC.", + ), + 5: ( + "The quadratic programming technique was unable to", + "find a feasible point.\n", + "Try changing or adding variables to IXC, or modify", + "their initial values (especially if only 1 optimisation", + "iteration was performed).", + ), + 6: ( + "The quadratic programming technique was restricted", + "by an artificial bound, or failed due to a singular", + "matrix.", + "Try changing the equations in ICC, or", + "adding new variables to IXC.", + ), + }.get(self.info, "Unknown Error code") + ) + + process_output.ocmmnt(constants.NOUT, strings) + print(strings) + class VmconBounded(Vmcon): """A solver that uses VMCON but checks x is in bounds before running""" diff --git a/process/main.py b/process/main.py index 4d3cf62342..ed5c903b54 100644 --- a/process/main.py +++ b/process/main.py @@ -55,6 +55,7 @@ from process.core.model import DataStructure, Model from process.core.process_output import OutputFileManager, oheadr from process.core.scan import Scan +from process.data_structure.blanket_variables import BlktModelTypes from process.data_structure.cost_variables import CostModels from process.data_structure.numerics import PROCESSRunMode from process.models.availability import Availability @@ -111,7 +112,7 @@ from process.models.stellarator.neoclassics import Neoclassics from process.models.stellarator.stellarator import Stellarator from process.models.structure import Structure -from process.models.tfcoil.base import TFCoil +from process.models.tfcoil.base import TFCoil, TFConductorModel from process.models.tfcoil.resistive import ( AluminiumTFCoil, CopperTFCoil, @@ -121,6 +122,7 @@ CICCSuperconductingTFCoil, CROCOSuperconductingTFCoil, SuperconductingTFCoil, + SuperconductingTFTurnType, ) from process.models.vacuum import Vacuum, VacuumVessel from process.models.water_use import WaterUse @@ -454,7 +456,6 @@ def run_scan(self): # ioptimz == 1: optimisation if self.data.numerics.ioptimz == PROCESSRunMode.OPTIMISATION: pass - # ioptimz == -2: evaluation elif self.data.numerics.ioptimz == PROCESSRunMode.EVALUATION: # No optimisation: # solve equality (consistency) constraints only using fsolve (HYBRD) @@ -841,6 +842,145 @@ def setup_data_structure(self): for model in self.models: model.data = self.data + def write(self, data, _outfile): + """Write the results to the main output file (OUT.DAT). + + Write the program results to a file, in a tidy format. + + Parameters + ---------- + self : process.main.Models + physics and engineering model objects + _outfile : int + Fortran output unit identifier + + """ + # ensure we are capturing warnings that occur in the 'output' stage as these are warnings + # that occur at our solution point. So we clear existing warnings + logging_model_handler.start_capturing() + logging_model_handler.clear_logs() + + # Call stellarator output routine instead if relevant + if data.stellarator.istell != 0: + self.stellarator.output() + return + + # Call IFE output routine instead if relevant + if data.ife.ife != 0: + self.ife.output() + return + + # Costs model + # Cost switch values + # No. | model + # ---- | ------ + # 0 | 1990 costs model + # 1 | 2015 Kovari model + # 2 | Custom model + self.costs.output() + + # Availability model + self.availability.output() + + # Physics model + self.physics.output() + + # Detailed physics, currently only done at final point as values are not used + # by any other functions + self.physics_detailed.output() + + # TODO what is this? Not in caller.py? + self.current_drive.output() + + # Pulsed reactor model + self.pulse.output() + + self.divertor.output() + + # Machine Build Model + self.build.output() + + # Cryostat build + self.cryostat.output() + + # Toroidal field coil copper model + if data.tfcoil.i_tf_sup == TFConductorModel.WATER_COOLED_COPPER: + self.copper_tf_coil.output() + + # Toroidal field coil superconductor model + if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: + tf_turn_type = SuperconductingTFTurnType( + data.superconducting_tfcoil.i_tf_turn_type + ) + if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT: + self.cicc_sctfcoil.output() + elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR: + self.croco_sctfcoil.output() + else: + raise ValueError( + "Unsupported superconducting TF turn type: " + f"{data.superconducting_tfcoil.i_tf_turn_type}" + ) + + # Toroidal field coil aluminium model + if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM: + self.aluminium_tf_coil.output() + + # Tight aspect ratio machine model + if ( + data.physics.itart == 1 + and data.tfcoil.i_tf_sup != TFConductorModel.SUPERCONDUCTING + ): + self.tfcoil.output() + + # Poloidal field coil model + self.pfcoil.output() + + # Structure Model + self.structure.output() + + # Blanket model + # Blanket switch values + # No. | model + # ---- | ------ + # 1 | CCFE HCPB model + # 2 | KIT HCPB model + # 3 | CCFE HCPB model with Tritium Breeding Ratio calculation + # 4 | KIT HCLL model + # 5 | DCLL model + + self.shield.output() + self.vacuum_vessel.output() + + # First wall geometry + self.fw.output() + + if data.fwbs.i_blanket_type == BlktModelTypes.CCFE_HCPB: + # CCFE HCPB model + self.ccfe_hcpb.output() + + elif data.fwbs.i_blanket_type == BlktModelTypes.DCLL: + # DCLL model + self.dcll.output() + + # FISPACT and LOCA model (not used)- removed + + # Power model + self.power.output() + + # Vacuum model + self.vacuum.output() + + # Buildings model + self.buildings.output() + + # Water usage in secondary cooling system + self.water_use.output() + + # stop capturing warnings so that Outfile does not end up with + # a lot of non-model logs + logging_model_handler.stop_capturing() + # setup handlers for writing to terminal (on warnings+) # or writing to the log file (on info+) From 65f3fc0e3a9115106cda33e33686bb7e8860e8dc Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:43:44 +0100 Subject: [PATCH 02/10] fix linting errors --- process/core/solver/constraints.py | 7 ++++--- process/core/solver/solver.py | 5 ----- process/main.py | 5 +++-- 3 files changed, 7 insertions(+), 10 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index dbef355604..0ce2c67b01 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -2066,13 +2066,14 @@ def constraints_output(data: DataStructure, solver_name: str): # Inequalities not necessarily satisfied when evaluating process_output.osubhd( constants.NOUT, - "Negative inequality constraint (normalised) residuals indicate a constraint is satisfied.", + "Negative inequality constraint (normalised) residuals " + "indicate a constraint is satisfied.", ) if solver_name == "fsolve": process_output.osubhd( constants.NOUT, - "This MFile was produced via an evaluation, not an optimisation, and so the constraints " - "might be violated.", + "This MFile was produced via an evaluation, not an optimisation, " + "and so the constraints might be violated.", ) for i in range( diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index dd3b0bd6f4..846e525e3d 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -301,11 +301,6 @@ def verror(self): This routine prints out relevant messages in the case of an unfeasible result from a VMCON (optimisation) run. - - Parameters - ---------- - ifail: int : - """ strings = "\n".join( { diff --git a/process/main.py b/process/main.py index ed5c903b54..2ca60408e7 100644 --- a/process/main.py +++ b/process/main.py @@ -855,8 +855,9 @@ def write(self, data, _outfile): Fortran output unit identifier """ - # ensure we are capturing warnings that occur in the 'output' stage as these are warnings - # that occur at our solution point. So we clear existing warnings + # ensure we are capturing warnings that occur in the 'output' stage + # as these are warnings that occur at our solution point. + # So we clear existing warnings logging_model_handler.start_capturing() logging_model_handler.clear_logs() From 597b838369f1e5a02b64f24af2a128c3677dd6b2 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:50:45 +0100 Subject: [PATCH 03/10] use new solver output condition --- process/core/scan.py | 121 +--------------------------------- process/core/solver/solver.py | 12 ++-- 2 files changed, 8 insertions(+), 125 deletions(-) diff --git a/process/core/scan.py b/process/core/scan.py index 2720cdb987..29be92fc6f 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -331,7 +331,7 @@ def post_optimise(self, ifail: int): # Error code handler for VMCON if self.solver == "vmcon": - self.verror(ifail) + self.solver_handler.solver.verror(ifail) process_output.oblnkl(constants.NOUT) process_output.oblnkl(constants.IOTTY) else: @@ -742,125 +742,6 @@ def post_optimise(self, ifail: int): ), ) - @staticmethod - def verror(ifail: int): - """Routine to print out relevant messages in the case of an - unfeasible result from a VMCON (optimisation) run - - ifail : input integer : error flag - This routine prints out relevant messages in the case of - an unfeasible result from a VMCON (optimisation) run. - - Parameters - ---------- - ifail: int : - - """ - if ifail == SolverOutputCondition.USER_TERMINATED: - process_output.ocmmnt(constants.NOUT, "User-terminated execution of VMCON.") - process_output.ocmmnt(constants.IOTTY, "User-terminated execution of VMCON.") - elif ifail == SolverOutputCondition.IMPROPER_INPUT: - process_output.ocmmnt( - constants.NOUT, "Improper input parameters to the VMCON routine." - ) - process_output.ocmmnt(constants.NOUT, "PROCESS coding must be checked.") - - process_output.ocmmnt( - constants.IOTTY, "Improper input parameters to the VMCON routine." - ) - process_output.ocmmnt(constants.IOTTY, "PROCESS coding must be checked.") - elif ifail == SolverOutputCondition.MAX_ITERATIONS: - process_output.ocmmnt( - constants.NOUT, - "The maximum number of calls has been reached without solution.", - ) - process_output.ocmmnt( - constants.NOUT, - "The code may be stuck in a minimum in the residual space that is " - "significantly above zero.", - ) - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, "There is either no solution possible, or the code" - ) - process_output.ocmmnt( - constants.NOUT, "is failing to escape from a deep local minimum." - ) - process_output.ocmmnt( - constants.NOUT, - "Try changing the variables in IXC, or modify their initial values.", - ) - - process_output.ocmmnt( - constants.IOTTY, - "The maximum number of calls has been reached without solution.", - ) - process_output.ocmmnt( - constants.IOTTY, - "The code may be stuck in a minimum in the residual space that is " - "significantly above zero.", - ) - process_output.oblnkl(constants.NOUT) - process_output.oblnkl(constants.IOTTY) - process_output.ocmmnt( - constants.IOTTY, "There is either no solution possible, or the code" - ) - process_output.ocmmnt( - constants.IOTTY, "is failing to escape from a deep local minimum." - ) - process_output.ocmmnt( - constants.IOTTY, - "Try changing the variables in IXC, or modify their initial values.", - ) - elif ifail == SolverOutputCondition.MAX_LINE_SEARCHES: - process_output.ocmmnt( - constants.NOUT, "The line search required the maximum of 10 calls." - ) - process_output.ocmmnt( - constants.NOUT, "A feasible solution may be difficult to achieve." - ) - process_output.ocmmnt( - constants.NOUT, "Try changing or adding variables to IXC." - ) - - process_output.ocmmnt( - constants.IOTTY, "The line search required the maximum of 10 calls." - ) - process_output.ocmmnt( - constants.IOTTY, "A feasible solution may be difficult to achieve." - ) - process_output.ocmmnt( - constants.IOTTY, "Try changing or adding variables to IXC." - ) - elif ifail == SolverOutputCondition.NO_SOLUTION: - process_output.ocmmnt( - constants.NOUT, "The quadratic programming technique was unable to" - ) - process_output.ocmmnt(constants.NOUT, "find a feasible point.") - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, "Try changing or adding variables to IXC, or modify" - ) - process_output.ocmmnt( - constants.NOUT, - "their initial values (especially if only 1 optimisation", - ) - process_output.ocmmnt(constants.NOUT, "iteration was performed).") - - process_output.ocmmnt( - constants.IOTTY, "The quadratic programming technique was unable to" - ) - process_output.ocmmnt(constants.IOTTY, "find a feasible point.") - process_output.oblnkl(constants.IOTTY) - process_output.ocmmnt( - constants.IOTTY, "Try changing or adding variables to IXC, or modify" - ) - process_output.ocmmnt( - constants.IOTTY, - "their initial values (especially if only 1 optimisation", - ) - process_output.ocmmnt(constants.IOTTY, "iteration was performed).") - def scan_1d(self): """Run a 1-D scan.""" # initialise dict which will contain ifail values for each scan point diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index 846e525e3d..fdabebd3dd 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -304,12 +304,14 @@ def verror(self): """ strings = "\n".join( { - -1: ("User-terminated execution of VMCON.",), - 0: ( + SolverOutputCondition.USER_TERMINATED: ( + "User-terminated execution of VMCON.", + ), + SolverOutputCondition.IMPROPER_INPUT: ( "Improper input parameters to the VMCON routine.", "PROCESS coding must be checked.", ), - 2: ( + SolverOutputCondition.MAX_ITERATIONS: ( "The maximum number of calls has been reached without solution.", ( "The code may be stuck in a minimum in the residual space that" @@ -319,7 +321,7 @@ def verror(self): "is failing to escape from a deep local minimum.", "Try changing the variables in IXC, or modify their initial values.", ), - 3: ( + SolverOutputCondition.MAX_LINE_SEARCHES: ( "The line search required the maximum of 10 calls.", "A feasible solution may be difficult to achieve.", "Try changing or adding variables to IXC.", @@ -329,7 +331,7 @@ def verror(self): "Try changing the equations in ICC, or", "adding new variables to IXC.", ), - 5: ( + SolverOutputCondition.NO_SOLUTION: ( "The quadratic programming technique was unable to", "find a feasible point.\n", "Try changing or adding variables to IXC, or modify", From 345d1bdc687bd0c3d73fc54b1ff41a08b6169059 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 10 Aug 2026 14:59:48 +0100 Subject: [PATCH 04/10] avoid value error --- process/core/scan.py | 2 +- process/main.py | 13 ++----------- process/models/tfcoil/superconducting.py | 9 +++++++++ 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/process/core/scan.py b/process/core/scan.py index 29be92fc6f..77dec536fd 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -331,7 +331,7 @@ def post_optimise(self, ifail: int): # Error code handler for VMCON if self.solver == "vmcon": - self.solver_handler.solver.verror(ifail) + self.solver_handler.solver.verror() process_output.oblnkl(constants.NOUT) process_output.oblnkl(constants.IOTTY) else: diff --git a/process/main.py b/process/main.py index 2ca60408e7..84cabd31d6 100644 --- a/process/main.py +++ b/process/main.py @@ -853,7 +853,6 @@ def write(self, data, _outfile): physics and engineering model objects _outfile : int Fortran output unit identifier - """ # ensure we are capturing warnings that occur in the 'output' stage # as these are warnings that occur at our solution point. @@ -912,16 +911,8 @@ def write(self, data, _outfile): if data.tfcoil.i_tf_sup == TFConductorModel.SUPERCONDUCTING: tf_turn_type = SuperconductingTFTurnType( data.superconducting_tfcoil.i_tf_turn_type - ) - if tf_turn_type == SuperconductingTFTurnType.CABLE_IN_CONDUIT: - self.cicc_sctfcoil.output() - elif tf_turn_type == SuperconductingTFTurnType.CROSS_CONDUCTOR: - self.croco_sctfcoil.output() - else: - raise ValueError( - "Unsupported superconducting TF turn type: " - f"{data.superconducting_tfcoil.i_tf_turn_type}" - ) + ).abbreviation.lower() + getattr(self, f"{tf_turn_type}_sctfcoil").output() # Toroidal field coil aluminium model if data.tfcoil.i_tf_sup == TFConductorModel.HELIUM_COOLED_ALUMINIUM: diff --git a/process/models/tfcoil/superconducting.py b/process/models/tfcoil/superconducting.py index 55a4737f89..9c14f7b21e 100644 --- a/process/models/tfcoil/superconducting.py +++ b/process/models/tfcoil/superconducting.py @@ -63,6 +63,15 @@ def full_name(self): """The full name for this superconductor type.""" return self._full_name_ + @classmethod + def _missing_(cls, value): + try: + return cls[value] + except KeyError: + raise ValueError( + f"Unsupported superconducting TF turn type: {value}" + ) from None + class SuperconductingTFWPShapeType(IntEnum): """Enum for the type of TF coil WP shape, which determines the geometry of the From f2e5849da0e046d84565cf98dc908e4371462c9b Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 08:47:03 +0100 Subject: [PATCH 05/10] move solver handler stuff to solver handler --- process/core/scan.py | 453 +------------------------- process/core/solver/solver_handler.py | 131 +++++++- 2 files changed, 130 insertions(+), 454 deletions(-) diff --git a/process/core/scan.py b/process/core/scan.py index 77dec536fd..fc346314c4 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -10,20 +10,13 @@ from typing import TYPE_CHECKING import numpy as np -from tabulate import tabulate from process.core import constants, process_output from process.core.caller import write_output_files from process.core.exceptions import ProcessValueError from process.core.io.data_structure_dicts import get_dicts from process.core.log import logging_model_handler, show_errors -from process.core.solver import constraints from process.core.solver.solver_handler import SolverHandler -from process.data_structure.numerics import ( - FiguresOfMerit, - PROCESSRunMode, - SolverOutputCondition, -) from process.data_structure.scan_variables import IPNSCNS, NOUTVARS, ScanData from process.models.availability import AvailabilityModel @@ -293,454 +286,10 @@ def run_scan(self): def doopt(self): """Run the optimiser or solver.""" ifail = self.solver_handler.run() - self.post_optimise(ifail) - - return ifail - - def post_optimise(self, ifail: int): - """Called after calling the optimising equation solver from Python. - - ifail : input integer : error flag - - Parameters - ---------- - ifail: int : - - """ self.data.numerics.sqsumsq = ( sum(r**2 for r in self.data.numerics.rcm[: self.data.numerics.neqns]) ** 0.5 ) - - process_output.oheadr(constants.NOUT, "Numerics") - if self.solver == "fsolve": - process_output.ocmmnt( - constants.NOUT, "PROCESS has performed an fsolve (evaluation) run." - ) - else: - process_output.ocmmnt( - constants.NOUT, "PROCESS has performed a VMCON (optimisation) run." - ) - if ifail != SolverOutputCondition.CONVERGED: - process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) - process_output.oheadr( - constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" - ) - process_output.oblnkl(constants.IOTTY) - - logger.critical("Solver returns with ifail /= 1.\nifail = %s", ifail) - - # Error code handler for VMCON - if self.solver == "vmcon": - self.solver_handler.solver.verror() - process_output.oblnkl(constants.NOUT) - process_output.oblnkl(constants.IOTTY) - else: - # Solution found - if self.solver != "fsolve": - process_output.ocmmnt( - constants.NOUT, "and found a feasible set of parameters." - ) - process_output.oheadr( - constants.IOTTY, "PROCESS found a feasible solution" - ) - else: - process_output.ocmmnt( - constants.NOUT, "and found a consistent set of parameters." - ) - process_output.oheadr( - constants.IOTTY, "PROCESS found a consistent solution" - ) - process_output.oblnkl(constants.NOUT) - process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) - - if self.data.numerics.sqsumsq >= 1.0e-2: - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.NOUT, - "WARNING: Constraint residues are HIGH; consider re-running", - ) - process_output.ocmmnt( - constants.NOUT, - " with lower values of EPSVMC to confirm convergence...", - ) - process_output.ocmmnt( - constants.NOUT, - " (should be able to get down to about 1.0E-8 okay)", - ) - process_output.oblnkl(constants.NOUT) - process_output.ocmmnt( - constants.IOTTY, - "WARNING: Constraint residues are HIGH; consider re-running", - ) - process_output.ocmmnt( - constants.IOTTY, - " with lower values of EPSVMC to confirm convergence...", - ) - process_output.ocmmnt( - constants.IOTTY, - " (should be able to get down to about 1.0E-8 okay)", - ) - process_output.oblnkl(constants.IOTTY) - - logger.warning( - f"High final constraint residues. {self.data.numerics.sqsumsq=}" - ) - - process_output.ovarre( - constants.NOUT, - "Number of iteration variables", - "(nvar)", - self.data.numerics.nvar, - ) - process_output.ovarre( - constants.NOUT, - "Number of constraints (total)", - "(neqns+nineqns)", - self.data.numerics.neqns + self.data.numerics.nineqns, - ) - process_output.ovarre( - constants.NOUT, - "Optimisation switch", - "(ioptimz)", - self.data.numerics.ioptimz, - ) - process_output.ocmmnt( - constants.NOUT, - f" {PROCESSRunMode(self.data.numerics.ioptimz).description}", - ) - # Objective function output: none for fsolve - if self.solver != "fsolve": - process_output.ovarre( - constants.NOUT, - "Figure of merit switch", - "(minmax)", - self.data.numerics.minmax, - ) - - objf_name = f'"{FiguresOfMerit(abs(self.data.numerics.minmax)).description}"' - - self.data.numerics.objf_name = objf_name - - process_output.ovarre( - constants.NOUT, - "Objective function name", - "(objf_name)", - self.data.numerics.objf_name, - ) - process_output.ovarre( - constants.NOUT, - "Normalised objective function", - "(norm_objf)", - self.data.numerics.norm_objf, - "OP ", - ) - - process_output.ovarre( - constants.NOUT, - "Square root of the sum of squares of the constraint residuals", - "(sqsumsq)", - self.data.numerics.sqsumsq, - "OP ", - ) - if self.solver != "fsolve": - process_output.ovarre( - constants.NOUT, - "VMCON convergence parameter", - "(convergence_parameter)", - self.data.globals.convergence_parameter, - "OP ", - ) - process_output.ovarre( - constants.NOUT, - "Number of optimising solver iterations", - "(nviter)", - self.data.numerics.nviter, - "OP ", - ) - process_output.oblnkl(constants.NOUT) - - if self.solver == "fsolve": - if ifail == SolverOutputCondition.CONVERGED: - msg = "PROCESS has solved using fsolve." - else: - msg = "PROCESS failed to solve using fsolve." - process_output.write( - constants.NOUT, - f"{msg}\n", - ) - else: - if ifail == SolverOutputCondition.CONVERGED: - string1 = "PROCESS has successfully optimised" - else: - string1 = "PROCESS has failed to optimise" - - string2 = "minimise" if self.data.numerics.minmax > 0 else "maximise" - - process_output.write( - constants.NOUT, - f"{string1} the optimisation parameters to {string2} " - f"the objective function: {objf_name}\n", - ) - - written_warning = False - - # Output optimisation parameters - solution_vector_table = [] - for i in range(self.data.numerics.nvar): - self.data.numerics.xcs[i] = ( - self.data.numerics.xcm[i] * self.data.numerics.scafc[i] - ) - - name = self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1] - solution_vector_table.append([ - name, - self.data.numerics.xcs[i], - self.data.numerics.xcm[i], - ]) - - xminn = 1.01 * self.data.numerics.itv_scaled_lower_bounds[i] - xmaxx = 0.99 * self.data.numerics.itv_scaled_upper_bounds[i] - - # Write to output file if close to optimisation parameter bounds - if self.data.numerics.xcm[i] < xminn or self.data.numerics.xcm[i] > xmaxx: - if not written_warning: - written_warning = True - process_output.ocmmnt( - constants.NOUT, - ( - "Certain operating limits have been reached," - "\n as shown by the following optimisation parameters" - " that are" - "\n at or near to the edge of their prescribed range :\n" - ), - ) - - xcval = self.data.numerics.xcm[i] * self.data.numerics.scafc[i] - - if self.data.numerics.xcm[i] < xminn: - location, bound = "below", "lower" - bounds = self.data.numerics.itv_scaled_lower_bounds - else: - location, bound = "above", "upper" - bounds = self.data.numerics.itv_scaled_upper_bounds - process_output.write( - constants.NOUT, - f" {name:<30}= {xcval} is at or {location} its {bound} bound:" - f" {bounds[i] * self.data.numerics.scafc[i]}", - ) - - # Write optimisation parameters to mfile - process_output.ovarre( - constants.MFILE, - self.data.numerics.lablxc[self.data.numerics.ixc[i] - 1], - f"(itvar{i + 1:03d})", - self.data.numerics.xcs[i], - ) - - if self.data.numerics.boundu[i] == self.data.numerics.boundl[i]: - xnorm = 1.0 - else: - xnorm = min( - max( - ( - self.data.numerics.xcm[i] - - self.data.numerics.itv_scaled_lower_bounds[i] - ) - / ( - self.data.numerics.itv_scaled_upper_bounds[i] - - self.data.numerics.itv_scaled_lower_bounds[i] - ), - 0.0, - ), - 1.0, - ) - - process_output.ovarre( - constants.MFILE, - f"{name} (final value/initial value)", - f"(xcm{i + 1:03d})", - self.data.numerics.xcm[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} (range normalised)", - f"(nitvar{i + 1:03d})", - xnorm, - ) - process_output.ovarre( - constants.MFILE, - f"{name} (upper bound)", - f"(boundu{i + 1:03d})", - self.data.numerics.itv_scaled_upper_bounds[i] - * self.data.numerics.scafc[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} (lower bound)", - f"(boundl{i + 1:03d})", - self.data.numerics.itv_scaled_lower_bounds[i] - * self.data.numerics.scafc[i], - ) - - # Write optimisation parameter headings to output file - process_output.osubhd( - constants.NOUT, "The solution vector is comprised as follows :" - ) - process_output.write( - constants.NOUT, - tabulate( - solution_vector_table, - headers=["", "Final value", "Final / initial"], - numalign="left", - ), - ) - - process_output.osubhd( - constants.NOUT, - "The following equality constraint residues should be close to zero:", - ) - - con1, con2, err, _, lab = constraints.constraint_eqns( - self.data.numerics.neqns + self.data.numerics.nineqns, -1, self.data - ) - - # Write equality constraints to mfile - equality_constraint_table = [] - for i in range(self.data.numerics.neqns): - name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1] - - equality_constraint_table.append([ - name, - "=", - f"{con2[i]} {lab[i]}", - f"{err[i]} {lab[i]}", - con1[i], - ]) - process_output.ovarre( - constants.MFILE, - f"{name:<33} normalised residue", - f"(eq_con{self.data.numerics.icc[i]:03d})", - con1[i], - ) - - process_output.ovarre( - constants.MFILE, - f"{name:<33} residual", - f"(res_eq_con{self.data.numerics.icc[i]:03d})", - err[i], - ) - process_output.ovarre( - constants.MFILE, - f"{name} constraint value", - f"(val_eq_con{self.data.numerics.icc[i]:03d})", - con2[i], - ) - - process_output.ovarre( - constants.MFILE, - f"{name} units", - f"(eq_units_con{self.data.numerics.icc[i]:03d})", - f"'{lab[i]}'", - ) - - # Write equality constraints to output file - process_output.write( - constants.NOUT, - tabulate( - equality_constraint_table, - headers=[ - "", - "", - "Physical constraint", - "Constraint residue", - "Normalised residue", - ], - numalign="left", - ), - ) - - # Write inequality constraints - if self.data.numerics.nineqns > 0: - inequality_constraint_table = [] - # Inequalities not necessarily satisfied when evaluating - process_output.osubhd( - constants.NOUT, - "Negative inequality constraint (normalised) residuals " - "indicate a constraint is satisfied.", - ) - if self.solver == "fsolve": - process_output.osubhd( - constants.NOUT, - "This MFile was produced via an evaluation, not an optimisation, " - "and so the constraints might be violated.", - ) - - for i in range( - self.data.numerics.neqns, - self.data.numerics.neqns + self.data.numerics.nineqns, - ): - name = self.data.numerics.lablcc[self.data.numerics.icc[i] - 1] - constraint = constraints.ConstraintManager.evaluate_constraint( - int(self.data.numerics.icc[i]), self.data - ) - - inequality_constraint_table.append([ - name, - f"{constraint.constraint_value} {constraint.units}", - constraint.symbol, - f"{constraint.constraint_bound} {constraint.units}", - f"{constraint.residual} {constraint.units}", - f"{constraint.normalised_residual}", - ]) - process_output.ovarre( - constants.MFILE, - f"{name} normalised residue", - f"(ineq_con{self.data.numerics.icc[i]:03d})", - -constraint.normalised_residual, - ) - process_output.ovarre( - constants.MFILE, - f"{name} physical value", - f"(ineq_value_con{self.data.numerics.icc[i]:03d})", - constraint.constraint_value, - ) - - process_output.ovarre( - constants.MFILE, - f"{name} symbol", - f"(ineq_symbol_con{self.data.numerics.icc[i]:03d})", - f"'{constraint.symbol}'", - ) - - process_output.ovarre( - constants.MFILE, - f"{name} units", - f"(ineq_units_con{self.data.numerics.icc[i]:03d})", - f"'{constraint.units}'", - ) - - process_output.ovarre( - constants.MFILE, - f"{name} physical bound", - f"(ineq_bound_con{self.data.numerics.icc[i]:03d})", - constraint.constraint_bound, - ) - - process_output.write( - constants.NOUT, - tabulate( - inequality_constraint_table, - headers=[ - "", - "Physical constraint", - "", - "Physical constraint bound", - "Constraint residue", - "Normalised residue", - ], - numalign="left", - ), - ) + return ifail def scan_1d(self): """Run a 1-D scan.""" diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index a999cc821b..37cbe5e7f3 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -1,12 +1,21 @@ """Module containing solver handler routines""" +import logging + +from process.core import constants, process_output from process.core.solver.evaluators import Evaluators from process.core.solver.iteration_variables import ( load_iteration_variables, load_scaled_bounds, ) from process.core.solver.solver import get_solver -from process.data_structure.numerics import SolverOutputCondition +from process.data_structure.numerics import ( + FiguresOfMerit, + PROCESSRunMode, + SolverOutputCondition, +) + +logger = logging.getLogger(__name__) class SolverHandler: @@ -100,7 +109,6 @@ def run(self): ifail = self.solver.solve() self.output() - return ifail def output(self): @@ -113,3 +121,122 @@ def output(self): # than required, size self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf + + nums = self.data.numerics + + process_output.oheadr(constants.NOUT, "Numerics") + process_output.ocmmnt( + constants.NOUT, + f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'}" + " (optimisation) run.", + ) + ifail = self.solver.info + if ifail != SolverOutputCondition.CONVERGED: + process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) + process_output.oheadr( + constants.IOTTY, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" + ) + print() + + logger.critical("Solver returns with ifail /= 1. %s", ifail) + + if self.solver_name == "vmcon": + self.solver.verror() + + process_output.oblnkl(constants.NOUT) + print() + else: + # Solution found + descr = "consistent" if self.solver == "fsolve" else "feasible" + process_output.ocmmnt( + constants.NOUT, f"and found a {descr} set of parameters." + ) + process_output.oheadr(constants.IOTTY, f"PROCESS found a {descr} solution") + process_output.oblnkl(constants.NOUT) + process_output.ovarre(constants.NOUT, "Error flag", "(ifail)", ifail) + + if nums.sqsumsq >= 1.0e-2: + string = ( + "WARNING: Constraint residues are HIGH; consider re-running\n" + " with lower values of EPSVMC to confirm convergence...\n" + " (should be able to get down to about 1.0E-8 okay)\n" + ) + process_output.ocmmnt(constants.NOUT, ("\n" + string)) + print(string) + + logger.warning(f"High final constraint residues. {nums.sqsumsq=}") + + for d, var, v in ( + ("Number of iteration variables", "(nvar)", nums.nvar), + ( + "Number of constraints (total)", + "(neqns+nineqns)", + nums.neqns + nums.nineqns, + ), + ("Optimisation switch", "(ioptimz)", nums.ioptimz), + ): + process_output.ovarre(constants.NOUT, d, var, v) + + process_output.ocmmnt( + constants.NOUT, + f" {PROCESSRunMode(nums.ioptimz).description}", + ) + + # Objective function output: none for fsolve + if self.solver_name != "fsolve": + process_output.ovarre( + constants.NOUT, + "Figure of merit switch", + "(minmax)", + nums.minmax, + ) + + nums.objf_name = f'"{FiguresOfMerit(abs(nums.minmax)).description}"' + + for d, var, v, o in ( + ("Objective function name", "(objf_name)", nums.objf_name, ""), + ("Normalised objective function", "(norm_objf)", nums.norm_objf, "OP "), + ( + "VMCON convergence parameter", + "(convergence_parameter)", + self.data.globals.convergence_parameter, + "OP ", + ), + ( + "Number of optimising solver iterations", + "(nviter)", + nums.nviter, + "OP ", + ), + ( + "Square root of the sum of squares of the constraint residuals", + "(sqsumsq)", + nums.sqsumsq, + "OP ", + ), + ): + process_output.ovarre(constants.NOUT, d, var, v, o) + + process_output.oblnkl(constants.NOUT) + + if self.solver_name == "fsolve": + process_output.write( + constants.NOUT, + "PROCESS has solved using fsolve.\n" + if ifail == SolverOutputCondition.CONVERGED + else "PROCESS failed to solve using fsolve.\n", + ) + else: + process_output.write( + constants.NOUT, + ( + ( + "PROCESS has successfully optimised" + if ifail == SolverOutputCondition.CONVERGED + else "PROCESS has failed to optimise" + ) + + " the optimisation parameters to" + + ("minimise" if nums.minmax > 0 else "maximise") + + f" the objective function: {nums.objf_name}\n" + ), + ) From f9e27bbf2ad97705f1baadda314334a1cb3d1bef Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 08:55:10 +0100 Subject: [PATCH 06/10] use constraints output function --- process/core/scan.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/process/core/scan.py b/process/core/scan.py index fc346314c4..1737245c52 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -16,6 +16,7 @@ from process.core.exceptions import ProcessValueError from process.core.io.data_structure_dicts import get_dicts 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.scan_variables import IPNSCNS, NOUTVARS, ScanData from process.models.availability import AvailabilityModel @@ -289,6 +290,8 @@ def doopt(self): self.data.numerics.sqsumsq = ( sum(r**2 for r in self.data.numerics.rcm[: self.data.numerics.neqns]) ** 0.5 ) + constraints.constraints_output(self.data, self.solver) + return ifail def scan_1d(self): From d9194dfd172188434a76c9278c7443f9ff3ca6f2 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 08:58:28 +0100 Subject: [PATCH 07/10] missed doc --- process/core/solver/constraints.py | 1 + 1 file changed, 1 insertion(+) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 0ce2c67b01..c531a28be7 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -2014,6 +2014,7 @@ def constraint_eqns(m: int, ieqn: int, data: DataStructure): def constraints_output(data: DataStructure, solver_name: str): + """Output constraints information to file""" nums = data.numerics process_output.osubhd( From 60f505b04070bec2fabe7dae43bf16761fb4869e Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 09:18:37 +0100 Subject: [PATCH 08/10] missed table --- process/core/solver/solver_handler.py | 121 ++++++++++++++++++++++++-- 1 file changed, 113 insertions(+), 8 deletions(-) diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index 37cbe5e7f3..4feb6cff9b 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -2,6 +2,8 @@ import logging +from tabulate import tabulate + from process.core import constants, process_output from process.core.solver.evaluators import Evaluators from process.core.solver.iteration_variables import ( @@ -125,10 +127,12 @@ def output(self): nums = self.data.numerics process_output.oheadr(constants.NOUT, "Numerics") + s_type = ( + "fsolve (evaluation)" if self.solver == "fsolve" else "VMCON (optimisation)" + ) process_output.ocmmnt( constants.NOUT, - f"PROCESS has performed a {'fsolve' if self.solver == 'fsolve' else 'VMCON'}" - " (optimisation) run.", + f"PROCESS has performed a {s_type} run", ) ifail = self.solver.info if ifail != SolverOutputCondition.CONVERGED: @@ -208,15 +212,17 @@ def output(self): nums.nviter, "OP ", ), - ( - "Square root of the sum of squares of the constraint residuals", - "(sqsumsq)", - nums.sqsumsq, - "OP ", - ), ): process_output.ovarre(constants.NOUT, d, var, v, o) + process_output.ovarre( + constants.NOUT, + "Square root of the sum of squares of the constraint residuals", + "(sqsumsq)", + nums.sqsumsq, + "OP ", + ) + process_output.oblnkl(constants.NOUT) if self.solver_name == "fsolve": @@ -240,3 +246,102 @@ def output(self): + f" the objective function: {nums.objf_name}\n" ), ) + + written_warning = False + + # Output optimisation parameters + solution_vector_table = [] + for i in range(nums.nvar): + nums.xcs[i] = nums.xcm[i] * nums.scafc[i] + + name = nums.lablxc[nums.ixc[i] - 1] + solution_vector_table.append([name, nums.xcs[i], nums.xcm[i]]) + + xminn = 1.01 * nums.itv_scaled_lower_bounds[i] + xmaxx = 0.99 * nums.itv_scaled_upper_bounds[i] + + # Write to output file if close to optimisation parameter bounds + if nums.xcm[i] < xminn or nums.xcm[i] > xmaxx: + if not written_warning: + written_warning = True + process_output.ocmmnt( + constants.NOUT, + ( + "Certain operating limits have been reached," + "\n as shown by the following optimisation parameters" + " that are" + "\n at or near to the edge of their prescribed range :\n" + ), + ) + + xcval = nums.xcm[i] * nums.scafc[i] + + if nums.xcm[i] < xminn: + location, bound = "below", "lower" + bounds = nums.itv_scaled_lower_bounds + else: + location, bound = "above", "upper" + bounds = nums.itv_scaled_upper_bounds + process_output.write( + constants.NOUT, + f" {name:<30}= {xcval} is at or {location} its {bound} bound:" + f" {bounds[i] * nums.scafc[i]}", + ) + + if nums.boundu[i] == nums.boundl[i]: + xnorm = 1.0 + else: + xnorm = min( + max( + (nums.xcm[i] - nums.itv_scaled_lower_bounds[i]) + / ( + nums.itv_scaled_upper_bounds[i] + - nums.itv_scaled_lower_bounds[i] + ), + 0.0, + ), + 1.0, + ) + + # Write optimisation parameters to mfile + for d, var, v in ( + ( + nums.lablxc[nums.ixc[i] - 1], + f"(itvar{i + 1:03d})", + nums.xcs[i], + ), + ( + f"{name} (final value/initial value)", + f"(xcm{i + 1:03d})", + nums.xcm[i], + ), + ( + f"{name} (range normalised)", + f"(nitvar{i + 1:03d})", + xnorm, + ), + ( + f"{name} (upper bound)", + f"(boundu{i + 1:03d})", + nums.itv_scaled_upper_bounds[i] * nums.scafc[i], + ), + ( + f"{name} (lower bound)", + f"(boundl{i + 1:03d})", + nums.itv_scaled_lower_bounds[i] * nums.scafc[i], + ), + ): + process_output.ovarre(constants.MFILE, d, var, v) + + # Write optimisation parameter headings to output file + process_output.osubhd( + constants.NOUT, "The solution vector is comprised as follows :" + ) + process_output.write( + constants.NOUT, + tabulate( + solution_vector_table, + headers=["", "Final value", "Final / initial"], + numalign="left", + ), + ) From f58780774a6b808f82655efd68d5431e045452aa Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 09:35:55 +0100 Subject: [PATCH 09/10] reduce complexity --- process/core/solver/solver_handler.py | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index 4feb6cff9b..2be0347c6a 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -114,7 +114,7 @@ def run(self): return ifail def output(self): - """Store results back in Fortran self.data.numerics module. + """Store results back in self.data.numerics module. Objective function value, solution vector and constraints vector. """ @@ -124,6 +124,10 @@ def output(self): self.data.numerics.xcm[: self.solver.x.shape[0]] = self.solver.x self.data.numerics.rcm[: self.solver.conf.shape[0]] = self.solver.conf + self._numerics_output() + self._optimisation_parameters_output() + + def _numerics_output(self): nums = self.data.numerics process_output.oheadr(constants.NOUT, "Numerics") @@ -247,6 +251,9 @@ def output(self): ), ) + def _optimisation_parameters_output(self): + nums = self.data.numerics + written_warning = False # Output optimisation parameters @@ -305,21 +312,13 @@ def output(self): # Write optimisation parameters to mfile for d, var, v in ( - ( - nums.lablxc[nums.ixc[i] - 1], - f"(itvar{i + 1:03d})", - nums.xcs[i], - ), + (nums.lablxc[nums.ixc[i] - 1], f"(itvar{i + 1:03d})", nums.xcs[i]), ( f"{name} (final value/initial value)", f"(xcm{i + 1:03d})", nums.xcm[i], ), - ( - f"{name} (range normalised)", - f"(nitvar{i + 1:03d})", - xnorm, - ), + (f"{name} (range normalised)", f"(nitvar{i + 1:03d})", xnorm), ( f"{name} (upper bound)", f"(boundu{i + 1:03d})", From b6fd9ec7176ba0bf3bd5285d6b5caa365201e8df Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Tue, 11 Aug 2026 15:16:34 +0100 Subject: [PATCH 10/10] sqsumsq and convergence param fix --- process/core/scan.py | 3 --- process/core/solver/solver.py | 3 +-- process/core/solver/solver_handler.py | 2 ++ 3 files changed, 3 insertions(+), 5 deletions(-) diff --git a/process/core/scan.py b/process/core/scan.py index 1737245c52..d62a84a769 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -287,9 +287,6 @@ def run_scan(self): def doopt(self): """Run the optimiser or solver.""" ifail = self.solver_handler.run() - self.data.numerics.sqsumsq = ( - sum(r**2 for r in self.data.numerics.rcm[: self.data.numerics.neqns]) ** 0.5 - ) constraints.constraints_output(self.data, self.solver) return ifail diff --git a/process/core/solver/solver.py b/process/core/solver/solver.py index fdabebd3dd..269177e338 100644 --- a/process/core/solver/solver.py +++ b/process/core/solver/solver.py @@ -35,7 +35,6 @@ def __init__(self, *, data: DataStructure): self.data = data self.tolerance = self.data.numerics.epsvmc self.b: float | None = None - self.convergence_parameter: float | None = None self.maxcal = self.data.globals.maxcal def set_evaluators(self, evaluators: Evaluators): @@ -198,7 +197,7 @@ def solve(self) -> SolverOutputCondition: def _solver_callback(i: int, _result, _x, convergence_param: float): self.data.numerics.nviter = i + 1 - self.convergence_parameter = convergence_param + self.data.globals.convergence_parameter = convergence_param print( f"{i + 1} | Convergence Parameter: {convergence_param:.3E}", end="\r", diff --git a/process/core/solver/solver_handler.py b/process/core/solver/solver_handler.py index 2be0347c6a..4901d5983a 100644 --- a/process/core/solver/solver_handler.py +++ b/process/core/solver/solver_handler.py @@ -130,6 +130,8 @@ def output(self): def _numerics_output(self): nums = self.data.numerics + nums.sqsumsq = sum(r**2 for r in nums.rcm[: nums.neqns]) ** 0.5 + process_output.oheadr(constants.NOUT, "Numerics") s_type = ( "fsolve (evaluation)" if self.solver == "fsolve" else "VMCON (optimisation)"