diff --git a/process/scan.py b/process/scan.py index f7d4fdffc2..553079b505 100644 --- a/process/scan.py +++ b/process/scan.py @@ -1,8 +1,142 @@ +from dataclasses import astuple, dataclass + import numpy as np +from tabulate import tabulate +import process.process_output as process_output from process.caller import write_output_files -from process.fortran import error_handling, numerics, scan_module +from process.exceptions import ProcessValueError +from process.fortran import ( + build_variables, + constants, + constraint_variables, + constraints, + cost_variables, + cs_fatigue_variables, + current_drive_variables, + divertor_variables, + error_handling, + fwbs_variables, + global_variables, + heat_transport_variables, + impurity_radiation_module, + numerics, + pfcoil_variables, + physics_variables, + rebco_variables, + scan_module, + tfcoil_variables, +) from process.optimiser import Optimiser +from process.utilities.f2py_string_patch import ( + f2py_compatible_to_string, + string_to_f2py_compatible, +) + + +@dataclass +class ScanVariable: + variable_name: str + variable_description: str + + def __iter__(self): + return iter(astuple(self)) + + +SCAN_VARIABLES = { + 1: ScanVariable("aspect", "Aspect_ratio"), + 2: ScanVariable("hldivlim", "Div_heat_limit_(MW/m2)"), + 3: ScanVariable("pnetelin", "Net_electric_power_(MW)"), + 4: ScanVariable("hfact", "Confinement_H_factor"), + 5: ScanVariable("oacdcp", "TF_inboard_leg_J_(MA/m2)"), + 6: ScanVariable("walalw", "Allow._wall_load_(MW/m2)"), + 7: ScanVariable("beamfus0", "Beam_bkgrd_multiplier"), + 8: ScanVariable("fqval", "Big_Q_f-value"), + 9: ScanVariable("te", "Electron_temperature_keV"), + 10: ScanVariable("boundu(15)", "Volt-second_upper_bound"), + 11: ScanVariable("beta_norm_max", "Beta_coefficient"), + 12: ScanVariable("bootstrap_current_fraction_max", "Bootstrap_fraction"), + 13: ScanVariable("boundu(10)", "H_factor_upper_bound"), + 14: ScanVariable("fiooic", "TFC_Iop_/_Icrit_f-value"), + 15: ScanVariable("fjprot", "TFC_Jprot_limit_f-value"), + 16: ScanVariable("rmajor", "Plasma_major_radius_(m)"), + 17: ScanVariable("bmxlim", "Max_toroidal_field_(T)"), + 18: ScanVariable("gammax", "Maximum_CD_gamma"), + 19: ScanVariable("boundl(16)", "CS_thickness_lower_bound"), + 20: ScanVariable("t_burn_min", "Minimum_burn_time_(s)"), + 22: ScanVariable("cfactr", "Plant_availability_factor"), + 23: ScanVariable("boundu(72)", "Ip/Irod_upper_bound"), + 24: ScanVariable("powfmax", "Fusion_power_limit_(MW)"), + 25: ScanVariable("kappa", "Plasma_elongation"), + 26: ScanVariable("triang", "Plasma_triangularity"), + 27: ScanVariable("tbrmin", "Min_tritium_breed._ratio"), + 28: ScanVariable("bt", "Tor._field_on_axis_(T)"), + 29: ScanVariable("coreradius", "Core_radius"), + 31: ScanVariable( + "f_alpha_energy_confinement_min", "t_alpha_confinement/taueff_lower_limit" + ), + 32: ScanVariable("epsvmc", "VMCON error tolerance"), + 38: ScanVariable("boundu(129)", " Neon upper limit"), + 39: ScanVariable("boundu(131)", " Argon upper limit"), + 40: ScanVariable("boundu(135)", " Xenon upper limit"), + 41: ScanVariable("dr_blkt_outboard", "Outboard blanket thick."), + 42: ScanVariable("fimp(9)", "Argon fraction"), + 44: ScanVariable("sig_tf_case_max", "Allowable_stress_in_tf_coil_case_Tresca_(pa)"), + 45: ScanVariable("tmargmin_tf", "Minimum_allowable_temperature_margin"), + 46: ScanVariable("boundu(152)", "Max allowable fgwsep"), + 48: ScanVariable("n_pancake", "TF Coil - n_pancake"), + 49: ScanVariable("n_layer", "TF Coil - n_layer"), + 50: ScanVariable("fimp(13)", "Xenon fraction"), + 51: ScanVariable("ftar", "lower_divertor_power_fraction"), + 52: ScanVariable("rad_fraction_sol", "SoL radiation fraction"), + 53: ScanVariable("boundu(157)", "Max allowable fvssu"), + 54: ScanVariable("Bc2(0K)", "GL_NbTi Bc2(0K)"), + 55: ScanVariable("dr_shld_inboard", "Inboard neutronic shield"), + 56: ScanVariable("crypmw_max", "max allowable crypmw"), + 57: ScanVariable("boundl(2)", "bt minimum"), + 58: ScanVariable("dr_fw_plasma_gap_inboard", "Inboard FW-plasma sep gap"), + 59: ScanVariable("dr_fw_plasma_gap_outboard", "Outboard FW-plasma sep gap"), + 60: ScanVariable( + "sig_tf_wp_max", "Allowable_stress_in_tf_coil_conduit_Tresca_(pa)" + ), + 61: ScanVariable("copperaoh_m2_max", "Max CS coil current / copper area"), + 62: ScanVariable("coheof", "CS coil current density at EOF (A/m2)"), + 63: ScanVariable("dr_cs", "CS coil thickness (m)"), + 64: ScanVariable("ohhghf", "CS height (m)"), + 65: ScanVariable("n_cycle_min", "CS stress cycles min"), + 66: ScanVariable("oh_steel_frac", "CS steel fraction"), + 67: ScanVariable("t_crack_vertical", "Initial crack vertical size (m)"), + 68: ScanVariable( + "inlet_temp_liq", "Inlet Temperature Liquid Metal Breeder/Coolant (K)" + ), + 69: ScanVariable( + "outlet_temp_liq", "Outlet Temperature Liquid Metal Breeder/Coolant (K)" + ), + 70: ScanVariable( + "blpressure_liq", "Blanket liquid metal breeder/coolant pressure (Pa)" + ), + 71: ScanVariable( + "n_liq_recirc", "Selected number of liquid metal breeder recirculations per day" + ), + 72: ScanVariable( + "bz_channel_conduct_liq", + "Conductance of liquid metal breeder duct walls (A V-1 m-1)", + ), + 73: ScanVariable( + "pnuc_fw_ratio_dcll", "Ratio of FW nuclear power as fraction of total (FW+BB)" + ), + 74: ScanVariable( + "f_nuc_pow_bz_struct", + "Fraction of BZ power cooled by primary coolant for dual-coolant balnket", + ), + 75: ScanVariable("dx_fw_module", "dx_fw_module of first wall cooling channels (m)"), + 76: ScanVariable("etath", "Thermal conversion eff."), + 77: ScanVariable("startupratio", "Gyrotron redundancy"), + 78: ScanVariable("fkind", "Multiplier for Nth of a kind costs"), + 79: ScanVariable("etaech", "ECH wall plug to injector efficiency"), + 80: ScanVariable("fcoolcp", "Coolant fraction of TF"), + 81: ScanVariable("n_tf_turn", "Number of turns in TF"), +} class Scan: @@ -54,50 +188,480 @@ def doopt(self): return None ifail = self.optimiser.run() - scan_module.post_optimise(ifail) + self.post_optimise(ifail) return ifail - def scan_1d(self): - """Run a 1-D scan.""" - # f90wrap requires that arrays output from Fortran subroutines are - # defined in Python and passed in as an argument, in similar style to - # an intent(inout) argument. They are modified, but not returned. - # Initialise intent(out) array outvar - outvar = np.zeros( - (scan_module.noutvars, scan_module.ipnscns), dtype=np.float64, order="F" + def post_optimise(self, ifail: int): + """Called after calling the optimising equation solver from Python. + author: P J Knight, CCFE, Culham Science Centre + ifail : input integer : error flag + """ + numerics.sqsumsq = (numerics.rcm[: numerics.neqns] ** 2).sum() ** 0.5 + + error_handling.errors_on = True + + process_output.oheadr(constants.nout, "Numerics") + process_output.ocmmnt( + constants.nout, "PROCESS has performed a VMCON (optimisation) run." + ) + if ifail != 1: + process_output.ovarin(constants.nout, "VMCON error flag", "(ifail)", ifail) + process_output.oheadr( + constants.iotty, "PROCESS COULD NOT FIND A FEASIBLE SOLUTION" + ) + process_output.oblnkl(constants.iotty) + + error_handling.idiags[0] = ifail + error_handling.report_error(132) + + self.verror(ifail) + process_output.oblnkl(constants.nout) + process_output.oblnkl(constants.iotty) + else: + process_output.ocmmnt( + constants.nout, "and found a feasible set of parameters." + ) + process_output.oblnkl(constants.nout) + process_output.ovarin(constants.nout, "VMCON error flag", "(ifail)", ifail) + process_output.oheadr(constants.iotty, "PROCESS found a feasible solution") + + if 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) + + error_handling.fdiags[0] = numerics.sqsumsq + error_handling.report_error(134) + + process_output.ovarin( + constants.nout, "Number of iteration variables", "(nvar)", numerics.nvar + ) + process_output.ovarin( + constants.nout, + "Number of constraints (total)", + "(neqns+nineqns)", + numerics.neqns + numerics.nineqns, + ) + process_output.ovarin( + constants.nout, "Optimisation switch", "(ioptimz)", numerics.ioptimz + ) + process_output.ovarin( + constants.nout, "Figure of merit switch", "(minmax)", numerics.minmax + ) + + objf_name = string_to_f2py_compatible( + numerics.objf_name, + f'"{f2py_compatible_to_string(numerics.lablmm[abs(numerics.minmax) - 1])}"', ) + numerics.objf_name = objf_name + + process_output.ovarst( + constants.nout, "Objective function name", "(objf_name)", numerics.objf_name + ) + process_output.ovarre( + constants.nout, + "Normalised objective function", + "(norm_objf)", + numerics.norm_objf, + "OP ", + ) + process_output.ovarre( + constants.nout, + "Square root of the sum of squares of the constraint residuals", + "(sqsumsq)", + numerics.sqsumsq, + "OP ", + ) + process_output.ovarre( + constants.nout, + "VMCON convergence parameter", + "(convergence_parameter)", + global_variables.convergence_parameter, + "OP ", + ) + process_output.ovarin( + constants.nout, + "Number of VMCON iterations", + "(nviter)", + numerics.nviter, + "OP ", + ) + process_output.oblnkl(constants.nout) + + if ifail == 1: + string1 = "PROCESS has successfully optimised" + else: + string1 = "PROCESS has tried to optimise" + + string2 = "minimise" if numerics.minmax > 0 else "maximise" + + process_output.write( + constants.nout, + f"{string1} the iteration variables to {string2} the figure of merit: {objf_name}\n", + ) + + written_warning = False + + solution_vector_table = [] + for i in range(numerics.nvar): + numerics.xcs[i] = numerics.xcm[i] * numerics.scafc[i] + + name = f2py_compatible_to_string(numerics.lablxc[numerics.ixc[i] - 1]) + solution_vector_table.append([name, numerics.xcs[i], numerics.xcm[i]]) + + xminn = 1.01 * numerics.bondl[i] + xmaxx = 0.99 * numerics.bondu[i] + + if numerics.xcm[i] < xminn or 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 iteration variables that are" + "\n at or near to the edge of their prescribed range :\n" + ), + ) + + xcval = numerics.xcm[i] * numerics.scafc[i] + + if numerics.xcm[i] < xminn: + location, bound = "below", "lower" + else: + location, bound = "above", "upper" + process_output.write( + constants.nout, + f" {name:<30}= {xcval} is at or {location} its {bound} bound:" + f" {numerics.bondu[i] * numerics.scafc[i]}", + ) + + process_output.ovarre( + constants.mfile, + numerics.lablxc[numerics.ixc[i] - 1], + f"(itvar{i + 1:03d})", + numerics.xcs[i], + ) + + if numerics.boundu[i] == numerics.boundl[i]: + xnorm = 1.0 + else: + xnorm = min( + max( + (numerics.xcm[i] - numerics.bondl[i]) + / (numerics.bondu[i] - numerics.bondl[i]), + 0.0, + ), + 1.0, + ) + + process_output.ovarre( + constants.mfile, + f"{name} (final value/initial value)", + f"(xcm{i + 1:03d})", + numerics.xcm[i], + ) + process_output.ovarre( + constants.mfile, + f"{name} (range normalised)", + f"(nitvar{i + 1:03d})", + xnorm, + ) + + 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, sym, lab = constraints.constraint_eqns( + numerics.neqns + numerics.nineqns, -1 + ) + + equality_constraint_table = [] + for i in range(numerics.neqns): + name = f2py_compatible_to_string(numerics.lablcc[numerics.icc[i] - 1]) + + equality_constraint_table.append([ + name, + sym[i], + f"{con2[i]} {f2py_compatible_to_string(lab[i])}", + f"{err[i]} {f2py_compatible_to_string(lab[i])}", + con1[i], + ]) + process_output.ovarre( + constants.mfile, + f"{name:<33} normalised residue", + f"(eq_con{numerics.icc[i]:03d})", + con1[i], + ) + + process_output.write( + constants.nout, + tabulate( + equality_constraint_table, + headers=[ + "", + "", + "Physical constraint", + "Constraint residue", + "Normalised residue", + ], + numalign="left", + ), + ) + + if numerics.nineqns > 0: + inequality_constraint_table = [] + process_output.osubhd( + constants.nout, + "The following inequality constraint residues should be " + "greater than or approximately equal to zero :", + ) + + for i in range(numerics.neqns, numerics.neqns + numerics.nineqns): + name = f2py_compatible_to_string(numerics.lablcc[numerics.icc[i] - 1]) + inequality_constraint_table.append([ + name, + sym[i], + f"{con2[i]} {f2py_compatible_to_string(lab[i])}", + f"{err[i]} {f2py_compatible_to_string(lab[i])}", + ]) + process_output.ovarre( + constants.mfile, + f"{name} normalised residue", + f"(ineq_con{numerics.icc[i]:03d})", + numerics.rcm[i], + ) + + process_output.write( + constants.nout, + tabulate( + inequality_constraint_table, + headers=[ + "", + "", + "Physical constraint", + "Constraint residue", + ], + numalign="left", + ), + ) + + def verror(self, ifail: int): + """Routine to print out relevant messages in the case of an + unfeasible result from a VMCON (optimisation) run + author: P J Knight, CCFE, Culham Science Centre + ifail : input integer : error flag + This routine prints out relevant messages in the case of + an unfeasible result from a VMCON (optimisation) run. + """ + if ifail == -1: + process_output.ocmmnt(constants.nout, "User-terminated execution of VMCON.") + process_output.ocmmnt( + constants.iotty, "User-terminated execution of VMCON." + ) + elif ifail == 0: + process_output.ocmmnt( + constants.nout, "Improper input parameters to the VMCON routine." + ) + process_output.ocmmnt(constants.nout, "PROCESS coding must be checked.") + + process_output.ocmmnt( + constants.iotty, "Improper input parameters to the VMCON routine." + ) + process_output.ocmmnt(constants.iotty, "PROCESS coding must be checked.") + elif ifail == 2: + process_output.ocmmnt( + constants.nout, + "The maximum number of calls has been reached without solution.", + ) + process_output.ocmmnt( + constants.nout, + "The code may be stuck in a minimum in the residual space that is significantly above zero.", + ) + process_output.oblnkl(constants.nout) + process_output.ocmmnt( + constants.nout, "There is either no solution possible, or the code" + ) + process_output.ocmmnt( + constants.nout, "is failing to escape from a deep local minimum." + ) + process_output.ocmmnt( + constants.nout, + "Try changing the variables in IXC, or modify their initial values.", + ) + + process_output.ocmmnt( + constants.iotty, + "The maximum number of calls has been reached without solution.", + ) + process_output.ocmmnt( + constants.iotty, + "The code may be stuck in a minimum in the residual space that is significantly above zero.", + ) + process_output.oblnkl(constants.nout) + process_output.oblnkl(constants.iotty) + process_output.ocmmnt( + constants.iotty, "There is either no solution possible, or the code" + ) + process_output.ocmmnt( + constants.iotty, "is failing to escape from a deep local minimum." + ) + process_output.ocmmnt( + constants.iotty, + "Try changing the variables in IXC, or modify their initial values.", + ) + elif ifail == 3: + process_output.ocmmnt( + constants.nout, "The line search required the maximum of 10 calls." + ) + process_output.ocmmnt( + constants.nout, "A feasible solution may be difficult to achieve." + ) + process_output.ocmmnt( + constants.nout, "Try changing or adding variables to IXC." + ) + + process_output.ocmmnt( + constants.iotty, "The line search required the maximum of 10 calls." + ) + process_output.ocmmnt( + constants.iotty, "A feasible solution may be difficult to achieve." + ) + process_output.ocmmnt( + constants.iotty, "Try changing or adding variables to IXC." + ) + elif ifail == 4: + process_output.ocmmnt( + constants.nout, "An uphill search direction was found." + ) + process_output.ocmmnt( + constants.nout, "Try changing the equations in ICC, or" + ) + process_output.ocmmnt(constants.nout, "adding new variables to IXC.") + + process_output.ocmmnt( + constants.iotty, "An uphill search direction was found." + ) + process_output.ocmmnt( + constants.iotty, "Try changing the equations in ICC, or" + ) + process_output.ocmmnt(constants.iotty, "adding new variables to IXC.") + elif ifail == 5: + process_output.ocmmnt( + constants.nout, "The quadratic programming technique was unable to" + ) + process_output.ocmmnt(constants.nout, "find a feasible point.") + process_output.oblnkl(constants.nout) + process_output.ocmmnt( + constants.nout, "Try changing or adding variables to IXC, or modify" + ) + process_output.ocmmnt( + constants.nout, + "their initial values (especially if only 1 optimisation", + ) + process_output.ocmmnt(constants.nout, "iteration was performed).") + + process_output.ocmmnt( + constants.iotty, "The quadratic programming technique was unable to" + ) + process_output.ocmmnt(constants.iotty, "find a feasible point.") + process_output.oblnkl(constants.iotty) + process_output.ocmmnt( + constants.iotty, "Try changing or adding variables to IXC, or modify" + ) + process_output.ocmmnt( + constants.iotty, + "their initial values (especially if only 1 optimisation", + ) + process_output.ocmmnt(constants.iotty, "iteration was performed).") + elif ifail == 6: + process_output.ocmmnt( + constants.nout, "The quadratic programming technique was restricted" + ) + process_output.ocmmnt( + constants.nout, "by an artificial bound, or failed due to a singular" + ) + process_output.ocmmnt(constants.nout, "matrix.") + process_output.ocmmnt( + constants.nout, "Try changing the equations in ICC, or" + ) + process_output.ocmmnt(constants.nout, "adding new variables to IXC.") + + process_output.ocmmnt( + constants.iotty, "The quadratic programming technique was restricted" + ) + process_output.ocmmnt( + constants.iotty, "by an artificial bound, or failed due to a singular" + ) + process_output.ocmmnt(constants.iotty, "matrix.") + process_output.ocmmnt( + constants.iotty, "Try changing the equations in ICC, or" + ) + process_output.ocmmnt(constants.iotty, "adding new variables to IXC.") + + def scan_1d(self): + """Run a 1-D scan.""" # initialise dict which will contain ifail values for each scan point scan_1d_ifail_dict = {} for iscan in range(1, scan_module.isweep + 1): - scan_module.scan_1d_write_point_header(iscan) + self.scan_1d_write_point_header(iscan) ifail = self.doopt() scan_1d_ifail_dict[iscan] = ifail write_output_files(models=self.models, ifail=ifail) - # outvar is an intent(out) of scan_1d_store_output() - outvar = scan_module.scan_1d_store_output( - iscan, - ifail, - scan_module.noutvars, - scan_module.ipnscns, - ) error_handling.show_errors() error_handling.init_error_handling() # outvar now contains results - scan_module.scan_1d_write_plot(iscan, outvar) + self.scan_1d_write_plot() print( " ****************************************** Scan Convergence Summary ****************************************** \n" ) sweep_values = scan_module.sweep[: scan_module.isweep] - nsweep_var_name, _ = scan_module.scan_select( + nsweep_var_name, _ = self.scan_select( scan_module.nsweep, scan_module.sweep, scan_module.isweep ) converged_count = 0 - nsweep_var_name = nsweep_var_name.decode("utf-8") # offsets for aligning the converged/unconverged column max_sweep_value_length = len(str(np.max(sweep_values)).replace(".", "")) offsets = [ @@ -124,13 +688,7 @@ def scan_1d(self): def scan_2d(self): """Run a 2-D scan.""" # Initialise intent(out) arrays - outvar = np.zeros( - (scan_module.noutvars, scan_module.ipnscns), dtype=np.float64, order="F" - ) - sweep_1_vals = np.ndarray(scan_module.ipnscns, dtype=np.float64, order="F") - sweep_2_vals = np.ndarray(scan_module.ipnscns, dtype=np.float64, order="F") - - scan_module.scan_2d_init() + self.scan_2d_init() iscan = 1 # initialise array which will contain ifail values for each scan point @@ -139,42 +697,29 @@ def scan_2d(self): ) for iscan_1 in range(1, scan_module.isweep + 1): for iscan_2 in range(1, scan_module.isweep_2 + 1): - iscan_r = scan_module.scan_2d_write_point_header( - iscan, iscan_1, iscan_2 - ) + self.scan_2d_write_point_header(iscan, iscan_1, iscan_2) ifail = self.doopt() write_output_files(models=self.models, ifail=ifail) - outvar, sweep_1_vals, sweep_2_vals = scan_module.scan_2d_store_output( - ifail, - iscan_1, - iscan_r, - iscan, - scan_module.noutvars, - scan_module.ipnscns, - ) error_handling.show_errors() error_handling.init_error_handling() scan_2d_ifail_list[iscan_1][iscan_2] = ifail iscan = iscan + 1 - scan_module.scan_2d_write_plot(iscan, outvar, sweep_1_vals, sweep_2_vals) print( " ****************************************** Scan Convergence Summary ****************************************** \n" ) sweep_1_values = scan_module.sweep[: scan_module.isweep] sweep_2_values = scan_module.sweep_2[: scan_module.isweep_2] - nsweep_var_name, _ = scan_module.scan_select( + nsweep_var_name, _ = self.scan_select( scan_module.nsweep, scan_module.sweep, scan_module.isweep ) - nsweep_2_var_name, _ = scan_module.scan_select( + nsweep_2_var_name, _ = self.scan_select( scan_module.nsweep_2, scan_module.sweep_2, scan_module.isweep_2 ) converged_count = 0 scan_point = 1 - nsweep_var_name = nsweep_var_name.decode("utf-8") - nsweep_2_var_name = nsweep_2_var_name.decode("utf-8") # offsets for aligning the converged/unconverged column max_sweep1_value_length = len(str(np.max(sweep_1_values)).replace(".", "")) max_sweep2_value_length = len(str(np.max(sweep_2_values)).replace(".", "")) @@ -211,3 +756,275 @@ def scan_2d(self): converged_count / (scan_module.isweep * scan_module.isweep_2) * 100 ) print(f"\nConvergence Percentage: {converged_percentage:.2f}%") + + def scan_2d_init(self): + process_output.ovarin( + constants.mfile, + "Number of first variable scan points", + "(isweep)", + scan_module.isweep, + ) + process_output.ovarin( + constants.mfile, + "Number of second variable scan points", + "(isweep_2)", + scan_module.isweep_2, + ) + process_output.ovarin( + constants.mfile, + "Scanning first variable number", + "(nsweep)", + scan_module.nsweep, + ) + process_output.ovarin( + constants.mfile, + "Scanning second variable number", + "(nsweep_2)", + scan_module.nsweep_2, + ) + process_output.ovarin( + constants.mfile, + "Scanning second variable number", + "(nsweep_2)", + scan_module.nsweep_2, + ) + process_output.ovarin( + constants.mfile, + "Scanning second variable number", + "(nsweep_2)", + scan_module.nsweep_2, + ) + + def scan_1d_write_point_header(self, iscan: int): + global_variables.iscan_global = iscan + global_variables.vlabel, global_variables.xlabel = self.scan_select( + scan_module.nsweep, scan_module.sweep, iscan + ) + + process_output.oblnkl(constants.nout) + process_output.ostars(constants.nout, 110) + + process_output.write( + constants.nout, + f"***** Scan point {iscan} of {scan_module.isweep} : {f2py_compatible_to_string(global_variables.xlabel)}" + f", {f2py_compatible_to_string(global_variables.vlabel)} = {scan_module.sweep[iscan - 1]} " + "*****", + ) + process_output.ostars(constants.nout, 110) + process_output.oblnkl(constants.mfile) + process_output.ovarin(constants.mfile, "Scan point number", "(iscan)", iscan) + + print( + f"Starting scan point {iscan} of {scan_module.isweep} : " + f"{f2py_compatible_to_string(global_variables.xlabel)} , {f2py_compatible_to_string(global_variables.vlabel)}" + f" = {scan_module.sweep[iscan - 1]}" + ) + + def scan_2d_write_point_header(self, iscan, iscan_1, iscan_2): + iscan_r = scan_module.isweep_2 - iscan_2 + 1 if iscan_1 % 2 == 0 else iscan_2 + + # Makes iscan available globally (read-only) + global_variables.iscan_global = iscan + + global_variables.vlabel, global_variables.xlabel = self.scan_select( + scan_module.nsweep, scan_module.sweep, iscan_1 + ) + global_variables.vlabel_2, global_variables.xlabel_2 = self.scan_select( + scan_module.nsweep_2, scan_module.sweep_2, iscan_r + ) + + process_output.oblnkl(constants.nout) + process_output.ostars(constants.nout, 110) + + process_output.write( + constants.nout, + f"***** 2D Scan point {iscan} of {scan_module.isweep * scan_module.isweep_2} : " + f"{f2py_compatible_to_string(global_variables.vlabel)} = {scan_module.sweep[iscan - 1]} and" + f" {f2py_compatible_to_string(global_variables.vlabel_2)} = {scan_module.sweep_2[iscan_r]} " + "*****", + ) + process_output.ostars(constants.nout, 110) + process_output.oblnkl(constants.mfile) + process_output.ovarin(constants.mfile, "Scan point number", "(iscan)", iscan) + + print( + f"Starting scan point {iscan}: {f2py_compatible_to_string(global_variables.xlabel)}, " + f"{f2py_compatible_to_string(global_variables.vlabel)} = {scan_module.sweep[iscan - 1]}" + f" and {f2py_compatible_to_string(global_variables.xlabel_2)}, " + f"{f2py_compatible_to_string(global_variables.vlabel_2)} = {scan_module.sweep_2[iscan_r]} " + ) + + return iscan_r + + def scan_1d_write_plot(self): + if scan_module.first_call_1d: + process_output.ovarin( + constants.mfile, "Number of scan points", "(isweep)", scan_module.isweep + ) + process_output.ovarin( + constants.mfile, + "Scanning variable number", + "(nsweep)", + scan_module.nsweep, + ) + + scan_module.first_call_1d = False + + def scan_select(self, nwp, swp, iscn): + match nwp: + case 1: + physics_variables.aspect = swp[iscn - 1] + case 2: + divertor_variables.hldivlim = swp[iscn - 1] + case 3: + constraint_variables.pnetelin = swp[iscn - 1] + case 4: + physics_variables.hfact = swp[iscn - 1] + case 5: + tfcoil_variables.oacdcp = swp[iscn - 1] + case 6: + constraint_variables.walalw = swp[iscn - 1] + case 7: + physics_variables.beamfus0 = swp[iscn - 1] + case 8: + constraint_variables.fqval = swp[iscn - 1] + case 9: + physics_variables.te = swp[iscn - 1] + case 10: + numerics.boundu[14] = swp[iscn - 1] + case 11: + physics_variables.beta_norm_max = swp[iscn - 1] + case 12: + current_drive_variables.bootstrap_current_fraction_max = swp[iscn - 1] + case 13: + numerics.boundu[9] = swp[iscn - 1] + case 14: + constraint_variables.fiooic = swp[iscn - 1] + case 15: + constraint_variables.fjprot = swp[iscn - 1] + case 16: + physics_variables.rmajor = swp[iscn - 1] + case 17: + constraint_variables.bmxlim = swp[iscn - 1] + case 18: + constraint_variables.gammax = swp[iscn - 1] + case 19: + numerics.boundl[15] = swp[iscn - 1] + case 20: + constraint_variables.t_burn_min = swp[iscn - 1] + case 22: + if cost_variables.iavail == 1: + raise ProcessValueError("Do not scan cfactr if iavail=1") + cost_variables.cfactr = swp[iscn - 1] + case 23: + numerics.boundu[71] = swp[iscn - 1] + case 24: + constraint_variables.powfmax = swp[iscn - 1] + case 25: + physics_variables.kappa = swp[iscn - 1] + case 26: + physics_variables.triang = swp[iscn - 1] + case 27: + constraint_variables.tbrmin = swp[iscn - 1] + case 28: + physics_variables.bt = swp[iscn - 1] + case 29: + impurity_radiation_module.coreradius = swp[iscn - 1] + case 31: + constraint_variables.f_alpha_energy_confinement_min = swp[iscn - 1] + case 32: + numerics.epsvmc = swp[iscn - 1] + case 38: + numerics.boundu[128] = swp[iscn - 1] + case 39: + numerics.boundu[130] = swp[iscn - 1] + case 40: + numerics.boundu[134] = swp[iscn - 1] + case 41: + build_variables.dr_blkt_outboard = swp[iscn - 1] + case 42: + impurity_radiation_module.fimp[8] = swp[iscn - 1] + impurity_radiation_module.impurity_arr_frac[8] = ( + impurity_radiation_module.fimp[8] + ) + case 44: + tfcoil_variables.sig_tf_case_max = swp[iscn - 1] + case 45: + tfcoil_variables.tmargmin_tf = swp[iscn - 1] + case 46: + numerics.boundu[151] = swp[iscn - 1] + case 48: + tfcoil_variables.n_pancake = int(swp[iscn - 1]) + case 49: + tfcoil_variables.n_layer = int(swp[iscn - 1]) + case 50: + impurity_radiation_module.fimp[12] = swp[iscn - 1] + impurity_radiation_module.impurity_arr_frac[12] = ( + impurity_radiation_module.fimp[12] + ) + case 51: + physics_variables.ftar = swp[iscn - 1] + case 52: + physics_variables.rad_fraction_sol = swp[iscn - 1] + case 53: + numerics.boundu[156] = swp[iscn - 1] + case 54: + tfcoil_variables.b_crit_upper_nbti = swp[iscn - 1] + case 55: + build_variables.dr_shld_inboard = swp[iscn - 1] + case 56: + heat_transport_variables.crypmw_max = swp[iscn - 1] + case 57: + numerics.boundl[1] = swp[iscn - 1] + case 58: + build_variables.dr_fw_plasma_gap_inboard = swp[iscn - 1] + case 59: + build_variables.dr_fw_plasma_gap_outboard = swp[iscn - 1] + case 60: + tfcoil_variables.sig_tf_wp_max = swp[iscn - 1] + case 61: + rebco_variables.copperaoh_m2_max = swp[iscn - 1] + case 62: + pfcoil_variables.coheof = swp[iscn - 1] + case 63: + build_variables.dr_cs = swp[iscn - 1] + case 64: + pfcoil_variables.ohhghf = swp[iscn - 1] + case 65: + cs_fatigue_variables.n_cycle_min = swp[iscn - 1] + case 66: + pfcoil_variables.oh_steel_frac = swp[iscn - 1] + case 67: + cs_fatigue_variables.t_crack_vertical = swp[iscn - 1] + case 68: + fwbs_variables.inlet_temp_liq = swp[iscn - 1] + case 69: + fwbs_variables.outlet_temp_liq = swp[iscn - 1] + case 70: + fwbs_variables.blpressure_liq = swp[iscn - 1] + case 71: + fwbs_variables.n_liq_recirc = swp[iscn - 1] + case 72: + fwbs_variables.bz_channel_conduct_liq = swp[iscn - 1] + case 73: + fwbs_variables.pnuc_fw_ratio_dcll = swp[iscn - 1] + case 74: + fwbs_variables.f_nuc_pow_bz_struct = swp[iscn - 1] + case 75: + fwbs_variables.dx_fw_module = swp[iscn - 1] + case 76: + heat_transport_variables.etath = swp[iscn - 1] + case 77: + cost_variables.startupratio = swp[iscn - 1] + case 78: + cost_variables.fkind = swp[iscn - 1] + case 79: + current_drive_variables.etaech = swp[iscn - 1] + case 80: + tfcoil_variables.fcoolcp = swp[iscn - 1] + case 81: + tfcoil_variables.n_tf_turn = swp[iscn - 1] + case _: + raise ProcessValueError("Illegal scan variable number", nwp=nwp) + + return SCAN_VARIABLES[int(nwp)] diff --git a/scripts/create_dicts.py b/scripts/create_dicts.py index 6b0a183512..8f126a7e51 100644 --- a/scripts/create_dicts.py +++ b/scripts/create_dicts.py @@ -27,6 +27,8 @@ import numpy as np from python_dicts import get_python_variables +from process.scan import SCAN_VARIABLES + output_dict = {} # Dict of nested dicts e.g. output_dict['DICT_DESCRIPTIONS'] = # {descriptions_dict} @@ -630,53 +632,25 @@ def remove_comments(line): def dict_ixc2nsweep(): """Returns a dict mapping ixc_no to nsweep_no, if both exist for a - particular variable. Looks in scan.f90 for mapping from nsweep_no to - iteration variable name, and uses ixc_full to map variable name to ixc_no. - - Example of a fragment we are looking for: - case (1) - aspect = swp(iscn) - vlabel = 'aspect = ' ; xlabel = 'Aspect_ratio' + particular variable. Example dictionary entry: DICT_IXC2NSWEEP['1'] = '1' """ - ixc2nsweep = {} - file = SOURCEDIR + "/scan.f90" - # slice the file to get the switch statement relating to nsweep - lines = slice_file(file, r"select case \(nwp\)", r"case default") - - # remove extra lines that aren't case(#) or varname = sweep(iscan) lines - modlines = [] - for line in lines[1:-1]: - if "case" in line or "swp(iscn)" in line: - line = remove_comments(line).replace(" ", "") - modlines.append(line) + name_to_nsweep = { + var.variable_name: nsweep for nsweep, var in SCAN_VARIABLES.items() + } # create a dictionary that maps iteration variable names to ixc_no ixc_full = dict_ixc_full() - ixc_simple_rev = {} - for num, value in ixc_full.items(): - ixc_simple_rev[value["name"]] = num - - for i in range(len(modlines)): - # get the number from the case statement - match = re.match(r"case\((\d+)\)", modlines[i]) - if match: - num = match.group(1) - # if the case statement matched, get the variable name - # from the next line - match_2 = re.match(r"(.*?)=swp\(iscn\)", modlines[i + 1]) - if not match_2: - logging.warning("Error in dict_ixc2nsweep\n") - else: - name = match_2.group(1) - if name in ixc_simple_rev: - ixcnum = ixc_simple_rev[name] - ixc2nsweep[ixcnum] = num + name_to_ixc = { + name: ixc + for ixc, data in ixc_full.items() + if (name := data["name"]) in name_to_nsweep + } - return ixc2nsweep + return {ixc: name_to_nsweep[name] for name, ixc in name_to_ixc.items()} def dict_nsweep2ixc(): @@ -801,31 +775,7 @@ def dict_input_bounds(): def dict_nsweep2varname(): - # This function creates the nsweep2varname dictionary from the fortran code - # It maps the sweep variable number to its variable name - - di = {} - file = SOURCEDIR + "/scan.f90" - - # slice the file to get the switch statement relating to nsweep - lines = slice_file(file, r"select case \(nwp\)", r"case default") - - # remove extra lines that aren't case(#) or varname = sweep(iscan) lines - modlines = [] - for line in lines[1:-1]: - if "case" in line or "swp(iscn)" in line: - line = remove_comments(line).replace(" ", "") - modlines.append(line) - - for i in range(len(modlines) // 2): - line1 = modlines[i * 2] - no = line1.replace("case(", "") - no = no.replace(")", "") - line2 = modlines[i * 2 + 1] - varname = line2.replace("=swp(iscn)", "") - di[no] = varname - - return di + return {str(nsweep): var.variable_name for nsweep, var in SCAN_VARIABLES.items()} def dict_ixc_full(): diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index e75af906dc..0f3d0d130d 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -72,8 +72,6 @@ module process_input #ifdef unit_test public :: parse_input_file #endif -! public :: upper_case - integer, parameter :: maxlen = 2000 ! maximum line length character(len=maxlen) :: line ! current line of text from input file integer :: linelen, lineno ! current line length, line number diff --git a/source/fortran/scan.f90 b/source/fortran/scan.f90 index e2b9bcb8a0..957befaf45 100644 --- a/source/fortran/scan.f90 +++ b/source/fortran/scan.f90 @@ -146,1091 +146,4 @@ subroutine init_scan_module first_call_1d = .true. first_call_2d = .true. end subroutine init_scan_module - - subroutine scan_1d_write_point_header(iscan) - use global_variables, only: iscan_global, xlabel, vlabel - use constants, only: mfile, nout - use process_output_fortran, only: ovarin, ostars, oblnkl - implicit none - integer, intent(in) :: iscan - !! Scan point number - - ! Makes iscan available globally (read-only) - iscan_global = iscan - - call scan_select(nsweep, sweep, iscan, vlabel, xlabel) - - ! Write banner to output file - call oblnkl(nout) - call ostars(nout,width) - write(nout,10) ' ***** Scan point ', iscan,' of ',isweep,': & - ',trim(xlabel),', ',trim(vlabel),' = ',sweep(iscan),' *****' -10 format(a,i2,a,i2,5a,1pe10.3,a) - call ostars(nout,width) - - ! Write additional information to mfile - call oblnkl(mfile) - call ovarin(mfile,'Scan point number','(iscan)',iscan) - - ! Call the optimization routine VMCON at this scan point - write(*,20)'Starting scan point ',iscan,' of ',isweep,': ', trim(xlabel),', & - ',trim(vlabel),' = ',sweep(iscan) -20 format(a,i2,a,i2,a,4a,1pe10.3) - end subroutine scan_1d_write_point_header - - subroutine scan_1d_store_output(iscan, ifail, noutvars_, ipnscns_, outvar) - use constraint_variables, only: f_alpha_energy_confinement_min - use cost_variables, only: cdirt, coe, coeoam, coefuelt, c222, ireactor, & - capcost, coecap, c221 - use current_drive_variables, only: pheat, pinjmw, bootstrap_current_fraction, beam_energy, bigq - use divertor_variables, only: hldiv - use error_handling, only: errors_on - use heat_transport_variables, only: pgrossmw, pinjwp, pnetelmw - use impurity_radiation_module, only: fimp - use pfcoil_variables, only: whtpf - use pf_power_variables, only: srcktpm - use process_output_fortran, only: oblnkl - use numerics, only: sqsumsq - use tfcoil_variables, only: tfareain, wwp2, sig_tf_wp, tfcmw, tcpmax, oacdcp, & - tfcpmw, fcutfsu, acond, fcoolcp, rcool, whttf, ppump, vcool, wwp1, n_tf_coils, & - dr_tf_wp, b_crit_upper_nbti - use fwbs_variables, only: temp_fw_peak - use physics_variables, only: q95, aspect, p_plasma_rad_mw, dene, fusion_power, btot, tesep, & - pdivt, f_nd_alpha_electron, ten, beta_poloidal, hfac, teped, alpha_power_beams, q95_min, rmajor, wallmw, & - beta, beta_max, bt, plasma_current - use global_variables, only: verbose, maxcal, runtitle, run_tests - use constants, only: nout - implicit none - - integer, intent(in) :: iscan - integer, intent(in) :: ifail - ! outvar - integer, intent(in) :: noutvars_, ipnscns_ - real(dp), dimension(noutvars_,ipnscns_), intent(out) :: outvar - - ! Turn off error reporting (until next output) - errors_on = .false. - - ! Store values for MFILE.DAT output - outvar( 1,iscan) = dble(ifail) - outvar( 2,iscan) = sqsumsq - outvar( 3,iscan) = coe - outvar( 4,iscan) = coecap - outvar( 5,iscan) = coefuelt - outvar( 6,iscan) = coeoam - outvar( 7,iscan) = capcost - outvar( 8,iscan) = c221 + c222 - outvar( 9,iscan) = cdirt / 1.0D3 - outvar(10,iscan) = rmajor - outvar(11,iscan) = aspect - outvar(12,iscan) = 1.0D-6 * plasma_current - outvar(13,iscan) = bt - outvar(14,iscan) = btot - outvar(15,iscan) = q95 - outvar(16,iscan) = q95_min - outvar(17,iscan) = beta - outvar(18,iscan) = beta_max - outvar(19,iscan) = beta_poloidal / aspect - outvar(20,iscan) = ten/10.0D0 - outvar(21,iscan) = dene/1.0D20 - outvar(22,iscan) = hfac(6) - outvar(23,iscan) = hfac(7) - outvar(24,iscan) = fusion_power - outvar(25,iscan) = alpha_power_beams * 5.0D0 - outvar(26,iscan) = wallmw - outvar(27,iscan) = pinjmw - outvar(28,iscan) = pinjwp - outvar(29,iscan) = pheat - outvar(30,iscan) = pinjmw - pheat - outvar(31,iscan) = bigq - outvar(32,iscan) = bootstrap_current_fraction - outvar(33,iscan) = beam_energy/1.0D3 - outvar(34,iscan) = hldiv - outvar(35,iscan) = tfcmw - outvar(36,iscan) = whttf - outvar(37,iscan) = sig_tf_wp - outvar(38,iscan) = oacdcp/1.0D6 - outvar(39,iscan) = tcpmax - outvar(40,iscan) = tfcpmw - outvar(41,iscan) = fcoolcp - outvar(42,iscan) = rcool - outvar(43,iscan) = vcool - outvar(44,iscan) = ppump/1.0D6 - outvar(45,iscan) = 1.0D-3 * srcktpm - outvar(46,iscan) = whtpf - outvar(47,iscan) = pgrossmw - outvar(48,iscan) = pnetelmw - if (ireactor == 1) then - outvar(49,iscan) = (pgrossmw-pnetelmw) / pgrossmw - else - outvar(49,iscan) = 0.0D0 - end if - outvar(50,iscan) = pdivt/rmajor - !outvar(51,iscan) = fimpvar #OBSOLETE - outvar(51,iscan) = 0.0d0 - outvar(52,iscan) = p_plasma_rad_mw - outvar(53,iscan) = temp_fw_peak - outvar(54,iscan) = fcutfsu - outvar(55,iscan) = (wwp1+wwp2)*dr_tf_wp - outvar(56,iscan) = acond - outvar(57,iscan) = tfareain/n_tf_coils - outvar(58,iscan) = f_alpha_energy_confinement_min - outvar(66,iscan) = f_nd_alpha_electron - outvar(69,iscan) = fimp(1) - outvar(70,iscan) = fimp(2) - outvar(71,iscan) = fimp(3) - outvar(72,iscan) = fimp(4) - outvar(73,iscan) = fimp(5) - outvar(74,iscan) = fimp(6) - outvar(75,iscan) = fimp(7) - outvar(76,iscan) = fimp(8) - outvar(77,iscan) = fimp(9) - outvar(78,iscan) = fimp(10) - outvar(79,iscan) = fimp(11) - outvar(80,iscan) = fimp(12) - outvar(81,iscan) = fimp(13) - outvar(82,iscan) = fimp(14) - outvar(83,iscan) = teped - end subroutine scan_1d_store_output - - subroutine scan_1d_write_plot(iscan, outvar) - use global_variables, only: icase, xlabel - use constants, only: nplot, mfile - use process_output_fortran, only: ovarin - implicit none - - integer, intent(inout) :: iscan - real(dp), dimension(:,:), intent(in) :: outvar - - character(len=48) :: tlabel - integer :: ivar - character(len=25), dimension(noutvars), save :: plabel - - tlabel = icase - - ! Set up labels for plotting output - ! Use underscores instead of spaces - if (first_call_1d) then - plabel( 1) = 'Ifail____________________' - plabel( 2) = 'Sqsumsq__________________' - plabel( 3) = 'Electric_cost_(mil/kwh)__' - plabel( 4) = 'Capital_cost_(mil/kwh)___' - plabel( 5) = 'Fuel_cost_(mil/kwh)______' - plabel( 6) = 'Operations_cost_(mil/kwh)' - plabel( 7) = 'Capital_cost_(millions)__' - plabel( 8) = 'Core_costs_(millions)____' - plabel( 9) = 'Direct_cost_(billions)___' - plabel(10) = 'Major_Radius_(m)_________' - plabel(11) = 'Aspect_Ratio_____________' - plabel(12) = 'Plasma_Current_(MA)______' - plabel(13) = 'B_Toroidal_Axis_(T)______' - plabel(14) = 'B_total_on_axis_(T)______' - plabel(15) = 'Safety_Factor____________' - plabel(16) = 'q95_min_(zero_if_i_plasma_geometry=0)__' - plabel(17) = 'Beta_____________________' - plabel(18) = 'Beta_Limit_______________' - plabel(19) = 'Epsilon_Beta_Poloidal____' - plabel(20) = 'Dens.weight_Te_(10keV)___' - plabel(21) = 'Average_Dens_(10^20/m^3)_' - plabel(22) = 'H-fact_Iter_Power________' - plabel(23) = 'H-fact_Iter_Offset_______' - plabel(24) = 'Fusion_Power_(MW)________' - plabel(25) = 'nb_Fusion_Power_(MW)_____' - plabel(26) = 'Wall_Load_(MW/m^2)_______' - plabel(27) = 'Injection_Power_(MW)_____' - plabel(28) = 'Inject_Pwr_Wall_Plug_(MW)' - plabel(29) = 'Heating_Power_(MW)_______' - plabel(30) = 'Current_Drive_(MW)_______' - plabel(31) = 'Big_Q____________________' - plabel(32) = 'Bootstrap_Fraction_______' - plabel(33) = 'Neutral_Beam_Energy_(MeV)' - plabel(34) = 'Divertor_Heat_(MW/m^2)___' - plabel(35) = 'TF_coil_Power_(MW)_______' - plabel(36) = 'TF_coil_weight_(kg)______' - plabel(37) = 'vM_stress_in_TF_case_(Pa)' - plabel(38) = 'J_TF_inboard_leg_(MA/m^2)' - plabel(39) = 'Centrepost_max_T_(TART)__' - plabel(40) = 'Res_TF_inbrd_leg_Pwr_(MW)' - plabel(41) = 'Coolant_Fraction_Ctr.____' - plabel(42) = 'C/P_coolant_radius_(m)___' - plabel(43) = 'C/P_coolant_velocity(m/s)' - plabel(44) = 'C/P_pump_power_(MW)______' - plabel(45) = 'PF_coil_Power_(MW)_______' - plabel(46) = 'PF_coil_weight_(kg)______' - plabel(47) = 'Gross_Elect_Pwr_(MW)_____' - plabel(48) = 'Net_electric_Pwr_(MW)____' - plabel(49) = 'Recirculating_Fraction___' - plabel(50) = 'Psep/R___________________' - plabel(51) = '' !OBSOLETE - plabel(52) = 'Tot._radiation_power_(MW)' - plabel(53) = 'First_wall_peak_temp_(K)_' - plabel(54) = 'Cu_frac_TFC_conductor____' - plabel(55) = 'Winding_pack_area_TFC(m2)' - plabel(56) = 'Conductor_area_TFC_(m2)__' - plabel(57) = 'Area_TF_inboard_leg_(m2)_' - plabel(58) = 't_alpha_confinement/taueff_lower_limit__' - plabel(59) = 'Plasma_temp_at_sep_[keV]_' - plabel(60) = 'SOL_density_at_OMP_______' - plabel(61) = 'Power_through__separatrix' - plabel(62) = 'neomp/nesep______________' - plabel(63) = 'qtargettotal_____________' - plabel(64) = 'Total_pressure_at_target_' - plabel(65) = 'Temperature_at_target____' - plabel(66) = 'Helium_fraction__________' - plabel(67) = 'Momentum_loss_factor_____' - plabel(68) = 'totalpowerlost_[W]_______' - plabel(69) = 'H__concentration_________' - plabel(70) = 'He_concentration_________' - plabel(71) = 'Be_concentration_________' - plabel(72) = 'C__concentration_________' - plabel(73) = 'N__concentration_________' - plabel(74) = 'O__concentration_________' - plabel(75) = 'Ne_concentration_________' - plabel(76) = 'Si_concentration_________' - plabel(77) = 'Ar_concentration_________' - plabel(78) = 'Fe_concentration_________' - plabel(79) = 'Ni_concentration_________' - plabel(80) = 'Kr_concentration_________' - plabel(81) = 'Xe_concentration_________' - plabel(82) = 'W__concentration_________' - plabel(83) = 'teped____________________' - plabel(84) = 'Max_field_on_TF_coil_____' - call ovarin(mfile,'Number of scan points','(isweep)',isweep) - call ovarin(mfile,'Scanning variable number','(nsweep)',nsweep) - - first_call_1d = .false. - end if - - end subroutine scan_1d_write_plot - - subroutine scan_2d_init - !! Routine to call 2-D scan - !! author: J Morris, UKAEA, Culham Science Centre - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - use constants, only: mfile - use process_output_fortran, only: ovarin - implicit none - - ! Set up labels for plotting output - ! Use underscores instead of spaces - call ovarin(mfile,'Number of first variable scan points','(isweep)',isweep) - call ovarin(mfile,'Number of second variable scan points','(isweep_2)',isweep_2) - call ovarin(mfile,'Scanning first variable number','(nsweep)',nsweep) - call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) - call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) - call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) - end subroutine scan_2d_init - - subroutine scan_2d_write_point_header(iscan, iscan_1, iscan_2, iscan_R) - use process_output_fortran, only: oblnkl, ostars, ovarin - use global_variables, only: vlabel, vlabel_2, xlabel, xlabel_2, iscan_global - use constants, only: nout, mfile - implicit none - - integer, intent(in) :: iscan - integer, intent(in) :: iscan_1 - integer, intent(in) :: iscan_2 - integer, intent(out) :: iscan_R - - integer :: ifail - - if (mod(iscan_1,2)==0) then - iscan_R = isweep_2 - iscan_2 + 1 - else - iscan_R = iscan_2 - end if - ! Makes iscan available globally (read-only) - iscan_global = iscan - - call scan_select(nsweep, sweep, iscan_1, vlabel, xlabel) - call scan_select(nsweep_2, sweep_2, iscan_R, vlabel_2, xlabel_2) - - ! Write banner to output file - call oblnkl(nout) - call ostars(nout,width) - write(nout,10) iscan, isweep*isweep_2, trim(vlabel), & - sweep(iscan_1), trim(vlabel_2), sweep_2(iscan_R) -! 10 format(a, i2, a, i2, 5a, 1pe10.3, a) -10 format(' ***** 2D scan point ', i3, ' of ', i3, ' : ', a, ' = ', & - 1pe10.3, ' and ', a, ' = ', 1pe10.3, ' *****') - call ostars(nout,width) - - ! Write additional information to mfile - call oblnkl(mfile) - call ovarin(mfile,'Scan point number','(iscan)',iscan) - - ! Call the optimization routine VMCON at this scan point - write(*,20) iscan, trim(xlabel), trim(vlabel), sweep(iscan_1), & - trim(xlabel_2), trim(vlabel_2), sweep_2(iscan_R) - ! 20 format(a,i2,a,4a,1pe10.3) -20 format('Starting scan point ', i3, ': ', a, ', ', a, ' = ', & - 1pe10.3, ' and ', a, ', ', a, ' = ', 1pe10.3) - end subroutine scan_2d_write_point_header - - subroutine scan_2d_store_output(ifail, iscan_1, iscan_R, iscan, noutvars_, ipnscns_, outvar, & - sweep_1_vals, sweep_2_vals) - implicit none - - integer, intent(in) :: ifail - integer, intent(in) :: iscan_1 - integer, intent(in) :: iscan_R - integer, intent(in) :: iscan - integer, intent(in) :: noutvars_, ipnscns_ - ! Required for shape of intent(out) arrays - real(dp), dimension(noutvars_,ipnscns_), intent(out) :: outvar - real(dp), dimension(ipnscns_), intent(out) :: sweep_1_vals, sweep_2_vals - - call scan_1d_store_output(iscan, ifail, noutvars_, ipnscns_, outvar) - - sweep_1_vals(iscan) = sweep(iscan_1) - sweep_2_vals(iscan) = sweep_2(iscan_R) - end subroutine scan_2d_store_output - - subroutine scan_2d_write_plot(iscan, outvar, sweep_1_vals, sweep_2_vals) - use constants, only: nplot - use global_variables, only: icase, xlabel, xlabel_2 - implicit none - - integer, intent(inout) :: iscan - real(dp), dimension(:,:), intent(in) :: outvar - real(dp), dimension(:), intent(in) :: sweep_1_vals, sweep_2_vals - - integer :: ivar - character(len=48) :: tlabel - character(len=25), dimension(noutvars), save :: plabel - - plabel( 1) = 'Ifail____________________' - plabel( 2) = 'Sqsumsq__________________' - plabel( 3) = 'Electric_cost_(mil/kwh)__' - plabel( 4) = 'Capital_cost_(mil/kwh)___' - plabel( 5) = 'Fuel_cost_(mil/kwh)______' - plabel( 6) = 'Operations_cost_(mil/kwh)' - plabel( 7) = 'Capital_cost_(millions)__' - plabel( 8) = 'Core_costs_(millions)____' - plabel( 9) = 'Direct_cost_(billions)___' - plabel(10) = 'Major_Radius_(m)_________' - plabel(11) = 'Aspect_Ratio_____________' - plabel(12) = 'Plasma_Current_(MA)______' - plabel(13) = 'B_Toroidal_Axis_(T)______' - plabel(14) = 'B_total_on_axis_(T)______' - plabel(15) = 'Safety_Factor____________' - plabel(16) = 'q95_min_(zero_if_i_plasma_geometry=0)__' - plabel(17) = 'Beta_____________________' - plabel(18) = 'Beta_Limit_______________' - plabel(19) = 'Epsilon_Beta_Poloidal____' - plabel(20) = 'Dens.weight_Te_(10keV)___' - plabel(21) = 'Average_Dens_(10^20/m^3)_' - plabel(22) = 'H-fact_Iter_Power________' - plabel(23) = 'H-fact_Iter_Offset_______' - plabel(24) = 'Fusion_Power_(MW)________' - plabel(25) = 'nb_Fusion_Power_(MW)_____' - plabel(26) = 'Wall_Load_(MW/m^2)_______' - plabel(27) = 'Injection_Power_(MW)_____' - plabel(28) = 'Inject_Pwr_Wall_Plug_(MW)' - plabel(29) = 'Heating_Power_(MW)_______' - plabel(30) = 'Current_Drive_(MW)_______' - plabel(31) = 'Big_Q____________________' - plabel(32) = 'Bootstrap_Fraction_______' - plabel(33) = 'Neutral_Beam_Energy_(MeV)' - plabel(34) = 'Divertor_Heat_(MW/m^2)___' - plabel(35) = 'TF_coil_Power_(MW)_______' - plabel(36) = 'TF_coil_weight_(kg)______' - plabel(37) = 'vM_stress_in_TF_cond_(Pa)' - plabel(38) = 'J_TF_inboard_leg_(MA/m^2)' - plabel(39) = 'Centrepost_max_T_(TART)__' - plabel(40) = 'Res_TF_inbrd_leg_Pwr_(MW)' - plabel(41) = 'Coolant_Fraction_Ctr.____' - plabel(42) = 'C/P_coolant_radius_(m)___' - plabel(43) = 'C/P_coolant_velocity(m/s)' - plabel(44) = 'C/P_pump_power_(MW)______' - plabel(45) = 'PF_coil_Power_(MW)_______' - plabel(46) = 'PF_coil_weight_(kg)______' - plabel(47) = 'Gross_Elect_Pwr_(MW)_____' - plabel(48) = 'Net_electric_Pwr_(MW)____' - plabel(49) = 'Recirculating_Fraction___' - plabel(50) = 'Psep/R___________________' - plabel(51) = '' !OBSOLETE - plabel(52) = 'Tot._radiation_power_(MW)' - plabel(53) = 'First_wall_peak_temp_(K)_' - plabel(54) = 'Cu_frac_TFC_conductor____' - plabel(55) = 'Winding_pack_area_TFC(m2)' - plabel(56) = 'Conductor_area_TFC_(m2)__' - plabel(57) = 'Area_TF_inboard_leg_(m2)_' - plabel(58) = 't_alpha_confinement/taueff_lower_limit__' - plabel(59) = 'Plasma_temp_at_sep_[keV]_' - plabel(60) = 'SOL_density_at_OMP_______' - plabel(61) = 'Power_through__separatrix' - plabel(62) = 'neomp/nesep______________' - plabel(63) = 'qtargettotal_____________' - plabel(64) = 'Total_pressure_at_target_' - plabel(65) = 'Temperature_at_target____' - plabel(66) = 'Helium_fraction__________' - plabel(67) = 'Momentum_loss_factor_____' - plabel(68) = 'totalpowerlost_[W]_______' - plabel(69) = 'H__concentration_________' - plabel(70) = 'He_concentration_________' - plabel(71) = 'Be_concentration_________' - plabel(72) = 'C__concentration_________' - plabel(73) = 'N__concentration_________' - plabel(74) = 'O__concentration_________' - plabel(75) = 'Ne_concentration_________' - plabel(76) = 'Si_concentration_________' - plabel(77) = 'Ar_concentration_________' - plabel(78) = 'Fe_concentration_________' - plabel(79) = 'Ni_concentration_________' - plabel(80) = 'Kr_concentration_________' - plabel(81) = 'Xe_concentration_________' - plabel(82) = 'W__concentration_________' - plabel(83) = 'teped____________________' - plabel(84) = 'Max_field_on_TF_coil_____' - - tlabel = icase - - end subroutine scan_2d_write_plot - - subroutine scan_select(nwp, swp, iscn, vlab, xlab) - !! Routine to select first scan case - !! author: J Morris, UKAEA, Culham Science Centre - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use build_variables, only: dr_blkt_outboard, dr_shld_inboard, dr_fw_plasma_gap_inboard, dr_fw_plasma_gap_outboard, dr_cs - use constraint_variables, only: fiooic, walalw, bmxlim, fqval, f_alpha_energy_confinement_min, & - gammax, t_burn_min, tbrmin, fjprot, pnetelin, powfmax - use cost_variables, only: cfactr, iavail, fkind, startupratio - use current_drive_variables, only: bootstrap_current_fraction_max, etaech - use divertor_variables, only: hldivlim - use error_handling, only: idiags, report_error - use fwbs_variables, only: inlet_temp_liq, outlet_temp_liq, blpressure_liq, & - n_liq_recirc, bz_channel_conduct_liq, pnuc_fw_ratio_dcll, f_nuc_pow_bz_struct, dx_fw_module - use impurity_radiation_module, only: fimp, coreradius, impurity_arr_frac - use physics_variables, only: kappa, beta_norm_max, te, aspect, ftar, bt, & - rad_fraction_sol, triang, rmajor, beamfus0, hfact - use numerics, only: epsvmc, boundu, boundl - use tfcoil_variables, only: tmargmin_tf, sig_tf_case_max, n_pancake, oacdcp, & - n_layer, b_crit_upper_nbti, sig_tf_wp_max, fcoolcp, n_tf_turn - use heat_transport_variables, only: crypmw_max, etath - use rebco_variables, only: copperaoh_m2_max - use pfcoil_variables, only: coheof, ohhghf, oh_steel_frac - use CS_fatigue_variables, only: n_cycle_min, t_crack_vertical - implicit none - - ! Arguments - integer, intent(in) :: nwp, iscn - real(dp), intent(in), dimension(:) :: swp - character(len=25), intent(out) :: vlab, xlab - - select case (nwp) - ! Use underscores instead of spaces in xlabel - ! MDK Remove the "=" from vlabel, to make it easier to compare with - ! list of iteration variables - - case (1) - aspect = swp(iscn) - vlab = 'aspect' ; xlab = 'Aspect_ratio' - case (2) - hldivlim = swp(iscn) - vlab = 'hldivlim' ; xlab = 'Div_heat_limit_(MW/m2)' - case (3) - pnetelin = swp(iscn) - vlab = 'pnetelin' ; xlab = 'Net_electric_power_(MW)' - case (4) - hfact = swp(iscn) - vlab = 'hfact' ; xlab = 'Confinement_H_factor' - case (5) - oacdcp = swp(iscn) - vlab = 'oacdcp' ; xlab = 'TF_inboard_leg_J_(MA/m2)' - case (6) - walalw = swp(iscn) - vlab = 'walalw' ; xlab = 'Allow._wall_load_(MW/m2)' - case (7) - beamfus0 = swp(iscn) - vlab = 'beamfus0' ; xlab = 'Beam_bkgrd_multiplier' - case (8) - fqval = swp(iscn) - vlab = 'fqval' ; xlab = 'Big_Q_f-value' - case (9) - te = swp(iscn) - vlab = 'te' ; xlab = 'Electron_temperature_keV' - case (10) - boundu(15) = swp(iscn) - vlab = 'boundu(15)' ; xlab = 'Volt-second_upper_bound' - case (11) - beta_norm_max = swp(iscn) - vlab = 'beta_norm_max' ; xlab = 'Beta_coefficient' - case (12) - bootstrap_current_fraction_max = swp(iscn) - vlab = 'bootstrap_current_fraction_max' ; xlab = 'Bootstrap_fraction' - case (13) - boundu(10) = swp(iscn) - vlab = 'boundu(10)' ; xlab = 'H_factor_upper_bound' - case (14) - fiooic = swp(iscn) - vlab = 'fiooic' ; xlab = 'TFC_Iop_/_Icrit_f-value' - case (15) - fjprot = swp(iscn) - vlab = 'fjprot' ; xlab = 'TFC_Jprot_limit_f-value' - case (16) - rmajor = swp(iscn) - vlab = 'rmajor' ; xlab = 'Plasma_major_radius_(m)' - case (17) - bmxlim = swp(iscn) - vlab = 'bmxlim' ; xlab = 'Max_toroidal_field_(T)' - case (18) - gammax = swp(iscn) - vlab = 'gammax' ; xlab = 'Maximum_CD_gamma' - case (19) - boundl(16) = swp(iscn) - vlab = 'boundl(16)' ; xlab = 'CS_thickness_lower_bound' - case (20) - t_burn_min = swp(iscn) - vlab = 't_burn_min' ; xlab = 'Minimum_burn_time_(s)' - case (21) - ! sigpfalw = swp(iscn) - vlab = 'obsolete' ; xlab = 'obsolete' - case (22) - if (iavail == 1) call report_error(95) - cfactr = swp(iscn) - vlab = 'cfactr' ; xlab = 'Plant_availability_factor' - case (23) - boundu(72) = swp(iscn) - vlab = 'boundu(72)' ; xlab = 'Ip/Irod_upper_bound' - case (24) - powfmax = swp(iscn) - vlab = 'powfmax' ; xlab = 'Fusion_power_limit_(MW)' - case (25) - kappa = swp(iscn) - vlab = 'kappa' ; xlab = 'Plasma_elongation' - case (26) - triang = swp(iscn) - vlab = 'triang' ; xlab = 'Plasma_triangularity' - case (27) - tbrmin = swp(iscn) - vlab = 'tbrmin' ; xlab = 'Min_tritium_breed._ratio' - case (28) - bt = swp(iscn) - vlab = 'bt' ; xlab = 'Tor._field_on_axis_(T)' - case (29) - coreradius = swp(iscn) - vlab = 'coreradius' ; xlab = 'Core_radius' - case (30) - !fimpvar = swp(iscn) - vlab = 'OBSOLETE' ; xlab = 'OBSOLETE' - case (31) - f_alpha_energy_confinement_min = swp(iscn) - vlab = 'f_alpha_energy_confinement_min' ; xlab = 't_alpha_confinement/taueff_lower_limit' - case (32) - epsvmc = swp(iscn) - vlab = 'epsvmc' ; xlab = 'VMCON error tolerance' - case (33, 34, 35, 36, 37, 47) - write(*,*) 'Kallenbach model has been removed, remove the kallenbach scan variables' - stop 1 - case (38) - boundu(129) = swp(iscn) - vlab = 'boundu(129)' ; xlab = ' Neon upper limit' - case (39) - boundu(131) = swp(iscn) - vlab = 'boundu(131)' ; xlab = ' Argon upper limit' - case (40) - boundu(135) = swp(iscn) - vlab = 'boundu(135)' ; xlab = ' Xenon upper limit' - case (41) - dr_blkt_outboard = swp(iscn) - vlab = 'dr_blkt_outboard' ; xlab = 'Outboard blanket thick.' - case (42) - fimp(9) = swp(iscn) - impurity_arr_frac(9) = fimp(9) - vlab = 'fimp(9)' ; xlab = 'Argon fraction' - case (43) - ! rho_ecrh = swp(iscn) - vlab = 'obsolete' ; xlab = 'obsolete' - case (44) - sig_tf_case_max = swp(iscn) - vlab = 'sig_tf_case_max' ; xlab = 'Allowable_stress_in_tf_coil_case_Tresca_(pa)' - case (45) - tmargmin_tf = swp(iscn) - vlab = 'tmargmin_tf' ; xlab = 'Minimum_allowable_temperature_margin' - case (46) - boundu(152) = swp(iscn) - vlab = 'boundu(152)' ; xlab = 'Max allowable fgwsep' - case (48) - n_pancake = int(swp(iscn)) - vlab = 'n_pancake' ; xlab = 'TF Coil - n_pancake' - case (49) - n_layer = int(swp(iscn)) - vlab = 'n_layer' ; xlab = 'TF Coil - n_layer' - case (50) - fimp(13) = swp(iscn) - impurity_arr_frac(13) = fimp(13) - vlab = 'fimp(13)' ; xlab = 'Xenon fraction' - case (51) - ftar = swp(iscn) - vlab = 'ftar' ; xlab = 'lower_divertor_power_fraction' - case (52) - rad_fraction_sol = swp(iscn) - vlab = 'rad_fraction_sol' ; xlab = 'SoL radiation fraction' - case (53) - boundu(157) = swp(iscn) - vlab = 'boundu(157)' ; xlab = 'Max allowable fvssu' - case (54) - b_crit_upper_nbti = swp(iscn) - vlab = 'Bc2(0K)' ; xlab = 'GL_NbTi Bc2(0K)' - case(55) - dr_shld_inboard = swp(iscn) - vlab = 'dr_shld_inboard' ; xlab = 'Inboard neutronic shield' - case(56) - crypmw_max = swp(iscn) - vlab = 'crypmw_max' ; xlab = 'max allowable crypmw' - case(57) - boundl(2) = swp(iscn) - vlab = 'boundl(2)' ; xlab = 'bt minimum' - case(58) - dr_fw_plasma_gap_inboard = swp(iscn) - vlab = 'dr_fw_plasma_gap_inboard' ; xlab = 'Inboard FW-plasma sep gap' - case(59) - dr_fw_plasma_gap_outboard = swp(iscn) - vlab = 'dr_fw_plasma_gap_outboard' ; xlab = 'Outboard FW-plasma sep gap' - case (60) - sig_tf_wp_max = swp(iscn) - vlab = 'sig_tf_wp_max' ; xlab = 'Allowable_stress_in_tf_coil_conduit_Tresca_(pa)' - case (61) - copperaoh_m2_max = swp(iscn) - vlab = 'copperaoh_m2_max' ; xlab = 'Max CS coil current / copper area' - case (62) - coheof = swp(iscn) - vlab = 'coheof' ; xlab = 'CS coil current density at EOF (A/m2)' - case (63) - dr_cs = swp(iscn) - vlab = 'dr_cs' ; xlab = 'CS coil thickness (m)' - case (64) - ohhghf = swp(iscn) - vlab = 'ohhghf' ; xlab = 'CS height (m)' - case (65) - n_cycle_min = swp(iscn) - vlab = 'n_cycle_min' ; xlab = 'CS stress cycles min' - case (66) - oh_steel_frac = swp(iscn) - vlab = 'oh_steel_frac' ; xlab = 'CS steel fraction' - case (67) - t_crack_vertical = swp(iscn) - vlab = 't_crack_vertical' ; xlab = 'Initial crack vertical size (m)' - case (68) - inlet_temp_liq = swp(iscn) - vlab = 'inlet_temp_liq' ; xlab = 'Inlet Temperature Liquid Metal Breeder/Coolant (K)' - case (69) - outlet_temp_liq = swp(iscn) - vlab = 'outlet_temp_liq' ; xlab = 'Outlet Temperature Liquid Metal Breeder/Coolant (K)' - case(70) - blpressure_liq = swp(iscn) - vlab = 'blpressure_liq' ; xlab = 'Blanket liquid metal breeder/coolant pressure (Pa)' - case(71) - n_liq_recirc = swp(iscn) - vlab = 'n_liq_recirc' ; xlab = 'Selected number of liquid metal breeder recirculations per day' - case(72) - bz_channel_conduct_liq = swp(iscn) - vlab = 'bz_channel_conduct_liq' ; xlab = 'Conductance of liquid metal breeder duct walls (A V-1 m-1)' - case(73) - pnuc_fw_ratio_dcll = swp(iscn) - vlab = 'pnuc_fw_ratio_dcll' ; xlab = 'Ratio of FW nuclear power as fraction of total (FW+BB)' - case(74) - f_nuc_pow_bz_struct = swp(iscn) - vlab = 'f_nuc_pow_bz_struct' ; xlab = 'Fraction of BZ power cooled by primary coolant for dual-coolant balnket' - case(75) - dx_fw_module = swp(iscn) - vlab = 'dx_fw_module' ; xlab = 'pitch of first wall cooling channels (m)' - case (76) - etath = swp(iscn) - vlab = 'etath' ; xlab = 'Thermal conversion eff.' - case (77) - startupratio = swp(iscn) - vlab = 'startupratio' ; xlab = 'Gyrotron redundancy' - case (78) - fkind = swp(iscn) - vlab = 'fkind' ; xlab = 'Multiplier for Nth of a kind costs' - case (79) - etaech = swp(iscn) - vlab = 'etaech' ; xlab = 'ECH wall plug to injector efficiency' - case (80) - fcoolcp = swp(iscn) - vlab = 'fcoolcp' ; xlab = 'Coolant fraction of TF' - case (81) - n_tf_turn = swp(iscn) - vlab = 'n_tf_turn' ; xlab = 'Number of turns in TF' - case default - idiags(1) = nwp ; call report_error(96) - - end select - - end subroutine scan_select - - subroutine post_optimise(ifail) - !! Called after calling the optimising equation solver from Python. - !! author: P J Knight, CCFE, Culham Science Centre - !! ifail : input integer : error flag - !! - use constraints - use error_handling - use numerics - use process_output_fortran - use utilities, only:upper_case - ! for ipedestal = 2 option - use global_variables, only: convergence_parameter - use constants, only: iotty, nout, mfile - use physics_variables, only: ipedestal - use define_iteration_variables, only: boundxc, loadxc - implicit none - - ! Arguments - integer, intent(in) :: ifail - - ! Local variables - integer :: ii,inn,iflag - real(dp) :: summ,xcval,xmaxx,xminn,f,xnorm - real(dp), dimension(ipeqns) :: con1, con2, err - character(len=1), dimension(ipeqns) :: sym - character(len=10), dimension(ipeqns) :: lab - character(len=60) :: string1, string2 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Check on accuracy of solution by summing the - ! squares of the residuals of the equality constraints - summ = 0.0D0 - do ii = 1,neqns - summ = summ + rcm(ii)*rcm(ii) - end do - sqsumsq = sqrt(summ) - - ! Turn on error reporting - errors_on = .true. - - ! Print out information on solution - call oheadr(nout,'Numerics') - call ocmmnt(nout,'PROCESS has performed a VMCON (optimisation) run.') - if (ifail /= 1) then - !call ocmmnt(nout,'but could not find a feasible set of parameters.') - ! call oheadr(nout,'PROCESS COULD NOT FIND A FEASIBLE SOLUTION') - ! call ovarin(iotty,'VMCON error flag (ifail)','',ifail) - call ovarin(nout,'VMCON error flag','(ifail)',ifail) - call oheadr(iotty,'PROCESS COULD NOT FIND A FEASIBLE SOLUTION') - call oblnkl(iotty) - - idiags(1) = ifail ; call report_error(132) - - else - call ocmmnt(nout,'and found a feasible set of parameters.') - call oblnkl(nout) - call ovarin(nout,'VMCON error flag','(ifail)',ifail) - call oheadr(iotty,'PROCESS found a feasible solution') - end if - - !call oblnkl(nout) - - ! If necessary, write out a relevant error message - if (ifail /= 1) then - call verror(ifail) - call oblnkl(nout) - call oblnkl(iotty) - else - ! Show a warning if the constraints appear high even if allegedly converged - if (sqsumsq >= 1.0D-2) then - call oblnkl(nout) - call ocmmnt(nout,'WARNING: Constraint residues are HIGH; consider re-running') - call ocmmnt(nout,' with lower values of EPSVMC to confirm convergence...') - call ocmmnt(nout,' (should be able to get down to about 1.0E-8 okay)') - call oblnkl(nout) - call ocmmnt(iotty,'WARNING: Constraint residues are HIGH; consider re-running') - call ocmmnt(iotty,' with lower values of EPSVMC to confirm convergence...') - call ocmmnt(iotty,' (should be able to get down to about 1.0E-8 okay)') - call oblnkl(iotty) - - fdiags(1) = sqsumsq ; call report_error(134) - - end if - end if - - call ovarin(nout,'Number of iteration variables','(nvar)',nvar) - call ovarin(nout,'Number of constraints (total)','(neqns+nineqns)',neqns+nineqns) - call ovarin(nout,'Optimisation switch','(ioptimz)',ioptimz) - call ovarin(nout,'Figure of merit switch','(minmax)',minmax) -! if (ifail /= 1) then -! call ovarin(nout,'VMCON error flag','(ifail)',ifail) -! end if - - objf_name = '"'//trim(lablmm(abs(minmax)))//'"' - ! Quotes required for string parsing in MFILE - call ovarst(nout,'Objective function name','(objf_name)',objf_name) - call ovarre(nout,'Normalised objective function','(norm_objf)',norm_objf, 'OP ') - call ovarre(nout,'Square root of the sum of squares of the constraint residuals','(sqsumsq)',sqsumsq, 'OP ') - call ovarre(nout,'VMCON convergence parameter','(convergence_parameter)',convergence_parameter, 'OP ') - call ovarin(nout,'Number of VMCON iterations','(nviter)',nviter, 'OP ') - call oblnkl(nout) - - if (ifail == 1) then - string1 = 'PROCESS has successfully optimised the iteration variables' - else - string1 = 'PROCESS has tried to optimise the iteration variables' - end if - - if (minmax > 0) then - string2 = ' to minimise the figure of merit: ' - else - string2 = ' to maximise the figure of merit: ' - end if - - call upper_case(objf_name) - write(nout,10) trim(string1) // trim(string2), trim(objf_name) -10 format(a90, t92, a22) - - call oblnkl(nout) - - ! Check which variables are at bounds - iflag = 0 - do ii = 1,nvar - xminn = 1.01D0*bondl(ii) - xmaxx = 0.99D0*bondu(ii) - - if (xcm(ii) < xminn) then - if (iflag == 0) then - call ocmmnt(nout, & - 'Certain operating limits have been reached,') - call ocmmnt(nout, & - 'as shown by the following iteration variables that are') - call ocmmnt(nout, & - 'at or near to the edge of their prescribed range :') - call oblnkl(nout) - iflag = 1 - end if - xcval = xcm(ii)*scafc(ii) - !write(nout,30) ii,lablxc(ixc(ii)),xcval,bondl(ii)*scafc(ii) - write(nout,30) lablxc(ixc(ii)),xcval,bondl(ii)*scafc(ii) - end if - - if (xcm(ii) > xmaxx) then - if (iflag == 0) then - call ocmmnt(nout, & - 'Certain operating limits have been reached,') - call ocmmnt(nout, & - 'as shown by the following iteration variables that are') - call ocmmnt(nout, & - 'at or near to the edge of their prescribed range :') - call oblnkl(nout) - iflag = 1 - end if - xcval = xcm(ii)*scafc(ii) - write(nout,40) lablxc(ixc(ii)),xcval,bondu(ii)*scafc(ii) - end if - end do - -!30 format(t4,'Variable ',i3,' (',a9, & -! ',',1pe12.4,') is at or below its lower bound:',1pe12.4) -30 format(t4, a30, '=',1pe12.4,' is at or below its lower bound:',1pe12.4) -40 format(t4, a30, '=',1pe12.4,' is at or above its upper bound:',1pe12.4) -!40 format(t4,'Variable ',i3,' (',a9, & -! ',',1pe12.4,') is at or above its upper bound:',1pe12.4) - - ! Print out information on numerics - call osubhd(nout,'The solution vector is comprised as follows :') -! write(nout,50) -! Remove Lagrange multipliers as no-one understands them. -! MFILE not changed -!50 format(t47,'lower',t59,'upper') - - write(nout,60) -!60 format(t23,'final',t33,'fractional',t46,'Lagrange',t58,'Lagrange') -60 format(t43,'final',t55,'final /') - - - write(nout,70) -!70 format(t5,'i',t23,'value',t35,'change',t45,'multiplier', & -! t57,'multiplier') -70 format(t5,'i',t43,'value',t55,'initial') - - call oblnkl(nout) - - do inn = 1,nvar - xcs(inn) = xcm(inn)*scafc(inn) -! write(nout,80) inn,lablxc(ixc(inn)),xcs(inn),xcm(inn), & -! vlam(neqns+nineqns+inn), vlam(neqns+nineqns+1+inn+nvar) - write(nout,80) inn,lablxc(ixc(inn)),xcs(inn),xcm(inn) -!80 format(t2,i4,t8,a9,t19,4(1pe12.4)) -!80 format(t2,i4,t8,a30,t39,2(1pe12.4)) -80 format(t2,i4,t8,a30,t39,1pe12.4, t52, 0pf10.4) -! MDK The 0p is needed because of a bizarre "feature"/bug in fortran: -! the 1p in the previous format continues until changed. - call ovarre(mfile,lablxc(ixc(inn)),'(itvar'//int_to_string3(inn)//')',xcs(inn)) - - ! 'Range-normalised' iteration variable values for MFILE: - ! 0.0 (at lower bound) to 1.0 (at upper bound) - if (bondl(inn) == bondu(inn)) then - xnorm = 1.0D0 - else - xnorm = (xcm(inn) - bondl(inn)) / (bondu(inn) - bondl(inn)) - xnorm = max(xnorm, 0.0D0) - xnorm = min(xnorm, 1.0D0) - end if - ! Added ratio final/initial to MFILE - call ovarre(mfile,trim(lablxc(ixc(inn)))//' (final value/initial value)', & - '(xcm'//int_to_string3(inn)//')',xcm(inn)) - call ovarre(mfile,trim(lablxc(ixc(inn)))//' (range normalised)', & - '(nitvar'//int_to_string3(inn)//')',xnorm) - end do - - - call osubhd(nout, & - 'The following equality constraint residues should be close to zero :') - - call constraint_eqns(neqns+nineqns,-1,con1,con2,err,sym,lab) - write(nout,90) -90 format(t48,'physical',t73,'constraint',t100,'normalised') - write(nout,100) -100 format(t47,'constraint',t74,'residue',t101,'residue') - call oblnkl(nout) - do inn = 1,neqns - write(nout,110) inn,lablcc(icc(inn)),sym(inn),con2(inn), & - lab(inn),err(inn),lab(inn),con1(inn) - call ovarre(mfile,lablcc(icc(inn))//' normalised residue', & - '(eq_con'//int_to_string3(icc(inn))//')',con1(inn)) - end do -110 format(t2,i4,t8,a33,t46,a1,t47,1pe12.4,t60,a10,t71,1pe12.4,t84,a10,t98,1pe12.4) - - if (nineqns > 0) then - call osubhd(nout, & - 'The following inequality constraint residues should be greater than or approximately equal to zero :') - - do inn = neqns+1,neqns+nineqns - !write(nout,120) inn,lablcc(icc(inn)),rcm(inn),vlam(inn) - write(nout,110) inn,lablcc(icc(inn)),sym(inn),con2(inn), & - lab(inn), err(inn), lab(inn) - call ovarre(mfile,lablcc(icc(inn)),'(ineq_con'//int_to_string3(icc(inn))//')',rcm(inn)) - end do - end if - -! 120 format(t2,i4,t8,a33,t45,1pe12.4,1pe12.4) - -end subroutine post_optimise - -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - -subroutine verror(ifail) - - !! Routine to print out relevant messages in the case of an - !! unfeasible result from a VMCON (optimisation) run - !! author: P J Knight, CCFE, Culham Science Centre - !! ifail : input integer : error flag - !! This routine prints out relevant messages in the case of - !! an unfeasible result from a VMCON (optimisation) run. - !!

The messages are written to units NOUT and IOTTY, which are - !! by default the output file and screen, respectively. - !!

If IFAIL=1 then a feasible solution has been - !! found and therefore no error message is required. - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: nout, iotty - use process_output_fortran, only: ocmmnt, oblnkl - implicit none - - ! Arguments - integer, intent(in) :: ifail - - ! Local variables - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - select case (ifail) - - case (:-1) - call ocmmnt(nout, 'User-terminated execution of VMCON.') - call ocmmnt(iotty,'User-terminated execution of VMCON.') - - case (0) - call ocmmnt(nout, 'Improper input parameters to the VMCON routine.') - call ocmmnt(nout, 'PROCESS coding must be checked.') - - call ocmmnt(iotty,'Improper input parameters to the VMCON routine.') - call ocmmnt(iotty,'PROCESS coding must be checked.') - - case (1) - continue - - case (2) - call ocmmnt(nout,'The maximum number of calls has been reached without solution.') - call ocmmnt(nout,'The code may be stuck in a minimum in the residual space that is significantly above zero.') - call oblnkl(nout) - call ocmmnt(nout,'There is either no solution possible, or the code') - call ocmmnt(nout,'is failing to escape from a deep local minimum.') - call ocmmnt(nout,'Try changing the variables in IXC, or modify their initial values.') - - call ocmmnt(iotty,'The maximum number of calls has been reached without solution.') - call ocmmnt(iotty,'The code may be stuck in a minimum in the residual space that is significantly above zero.') - call oblnkl(nout) - call oblnkl(iotty) - call ocmmnt(iotty,'There is either no solution possible, or the code') - call ocmmnt(iotty,'is failing to escape from a deep local minimum.') - call ocmmnt(iotty,'Try changing the variables in IXC, or modify their initial values.') - - case (3) - call ocmmnt(nout,'The line search required the maximum of 10 calls.') - call ocmmnt(nout,'A feasible solution may be difficult to achieve.') - call ocmmnt(nout,'Try changing or adding variables to IXC.') - - call ocmmnt(iotty,'The line search required the maximum of 10 calls.') - call ocmmnt(iotty,'A feasible solution may be difficult to achieve.') - call ocmmnt(iotty,'Try changing or adding variables to IXC.') - - case (4) - call ocmmnt(nout,'An uphill search direction was found.') - call ocmmnt(nout,'Try changing the equations in ICC, or') - call ocmmnt(nout,'adding new variables to IXC.') - - call ocmmnt(iotty,'An uphill search direction was found.') - call ocmmnt(iotty,'Try changing the equations in ICC, or') - call ocmmnt(iotty,'adding new variables to IXC.') - - case (5) - call ocmmnt(nout, & - 'The quadratic programming technique was unable to') - call ocmmnt(nout,'find a feasible point.') - call oblnkl(nout) - call ocmmnt(nout,'Try changing or adding variables to IXC, or modify') - call ocmmnt(nout,'their initial values (especially if only 1 optimisation') - call ocmmnt(nout,'iteration was performed).') - - call ocmmnt(iotty, & - 'The quadratic programming technique was unable to') - call ocmmnt(iotty,'find a feasible point.') - call oblnkl(iotty) - call ocmmnt(iotty,'Try changing or adding variables to IXC, or modify') - call ocmmnt(iotty,'their initial values (especially if only 1 optimisation') - call ocmmnt(iotty,'iteration was performed).') - - case (6) - call ocmmnt(nout, & - 'The quadratic programming technique was restricted') - call ocmmnt(nout, & - 'by an artificial bound, or failed due to a singular') - call ocmmnt(nout,'matrix.') - call ocmmnt(nout,'Try changing the equations in ICC, or') - call ocmmnt(nout,'adding new variables to IXC.') - - call ocmmnt(iotty, & - 'The quadratic programming technique was restricted') - call ocmmnt(iotty, & - 'by an artificial bound, or failed due to a singular') - call ocmmnt(iotty,'matrix.') - call ocmmnt(iotty,'Try changing the equations in ICC, or') - call ocmmnt(iotty,'adding new variables to IXC.') - - case default - call ocmmnt(nout,'This value of IFAIL should not be possible...') - call ocmmnt(nout,'See source code for details.') - - call ocmmnt(iotty,'This value of IFAIL should not be possible...') - call ocmmnt(iotty,'See source code for details.') - - end select - -end subroutine verror - end module scan_module diff --git a/source/fortran/utilities.f90 b/source/fortran/utilities.f90 deleted file mode 100644 index b918c6f9fd..0000000000 --- a/source/fortran/utilities.f90 +++ /dev/null @@ -1,62 +0,0 @@ -module utilities - -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - implicit none - -contains -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - subroutine upper_case(string,start,finish) - - !! Routine that converts a (sub-)string to uppercase - !! author: P J Knight, CCFE, Culham Science Centre - !! string : input string : character string of interest - !! start : optional input integer : starting character for conversion - !! finish : optional input integer : final character for conversion - !! This routine converts the specified section of a string - !! to uppercase. By default, the whole string will be converted. - !! None - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - implicit none - - ! Arguments - - character(len=*), intent(inout) :: string - integer, optional, intent(in) :: start,finish - - ! Local variables - - character(len=1) :: letter - character(len=27), parameter :: upptab = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ_' - integer :: loop, i - - integer :: first, last - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - if (present(start)) then - first = start - else - first = 1 - end if - - if (present(finish)) then - last = finish - else - last = len(string) - end if - - if (first <= last) then - do loop = first,last - letter = string(loop:loop) - i = index('abcdefghijklmnopqrstuvwxyz_',letter) - if (i > 0) string(loop:loop) = upptab(i:i) - end do - end if - - end subroutine upper_case - -end module utilities