From 1037019c2b5e964740f75fb30ec6e9c24f85441e Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Thu, 14 Nov 2024 16:06:13 +0000 Subject: [PATCH 1/3] Convert PROCESS and module variable initialisation to Python --- CMakeLists.txt | 1 + process/init.py | 252 ++++++++ process/main.py | 9 +- scripts/vardes.py | 3 +- source/fortran/init_module.f90 | 910 ++++++++++++++++++++++++----- source/fortran/initial.f90 | 977 -------------------------------- tests/integration/test_vmcon.py | 9 +- tests/unit/conftest.py | 5 +- tests/unit/test_availability.py | 12 +- tests/unit/test_input.py | 3 +- 10 files changed, 1042 insertions(+), 1139 deletions(-) create mode 100644 process/init.py delete mode 100755 source/fortran/initial.f90 diff --git a/CMakeLists.txt b/CMakeLists.txt index a2be739107..000b25dde8 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -117,6 +117,7 @@ LIST(APPEND PROCESS_SRCS stellarator_variables.f90 stellarator.f90 stellarator_configuration.f90 + input.f90 ) PREPROCESS() diff --git a/process/init.py b/process/init.py new file mode 100644 index 0000000000..f58a51c793 --- /dev/null +++ b/process/init.py @@ -0,0 +1,252 @@ +import process.fortran as fortran + + +def init_process(): + """Routine that calls the initialisation routines + author: P J Knight, CCFE, Culham Science Centre + None + This routine calls the main initialisation routines that set + the default values for the global variables, reads in data from + the input file, and checks the run parameters for consistency. + """ + # Initialise error handling + fortran.error_handling.initialise_error_list() + + # Initialise the program variables + initialise_iterative_variables() + + # Initialise the Fortran file specifiers + # (creating and opening the files in the process) + fortran.init_module.open_files() + + # Input any desired new initial values + fortran.process_input.input() + + # Initialise the Stellarator + fortran.stellarator_module.stinit() + + # Check input data for errors/ambiguities + fortran.init_module.check() + + fortran.main_module.run_summary() + + +def init_all_module_vars(): + """Initialise all module variables + This is vital to ensure a 'clean' state of Process before a new run starts, + otherwise components of the previous run's state can persist into the new + run. This matters ever since Process is used as a shared library, rather + than a 'run-once' executable. + """ + fortran.numerics.init_numerics() + fortran.process_input.init_input() + fortran.buildings_variables.init_buildings_variables() + fortran.cost_variables.init_cost_variables() + fortran.divertor_variables.init_divertor_variables() + fortran.error_handling.init_error_handling() + fortran.fwbs_variables.init_fwbs_variables() + fortran.global_variables.init_global_variables() + fortran.ccfe_hcpb_module.init_ccfe_hcpb_module() + fortran.heat_transport_variables.init_heat_transport_variables() + fortran.ife_variables.init_ife_variables() + fortran.impurity_radiation_module.init_impurity_radiation_module() + fortran.pfcoil_module.init_pfcoil_module() + fortran.physics_module.init_physics_module() + fortran.physics_variables.init_physics_variables() + fortran.scan_module.init_scan_module() + fortran.sctfcoil_module.init_sctfcoil_module() + fortran.stellarator_module.init_stellarator_module() + fortran.stellarator_variables.init_stellarator_variables() + fortran.tfcoil_variables.init_tfcoil_variables() + fortran.times_variables.init_times_variables() + fortran.constants.init_constants() + fortran.current_drive_variables.init_current_drive_variables() + fortran.primary_pumping_variables.init_primary_pumping_variables() + fortran.pfcoil_variables.init_pfcoil_variables() + fortran.structure_variables.init_structure_variables() + fortran.vacuum_variables.init_vacuum_variables() + fortran.pf_power_variables.init_pf_power_variables() + fortran.build_variables.init_build_variables() + fortran.constraint_variables.init_constraint_variables() + fortran.pulse_variables.init_pulse_variables() + fortran.rebco_variables.init_rebco_variables() + fortran.reinke_variables.init_reinke_variables() + fortran.define_iteration_variables.init_define_iteration_variables() + fortran.reinke_module.init_reinke_module() + fortran.water_usage_variables.init_watuse_variables() + fortran.cs_fatigue_variables.init_cs_fatigue_variables() + fortran.blanket_library.init_blanket_library() + fortran.dcll_module.init_dcll_module() + + fortran.init_module.init_fortran_modules() + + +def initialise_iterative_variables(): + """Initialise each of the iteration variables""" + fortran.define_iteration_variables.init_itv_1() + fortran.define_iteration_variables.init_itv_2() + fortran.define_iteration_variables.init_itv_3() + fortran.define_iteration_variables.init_itv_4() + fortran.define_iteration_variables.init_itv_5() + fortran.define_iteration_variables.init_itv_6() + fortran.define_iteration_variables.init_itv_7() + fortran.define_iteration_variables.init_itv_8() + fortran.define_iteration_variables.init_itv_9() + fortran.define_iteration_variables.init_itv_10() + fortran.define_iteration_variables.init_itv_11() + fortran.define_iteration_variables.init_itv_12() + fortran.define_iteration_variables.init_itv_13() + fortran.define_iteration_variables.init_itv_14() + fortran.define_iteration_variables.init_itv_15() + fortran.define_iteration_variables.init_itv_16() + fortran.define_iteration_variables.init_itv_17() + fortran.define_iteration_variables.init_itv_18() + fortran.define_iteration_variables.init_itv_19() + fortran.define_iteration_variables.init_itv_20() + fortran.define_iteration_variables.init_itv_21() + + fortran.define_iteration_variables.init_itv_23() + + fortran.define_iteration_variables.init_itv_25() + fortran.define_iteration_variables.init_itv_26() + fortran.define_iteration_variables.init_itv_27() + fortran.define_iteration_variables.init_itv_28() + fortran.define_iteration_variables.init_itv_29() + fortran.define_iteration_variables.init_itv_30() + fortran.define_iteration_variables.init_itv_31() + fortran.define_iteration_variables.init_itv_32() + fortran.define_iteration_variables.init_itv_33() + fortran.define_iteration_variables.init_itv_34() + fortran.define_iteration_variables.init_itv_35() + fortran.define_iteration_variables.init_itv_36() + fortran.define_iteration_variables.init_itv_37() + fortran.define_iteration_variables.init_itv_38() + fortran.define_iteration_variables.init_itv_39() + fortran.define_iteration_variables.init_itv_40() + fortran.define_iteration_variables.init_itv_41() + fortran.define_iteration_variables.init_itv_42() + + fortran.define_iteration_variables.init_itv_44() + fortran.define_iteration_variables.init_itv_45() + fortran.define_iteration_variables.init_itv_46() + fortran.define_iteration_variables.init_itv_47() + fortran.define_iteration_variables.init_itv_48() + fortran.define_iteration_variables.init_itv_49() + fortran.define_iteration_variables.init_itv_50() + fortran.define_iteration_variables.init_itv_51() + + fortran.define_iteration_variables.init_itv_53() + fortran.define_iteration_variables.init_itv_54() + + fortran.define_iteration_variables.init_itv_56() + fortran.define_iteration_variables.init_itv_57() + fortran.define_iteration_variables.init_itv_58() + fortran.define_iteration_variables.init_itv_59() + fortran.define_iteration_variables.init_itv_60() + fortran.define_iteration_variables.init_itv_61() + fortran.define_iteration_variables.init_itv_62() + fortran.define_iteration_variables.init_itv_63() + fortran.define_iteration_variables.init_itv_64() + fortran.define_iteration_variables.init_itv_65() + fortran.define_iteration_variables.init_itv_66() + fortran.define_iteration_variables.init_itv_67() + fortran.define_iteration_variables.init_itv_68() + fortran.define_iteration_variables.init_itv_69() + fortran.define_iteration_variables.init_itv_70() + fortran.define_iteration_variables.init_itv_71() + fortran.define_iteration_variables.init_itv_72() + fortran.define_iteration_variables.init_itv_73() + fortran.define_iteration_variables.init_itv_74() + fortran.define_iteration_variables.init_itv_75() + + fortran.define_iteration_variables.init_itv_79() + + fortran.define_iteration_variables.init_itv_81() + fortran.define_iteration_variables.init_itv_82() + fortran.define_iteration_variables.init_itv_83() + + fortran.define_iteration_variables.init_itv_85() + fortran.define_iteration_variables.init_itv_86() + + fortran.define_iteration_variables.init_itv_89() + fortran.define_iteration_variables.init_itv_90() + fortran.define_iteration_variables.init_itv_91() + fortran.define_iteration_variables.init_itv_92() + fortran.define_iteration_variables.init_itv_93() + fortran.define_iteration_variables.init_itv_94() + fortran.define_iteration_variables.init_itv_95() + fortran.define_iteration_variables.init_itv_96() + fortran.define_iteration_variables.init_itv_97() + fortran.define_iteration_variables.init_itv_98() + + fortran.define_iteration_variables.init_itv_103() + fortran.define_iteration_variables.init_itv_104() + fortran.define_iteration_variables.init_itv_105() + fortran.define_iteration_variables.init_itv_106() + fortran.define_iteration_variables.init_itv_107() + fortran.define_iteration_variables.init_itv_108() + fortran.define_iteration_variables.init_itv_109() + fortran.define_iteration_variables.init_itv_110() + fortran.define_iteration_variables.init_itv_111() + fortran.define_iteration_variables.init_itv_112() + fortran.define_iteration_variables.init_itv_113() + fortran.define_iteration_variables.init_itv_114() + fortran.define_iteration_variables.init_itv_115() + fortran.define_iteration_variables.init_itv_116() + fortran.define_iteration_variables.init_itv_117() + fortran.define_iteration_variables.init_itv_118() + fortran.define_iteration_variables.init_itv_119() + fortran.define_iteration_variables.init_itv_120() + fortran.define_iteration_variables.init_itv_121() + fortran.define_iteration_variables.init_itv_122() + fortran.define_iteration_variables.init_itv_123() + fortran.define_iteration_variables.init_itv_124() + fortran.define_iteration_variables.init_itv_125() + fortran.define_iteration_variables.init_itv_126() + fortran.define_iteration_variables.init_itv_127() + fortran.define_iteration_variables.init_itv_128() + fortran.define_iteration_variables.init_itv_129() + fortran.define_iteration_variables.init_itv_130() + fortran.define_iteration_variables.init_itv_131() + fortran.define_iteration_variables.init_itv_132() + fortran.define_iteration_variables.init_itv_133() + fortran.define_iteration_variables.init_itv_134() + fortran.define_iteration_variables.init_itv_135() + fortran.define_iteration_variables.init_itv_136() + fortran.define_iteration_variables.init_itv_137() + fortran.define_iteration_variables.init_itv_138() + fortran.define_iteration_variables.init_itv_139() + fortran.define_iteration_variables.init_itv_140() + fortran.define_iteration_variables.init_itv_141() + fortran.define_iteration_variables.init_itv_142() + fortran.define_iteration_variables.init_itv_143() + fortran.define_iteration_variables.init_itv_144() + fortran.define_iteration_variables.init_itv_145() + fortran.define_iteration_variables.init_itv_146() + fortran.define_iteration_variables.init_itv_147() + fortran.define_iteration_variables.init_itv_148() + fortran.define_iteration_variables.init_itv_149() + fortran.define_iteration_variables.init_itv_152() + fortran.define_iteration_variables.init_itv_153() + fortran.define_iteration_variables.init_itv_154() + fortran.define_iteration_variables.init_itv_155() + fortran.define_iteration_variables.init_itv_156() + fortran.define_iteration_variables.init_itv_157() + fortran.define_iteration_variables.init_itv_158() + fortran.define_iteration_variables.init_itv_159() + fortran.define_iteration_variables.init_itv_160() + fortran.define_iteration_variables.init_itv_161() + fortran.define_iteration_variables.init_itv_162() + fortran.define_iteration_variables.init_itv_163() + fortran.define_iteration_variables.init_itv_164() + fortran.define_iteration_variables.init_itv_165() + fortran.define_iteration_variables.init_itv_166() + fortran.define_iteration_variables.init_itv_167() + fortran.define_iteration_variables.init_itv_168() + fortran.define_iteration_variables.init_itv_169() + fortran.define_iteration_variables.init_itv_170() + fortran.define_iteration_variables.init_itv_171() + fortran.define_iteration_variables.init_itv_172() + fortran.define_iteration_variables.init_itv_173() + fortran.define_iteration_variables.init_itv_174() + fortran.define_iteration_variables.init_itv_175() diff --git a/process/main.py b/process/main.py index be7052a465..94cabae036 100644 --- a/process/main.py +++ b/process/main.py @@ -73,6 +73,7 @@ from process.current_drive import CurrentDrive from process.impurity_radiation import initialise_imprad from process.caller import write_output_files +import process.init as init import process @@ -298,8 +299,8 @@ def run(self): config = RunProcessConfig(self.config_file) config.setup() - fortran.init_module.init_all_module_vars() - fortran.init_module.init() + init.init_all_module_vars() + init.init_process() neqns, itervars = get_neqns_itervars() lbs, ubs = get_variable_range(itervars, config.factor) @@ -399,7 +400,7 @@ def init_module_vars(): This "resets" all module variables to their initialised values, so each new run doesn't have any side-effects from previous runs. """ - fortran.init_module.init_all_module_vars() + init.init_all_module_vars() def set_filenames(self): """Validate the input filename and create other filenames from it.""" @@ -455,7 +456,7 @@ def initialise(): """Run the init module to call all initialisation routines.""" initialise_imprad() # Reads in input file - fortran.init_module.init() + init.init_process() # Order optimisation parameters (arbitrary order in input file) # Ensures consistency and makes output comparisons more straightforward diff --git a/scripts/vardes.py b/scripts/vardes.py index ab3cdf8831..baee244c33 100644 --- a/scripts/vardes.py +++ b/scripts/vardes.py @@ -7,6 +7,7 @@ import numpy as np import jinja2 +from process.init import init_all_module_vars from process import fortran @@ -81,7 +82,7 @@ def get_input_output_variables(variables: List[FortranVariable]): except AttributeError: continue - fortran.init_module.init_all_module_vars() + init_all_module_vars() for var in variables: current_values_entry = f"{var.module}.{var.name}" diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90 index 1b49da3b01..ee4d42462a 100644 --- a/source/fortran/init_module.f90 +++ b/source/fortran/init_module.f90 @@ -1,133 +1,31 @@ -#ifndef INSTALLDIR -#error INSTALLDIR not defined! -#endif - module init_module +#ifndef dp +use, intrinsic :: iso_fortran_env, only: dp=>real64 +#endif + implicit none contains - subroutine init_all_module_vars - !! Initialise all module variables - !! This is vital to ensure a 'clean' state of Process before a new run starts, - !! otherwise components of the previous run's state can persist into the new - !! run. This matters ever since Process is used as a shared library, rather - !! than a 'run-once' executable. - use numerics, only: init_numerics - use process_input, only: init_input - use buildings_variables, only: init_buildings_variables - use cost_variables, only: init_cost_variables - use divertor_variables, only: init_divertor_variables - use error_handling, only: init_error_handling + subroutine init_fortran_modules + !! Temporary routine to call initialisation routines for Fortran modules + !! that are not wrapped by f2py and thus cannot be called from Python. + use fson_library, only: init_fson_library - use fwbs_variables, only: init_fwbs_variables - use global_variables, only: init_global_variables - use ccfe_hcpb_module, only: init_ccfe_hcpb_module - use heat_transport_variables, only: init_heat_transport_variables - use ife_variables, only: init_ife_variables - use impurity_radiation_module, only: init_impurity_radiation_module - use pfcoil_module, only: init_pfcoil_module - use physics_module, only: init_physics_module - use physics_variables, only: init_physics_variables - use scan_module, only: init_scan_module - use sctfcoil_module, only: init_sctfcoil_module - use stellarator_module, only: init_stellarator_module - use stellarator_variables, only: init_stellarator_variables - use tfcoil_variables, only: init_tfcoil_variables - use times_variables, only: init_times_variables - use constants, only: init_constants - use current_drive_variables, only: init_current_drive_variables - use primary_pumping_variables, only: init_primary_pumping_variables - use pfcoil_variables, only: init_pfcoil_variables - use structure_variables, only: init_structure_variables - use vacuum_variables, only: init_vacuum_variables - use pf_power_variables, only: init_pf_power_variables - use build_variables, only: init_build_variables - use constraint_variables, only: init_constraint_variables - use pulse_variables, only: init_pulse_variables - use rebco_variables, only: init_rebco_variables - use reinke_variables, only: init_reinke_variables - use define_iteration_variables, only: init_define_iteration_variables - use reinke_module, only: init_reinke_module - use water_usage_variables, only: init_watuse_variables - use CS_fatigue_variables, only: init_CS_fatigue_variables - use blanket_library, only: init_blanket_library - use dcll_module, only: init_dcll_module - - call init_numerics - call init_input - call init_buildings_variables - call init_cost_variables - call init_divertor_variables - call init_error_handling - call init_fson_library - call init_fwbs_variables - call init_global_variables - call init_ccfe_hcpb_module - call init_heat_transport_variables - call init_ife_variables - call init_impurity_radiation_module - call init_pfcoil_module - call init_physics_module - call init_physics_variables - call init_scan_module - call init_sctfcoil_module - call init_stellarator_module - call init_stellarator_variables - call init_tfcoil_variables - call init_times_variables - call init_constants - call init_current_drive_variables - call init_primary_pumping_variables - call init_pfcoil_variables - call init_structure_variables - call init_vacuum_variables - call init_pf_power_variables - call init_build_variables - call init_constraint_variables - call init_pulse_variables - call init_rebco_variables - call init_reinke_variables - call init_define_iteration_variables - call init_reinke_module - call init_watuse_variables - call init_CS_fatigue_variables - call init_blanket_library - call init_dcll_module - end subroutine init_all_module_vars - - subroutine init - - !! Routine that calls the initialisation routines - !! author: P J Knight, CCFE, Culham Science Centre - !! None - !! This routine calls the main initialisation routines that set - !! the default values for the global variables, reads in data from - !! the input file, and checks the run parameters for consistency. - use global_variables, only: verbose, fileprefix, output_prefix - use main_module, only: run_summary - use constants, only: opt_file, vfile, nout, nplot, mfile, sig_file - use error_handling, only: initialise_error_list - use numerics, only: ixc , lablxc, nvar - use process_input, only: nin, input - use stellarator_module, only: stinit implicit none - ! Arguments - - ! Local variables - integer :: i - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + call init_fson_library - ! Initialise error handling + end subroutine init_fortran_modules - call initialise_error_list + subroutine open_files + use global_variables, only: verbose, fileprefix, output_prefix + use constants, only: nout, mfile + use process_input, only: nin - ! Initialise the program variables - call initial + implicit none ! Open the input/output external files if (trim(fileprefix) == "") then @@ -135,36 +33,11 @@ subroutine init else open(unit=nin,file=trim(fileprefix),status='old') end if - ! open(unit=nin,file=trim(fileprefix)//'IN.DAT',status='old') open(unit=nout ,file=trim(output_prefix)//'OUT.DAT' ,status='unknown') open(unit=mfile ,file=trim(output_prefix)//'MFILE.DAT' ,status='unknown') - ! Input any desired new initial values - call input - - ! Initialise stellarator parameters if necessary - ! This overrides some of the bounds of the tokamak parameters - call stinit - - ! Check input data for errors/ambiguities - call check - - ! Write to the output file certain relevant details about this run - call run_summary - - ! Open verbose diagnostics file - if (verbose == 1) then - open(unit=vfile,file=trim(output_prefix)//'VFILE.DAT',status='unknown') - write(vfile,'(a80)') 'nviter = number of VMCON iterations.' - write(vfile,'(a80)') '(1-mod(ifail,7))=1 indicates that there has '// & - 'been an escape from a failed line search.' - write(vfile,'(a80)') 'odd/even is a convenient plotting bit.' - write(vfile,'(100a13)') 'nviter','escape', 'odd/even', 'te','coe','rmajor', & - 'fusion_power','bt','t_burn','sqsumsq', (lablxc(ixc(i)),i=1,nvar) - end if - - end subroutine init + end subroutine open_files subroutine open_idempotence_files ! Open new output file and mfile to write output to @@ -224,4 +97,755 @@ subroutine finish if (verbose == 1) close(unit = vfile) end subroutine finish + subroutine check + + !! Routine to reset specific variables if certain options are + !! being used + !! author: P J Knight, CCFE, Culham Science Centre + !! None + !! This routine performs a sanity check of the input variables + !! and ensures other dependent variables are given suitable values. + + !! ! + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + use build_variables, only: blnkith, bore, gapoh, ohcth, precomp, iprecomp, & + i_r_cp_top, r_cp_top, vgaptop, vgap_xpoint_divertor, shldtth, shldlth, d_vv_top, d_vv_bot, tf_in_cs + use buildings_variables, only: esbldgm3, triv + use current_drive_variables, only: gamcd, iefrf, irfcd + use error_handling, only: errors_on, idiags, fdiags, report_error + use fwbs_variables, only: breeder_multiplier, iblanket, vfcblkt, vfpblkt, & + iblnkith + use global_variables, only: icase + use heat_transport_variables, only: trithtmw + use ife_variables, only: ife + use impurity_radiation_module, only: nimp, impurity_arr_frac, fimp + use numerics, only: ixc, icc, ioptimz, neqns, nineqns, nvar, boundl, & + boundu + use pfcoil_variables, only: ipfres, ngrp, pfclres, ipfloc, ncls, isumatoh + use physics_variables, only: aspect, f_deuterium, fgwped, f_helium3, & + fgwsep, f_tritium, i_bootstrap_current, i_single_null, i_plasma_current, idivrt, ishape, & + iradloss, isc, ipedestal, ilhthresh, itart, nesep, rhopedn, rhopedt, & + rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar, i_diamagnetic_current + use pulse_variables, only: lpulse + use reinke_variables, only: fzactual, impvardiv + use tfcoil_variables, only: casthi, casthi_is_fraction, casths, i_tf_sup, & + tcoolin, tcpav, tfc_sidewall_is_fraction, tmargmin, tmargmin_cs, & + tmargmin_tf, eff_tf_cryo, eyoung_ins, i_tf_bucking, i_tf_shape, & + n_tf_graded_layers, n_tf_stress_layers, tlegav, i_tf_stress_model, & + i_tf_sc_mat, i_tf_wp_geom, i_tf_turns_integer, tinstf, thwcndut, & + tfinsgap, rcool, dhecoil, thicndut, i_cp_joints, t_turn_tf_is_input, & + t_turn_tf, tftmp, t_cable_tf, t_cable_tf_is_input, tftmp, tmpcry, & + i_tf_cond_eyoung_axial, eyoung_cond_axial, eyoung_cond_trans, & + i_tf_cond_eyoung_trans, i_str_wp + use stellarator_variables, only: istell + use sctfcoil_module, only: initialise_cables + use vacuum_variables, only: vacuum_model + use, intrinsic :: iso_fortran_env, only: dp=>real64 + + implicit none + + ! Local variables + + integer :: i,j,k,imp + real(dp) :: fsum + + real(dp) :: dr_tf_wp_min + !! Minimal WP or conductor layer thickness [m] + ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + errors_on = .true. + + ! Check that there are sufficient iteration variables + if (nvar < neqns) then + idiags(1) = nvar ; idiags(2) = neqns + call report_error(137) + end if + + ! Check that sufficient elements of ixc and icc have been specified + if ( any(ixc(1:nvar) == 0) ) then + idiags(1) = nvar + call report_error(139) + end if + + + if ( any(icc(1:neqns+nineqns) == 0) ) then + idiags(1) = neqns ; idiags(2) = nineqns + call report_error(140) + end if + + ! Deprecate constraints 3 and 4 + if ( any(icc(1:neqns+nineqns) == 3) ) then + call report_error(162) + write(*,*) 'PROCESS stopping' + stop 1 + end if + + if ( any(icc(1:neqns+nineqns) == 4) ) then + call report_error(163) + write(*,*) 'PROCESS stopping' + stop 1 + end if + + + ! MDK Report error if constraint 63 is used with old vacuum model + if (any(icc(1:neqns+nineqns) == 63).and.(vacuum_model.ne.'simple') ) then + write(*,*) 'Constraint 63 is requested without the correct vacuum model ("simple").' + write(*,*) 'vacuum_model = ', vacuum_model + write(*,*) 'PROCESS stopping' + stop 1 + end if + + if ( any(icc(1:neqns+nineqns) == 74) ) then + write(*,*)'Constraint 74 (TF coil quench temperature for Croco HTS conductor) is not yet implemented' + write(*,*) 'PROCESS stopping' + stop 1 + end if + + ! Fuel ion fractions must add up to 1.0 + if (abs(1.0D0 - f_deuterium - f_tritium - f_helium3) > 1.0D-6) then + fdiags(1) = f_deuterium; fdiags(2) = f_tritium ; fdiags(3) = f_helium3 + call report_error(36) + end if + + if (f_tritium < 1.0D-3) then ! tritium fraction is negligible + triv = 0.0D0 + trithtmw = 0.0D0 + end if + + if (fimp(2) .ne. 0.1D0) then + write(*,*)'The thermal alpha/electron density ratio should be controlled using ralpne (itv 109) and not fimp(2).' + write(*,*)'fimp(2) should be removed from the input file, or set to the default value 0.1D0.' + stop 1 + end if + + ! Impurity fractions + do imp = 1,nimp + impurity_arr_frac(imp) = fimp(imp) + end do + + ! The 1/R B field dependency constraint variable is being depreciated + ! Stop the run if the constraint 10 is used + if ( any( icc == 10 ) ) then + call report_error(236) + stop 1 + end if + + ! Stop the run if oacdcp is used as an optimisation variable + ! As the current density is now calculated from bt without constraint 10 + if ( any( ixc == 12 ) ) then + call report_error(236) + stop 1 + end if + + ! Warn if ion power balance equation is being used with the new radiation model + if (any(icc == 3)) then + call report_error(138) + end if + + ! Plasma profile consistency checks + if (ife /= 1) then + if (ipedestal == 1) then + + ! Temperature checks + if (teped < tesep) then + fdiags(1) = teped ; fdiags(2) = tesep + call report_error(146) + end if + + if ((abs(rhopedt-1.0D0) <= 1.0D-7).and.((teped-tesep) >= 1.0D-7)) then + fdiags(1) = rhopedt ; fdiags(2) = teped ; fdiags(3) = tesep + call report_error(147) + end if + + ! Core temperature should always be calculated (later) as being + ! higher than the pedestal temperature, if and only if the + ! volume-averaged temperature never drops below the pedestal + ! temperature. Prevent this by adjusting te, and its lower bound + ! (which will only have an effect if this is an optimisation run) + if (te <= teped) then + fdiags(1) = te ; fdiags(2) = teped + te = teped*1.001D0 + call report_error(149) + end if + + if ((ioptimz >= 0).and.(any(ixc == 4)).and.(boundl(4) < teped*1.001D0)) then + call report_error(150) + boundl(4) = teped*1.001D0 + boundu(4) = max(boundu(4), boundl(4)) + end if + + ! Density checks + ! Case where pedestal density is set manually + ! --------------- + if ( (fgwped < 0) .or. (.not.any(ixc==145)) ) then + + ! Issue #589 Pedestal density is set manually using neped but it is less than nesep. + if ( neped < nesep ) then + fdiags(1) = neped ; fdiags(2) = nesep + call report_error(151) + end if + + ! Issue #589 Pedestal density is set manually using neped, + ! but pedestal width = 0. + if ( (abs(rhopedn-1.0D0) <= 1.0D-7).and.((neped-nesep) >= 1.0D-7) ) then + fdiags(1) = rhopedn ; fdiags(2) = neped ; fdiags(3) = nesep + call report_error(152) + end if + end if + + ! Issue #862 : Variable ne0/neped ratio without constraint eq 81 (ne0>neped) + ! -> Potential hollowed density profile + if ( (ioptimz >= 0) .and. (.not.any(icc==81)) ) then + if ( any(ixc == 145 )) call report_error(154) + if ( any(ixc == 6 )) call report_error(155) + end if + end if + end if + ! --------------- + + + ! Cannot use Psep/R and PsepB/qAR limits at the same time + if(any(icc == 68) .and. any(icc == 56)) then + call report_error(178) + endif + + if ((any(ixc==145)) .and. (boundl(145) < fgwsep)) then !if lower bound of fgwped < fgwsep + fdiags(1) = boundl(145); fdiags(2) = fgwsep + call report_error(186) + end if + + if (any(icc == 78)) then + + !If Reinke criterion is used tesep is calculated and cannot be an + !iteration variable + if (any(ixc == 119)) then + call report_error(219) + endif + + !If Reinke criterion is used need to enforce LH-threshold + !using Martin scaling for consistency + if (.not. ilhthresh == 6) then + call report_error(218) + endif + if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then + call report_error(218) + endif + + + endif + + if (any(icc == 78)) then + + !If Reinke criterion is used tesep is calculated and cannot be an + !iteration variable + if (any(ixc == 119)) then + call report_error(219) + endif + + !If Reinke criterion is used need to enforce LH-threshold + !using Martin scaling for consistency + if (.not. ilhthresh == 6) then + call report_error(218) + endif + if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then + call report_error(218) + endif + + + endif + + if (i_single_null == 0) then + idivrt = 2 + vgaptop = vgap_xpoint_divertor + shldtth = shldlth + d_vv_top = d_vv_bot + call report_error(272) + else ! i_single_null == 1 + idivrt = 1 + end if + + + ! Tight aspect ratio options (ST) + ! -------------------------------- + if ( itart == 1 ) then + + icase = 'Tight aspect ratio tokamak model' + + ! Disabled Forcing that no inboard breeding blanket is used + ! Disabled iblnkith = 0 + + ! Check if the choice of plasma current is addapted for ST + ! 2 : Peng Ip scaling (See STAR code documentation) + ! 9 : Fiesta Ip scaling + if (i_plasma_current /= 2 .and. i_plasma_current /= 9) then + idiags(1) = i_plasma_current ; call report_error(37) + end if + + !! If using Peng and Strickler (1986) model (itartpf == 0) + ! Overwrite the location of the TF coils + ! 2 : PF coil on top of TF coil + ! 3 : PF coil outside of TF coil + if (itartpf == 0) then + ipfloc(1) = 2 + ipfloc(2) = 3 + ipfloc(3) = 3 + end if + + ! Water cooled copper magnets initalisation / checks + if ( i_tf_sup == 0 ) then + ! Check if the initial centrepost coolant loop adapted to the magnet technology + ! Ice cannot flow so tcoolin > 273.15 K + if ( tcoolin < 273.15D0 ) call report_error(234) + + ! Temperature of the TF legs cannot be cooled down + if ( abs(tlegav+1.0D0) > epsilon(tlegav) .and. tlegav < 273.15D0 ) call report_error(239) + + ! Check if conductor upper limit is properly set to 50 K or below + if ( any(ixc == 20 ) .and. boundu(20) < 273.15D0 ) call report_error(241) + + ! Call a lvl 3 error if superconductor magnets are used + else if ( i_tf_sup == 1 ) then + call report_error(233) + + ! Aluminium magnets initalisation / checks + ! Initialize the CP conductor temperature to cryogenic temperature for cryo-al magnets (20 K) + else if ( i_tf_sup == 2 ) then + + ! Call a lvl 3 error if the inlet coolant temperature is too large + ! Motivation : ill-defined aluminium resistivity fit for T > 40-50 K + if ( tcoolin > 40.0D0 ) call report_error(235) + + ! Check if the leg average temperature is low enough for the resisitivity fit + if ( tlegav > 50.0D0 ) call report_error(238) + + ! Check if conductor upper limit is properly set to 50 K or below + if ( any(ixc == 20 ) .and. boundu(20) > 50.0D0 ) call report_error(240) + + ! Otherwise intitialise the average conductor temperature at + tcpav = tcoolin + + end if + + ! Check if the boostrap current selection is addapted to ST + if (i_bootstrap_current == 1) call report_error(38) + + ! Check if a single null divertor is used in double null machine + if (i_single_null == 0 .and. (ftar == 1.0 .or. ftar == 0.0)) then + call report_error(39) + end if + + ! Set the TF coil shape to picture frame (if default value) + if ( i_tf_shape == 0 ) i_tf_shape = 2 + + ! Warning stating that the CP fast neutron fluence calculation + ! is not addapted for cryoaluminium calculations yet + if ( i_tf_sup == 2 .and. any( icc == 85 ) .and. itart == 1 ) then + call report_error(260) + end if + + ! Setting the CP joints default options : + ! 0 : No joints for superconducting magents (i_tf_sup = 1) + ! 1 : Sliding joints for resistive magnets (i_tf_sup = 0, 2) + if ( i_cp_joints == -1 ) then + if ( i_tf_sup == 1 ) then + i_cp_joints = 0 + else + i_cp_joints = 1 + end if + end if + + ! Checking the CP TF top radius + if ( ( abs(r_cp_top) > epsilon(r_cp_top) .or. any(ixc(1:nvar) == 174) ) & + .and. i_r_cp_top /= 1 ) then + call report_error(267) + end if + ! -------------------------------- + + + ! Conventionnal aspect ratios specific + ! ------------------------------------ + else + + if (i_plasma_current == 2 .or. i_plasma_current == 9) call report_error(40) + + ! Set the TF coil shape to PROCESS D-shape (if default value) + if ( i_tf_shape == 0 ) i_tf_shape = 1 + + ! Check PF coil configurations + j = 0 ; k = 0 + do i = 1, ngrp + if ((ipfloc(i) /= 2).and.(ncls(i) /= 2)) then + idiags(1) = i ; idiags(2) = ncls(i) + call report_error(41) + end if + + if (ipfloc(i) == 2) then + j = j + 1 + k = k + ncls(i) + end if + end do + + if (k == 1) call report_error(42) + if (k > 2) call report_error(43) + if ((i_single_null == 1).and.(j < 2)) call report_error(44) + + ! Constraint 10 is dedicated to ST designs with demountable joints + if ( any(icc(1:neqns+nineqns) == 10 ) ) call report_error(259) + + end if + ! ------------------------------------ + + ! Pulsed power plant model + if (lpulse == 1) then + icase = 'Pulsed tokamak model' + else + esbldgm3 = 0.0D0 + end if + + ! Ensure minimum cycle time constraint is turned off + ! (not currently available, as routine thrmal has been commented out) + if ( any(icc == 42) ) then + call report_error(164) + end if + + + + ! TF coil + ! ------- + ! TF stress model not defined of r_tf_inboard = 0 + ! Unless i_tf_stress_model == 2 + ! -> If bore + gapoh + ohcth = 0 and fixed and stress constraint is used + ! Generate a lvl 3 error proposing not to use any stress constraints + if ( ( .not. ( any(ixc == 16 ) .or. any(ixc == 29 ) .or. any(ixc == 42 ) ) ) & ! No bore,gapoh, ohcth iteration + .and. ( abs(bore + gapoh + ohcth + precomp) < epsilon(bore) ) & ! bore + gapoh + ohcth = 0 + .and. ( any(icc == 31) .or. any(icc == 32) ) & ! Stress constraints (31 or 32) is used + .and. ( i_tf_stress_model /= 2 ) ) then ! TF stress model can't handle no bore + + call report_error(246) + stop 1 + end if + + ! Make sure that plane stress model is not used for resistive magnets + if ( i_tf_stress_model == 1 .and. i_tf_sup /= 1 ) call report_error(253) + + ! bucking cylinder default option setting + ! - bucking (casing) for SC i_tf_bucking ( i_tf_bucking = 1 ) + ! - No bucking for copper magnets ( i_tf_bucking = 0 ) + ! - Bucking for aluminium magnets ( i_tf_bucking = 1 ) + if ( i_tf_bucking == -1 ) then + if ( i_tf_sup == 0 ) then + i_tf_bucking = 0 + else + i_tf_bucking = 1 + end if + end if + + ! Ensure that the TF isnt placed against the + ! CS which is now outside it + if ( i_tf_bucking >= 2 .and. tf_in_cs == 1 ) then + call report_error(281) + end if + ! Ensure that no pre-compression structure + ! is used for bucked and wedged design + if ( i_tf_bucking >= 2 .and. iprecomp == 1 ) then + call report_error(252) + end if + + ! Number of stress calculation layers + ! +1 to add in the inboard TF coil case on the plasma side, per Issue #1509 + n_tf_stress_layers = i_tf_bucking + n_tf_graded_layers + 1 + + ! If TFC sidewall has not been set by user + if ( casths < 0.1d-10 ) tfc_sidewall_is_fraction = .true. + + ! If inboard TF coil case plasma side thickness has not been set by user + if( casthi < 0.1d-10 ) casthi_is_fraction = .true. + + ! Setting the default cryo-plants efficiencies + !-! + if ( abs(eff_tf_cryo + 1.0D0) < epsilon(eff_tf_cryo) ) then + + ! The ITER cyoplant efficiency is used for SC + if ( i_tf_sup == 1 ) then + eff_tf_cryo = 0.13D0 + + ! Strawbrige plot extrapolation is used for Cryo-Al + else if ( i_tf_sup == 2 ) then + eff_tf_cryo = 0.40D0 + end if + + ! Cryo-plane efficiency must be in [0-1.0] + else if ( eff_tf_cryo > 1.0D0 .or. eff_tf_cryo < 0.0D0 ) then + call report_error(248) + stop 1 + end if + !-! + + ! Integer turns option not yet available for REBCO taped turns + !-! + if ( i_tf_sc_mat == 6 .and. i_tf_turns_integer == 1 ) then + call report_error(254) + stop 1 + end if + !-! + + + ! Setting up insulation layer young modulae default values [Pa] + !-! + if ( abs(eyoung_ins - 1.0D8 ) < epsilon(eyoung_ins) ) then + + ! Copper magnets, no insulation material defined + ! But use the ITER design by default + if ( i_tf_sup == 0 ) then + eyoung_ins = 20.0D9 + + ! SC magnets + ! Value from DDD11-2 v2 2 (2009) + else if ( i_tf_sup == 1 ) then + eyoung_ins = 20.0D9 + + ! Cryo-aluminum magnets (Kapton polymer) + else if ( i_tf_sup == 2 ) then + eyoung_ins = 2.5D9 + end if + end if + !-! + + !-! Setting the default WP geometry + !-! + if ( i_tf_wp_geom == -1 ) then + if ( i_tf_turns_integer == 0 ) i_tf_wp_geom = 1 + if ( i_tf_turns_integer == 1 ) i_tf_wp_geom = 0 + end if + !-! + + !-! Setting the TF coil conductor elastic properties + !-! + if ( i_tf_cond_eyoung_axial == 0 ) then + ! Conductor stiffness is not considered + eyoung_cond_axial = 0 + eyoung_cond_trans = 0 + else if ( i_tf_cond_eyoung_axial == 2 ) then + ! Select sensible defaults from the literature + select case (i_tf_sc_mat) + case (1,4,5) + ! Nb3Sn: Nyilas, A et. al, Superconductor Science and Technology 16, no. 9 (2003): 1036–42. https://doi.org/10.1088/0953-2048/16/9/313. + eyoung_cond_axial = 32D9 + case (2) + ! Bi-2212: Brown, M. et al, IOP Conference Series: Materials Science and Engineering 279 (2017): 012022. https://doi.org/10.1088/1757-899X/279/1/012022. + eyoung_cond_axial = 80D9 + case (3,7) + ! NbTi: Vedrine, P. et. al, IEEE Transactions on Applied Superconductivity 9, no. 2 (1999): 236–39. https://doi.org/10.1109/77.783280. + eyoung_cond_axial = 6.8D9 + case (6,8,9) + ! REBCO: Fujishiro, H. et. al, Physica C: Superconductivity, 426–431 (2005): 699–704. https://doi.org/10.1016/j.physc.2005.01.045. + eyoung_cond_axial = 145D9 + end select + + if ( i_tf_cond_eyoung_trans == 0) then + ! Transverse stiffness is not considered + eyoung_cond_trans = 0 + else + ! Transverse stiffness is significant + eyoung_cond_trans = eyoung_cond_axial + end if + end if + !-! + + ! Check if the WP/conductor radial thickness (dr_tf_wp) is large enough + ! To contains the insulation, cooling and the support structure + ! Rem : Only verified if the WP thickness is used + if ( any(ixc(1:nvar) == 140) ) then + + ! Minimal WP thickness + if ( i_tf_sup == 1 ) then + dr_tf_wp_min = 2.0D0 * ( tinstf + tfinsgap + thicndut + dhecoil ) + + ! Steel conduit thickness (can be an iteration variable) + if ( any(ixc(1:nvar) == 58 ) ) then + dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * boundl(58) + else + dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * thwcndut + end if + + ! Minimal conductor layer thickness + else if ( i_tf_sup == 0 .or. i_tf_sup == 2 ) then + dr_tf_wp_min = 2.0D0 * ( thicndut + tinstf ) + 4.0D0 * rcool + end if + + if ( boundl(140) < dr_tf_wp_min ) then + fdiags(1) = dr_tf_wp_min + call report_error(255) + end if + end if + + ! Setting t_turn_tf_is_input to true if t_turn_tf is an input + if ( abs(t_turn_tf) < epsilon(t_turn_tf) ) then + t_turn_tf_is_input = .false. + else + t_turn_tf_is_input = .true. + end if + + ! Impossible to set the turn size of integer turn option + if ( t_turn_tf_is_input .and. i_tf_turns_integer == 1 ) then + call report_error(269) + end if + + if ( i_tf_wp_geom /= 0 .and. i_tf_turns_integer == 1 ) then + call report_error(283) + end if + + if ( i_bootstrap_current == 5 .and. i_diamagnetic_current /= 0 ) then + call report_error(284) + end if + + ! Setting t_cable_tf_is_input to true if t_cable_tf is an input + if ( abs(t_cable_tf) < epsilon(t_cable_tf) ) then + t_cable_tf_is_input = .false. + else + t_cable_tf_is_input = .true. + end if + + ! Impossible to set the cable size of integer turn option + if ( t_cable_tf_is_input .and. i_tf_turns_integer == 1 ) then + call report_error(269) + end if + + ! Impossible to set both the TF coil turn and the cable dimension + if ( t_turn_tf_is_input .and. t_cable_tf_is_input ) then + call report_error(271) + end if + + ! Checking the SC temperature for LTS + if ( ( i_tf_sc_mat == 1 .or. & + i_tf_sc_mat == 3 .or. & + i_tf_sc_mat == 4 .or. & + i_tf_sc_mat == 5 ) .and. tftmp > 10.0D0 ) then + call report_error(270) + end if + ! ------- + + + + ! PF coil resistivity is zero if superconducting + if (ipfres == 0) pfclres = 0.0D0 + + ! If there is no NBI, then hot beam density should be zero + if (irfcd == 1) then + if ((iefrf /= 5).and.(iefrf /= 8)) rnbeam = 0.0D0 + else + rnbeam = 0.0D0 + end if + + ! Set inboard blanket thickness to zero if no inboard blanket switch + ! used (Issue #732) + if (iblnkith == 0) blnkith = 0.0D0 + + ! Solid breeder assumed if ipowerflow=0 + + !if (ipowerflow == 0) blkttype = 3 + + ! Set coolant fluid type + + !if ((blkttype == 1).or.(blkttype == 2)) then + ! coolwh = 2 ! water + !else + ! coolwh = 1 ! helium + !end if + + ! But... set coolant to water if blktmodel > 0 + ! Although the *blanket* is by definition helium-cooled in this case, + ! the shield etc. are assumed to be water-cooled, and since water is + ! heavier (and the unit cost of pumping it is higher), the calculation + ! for coolmass is better done with coolwh=2 if blktmodel > 0 to give + ! slightly pessimistic results. + + !if (blktmodel > 0) then + ! secondary_cycle = 0 + ! blkttype = 3 ! HCPB + ! coolwh = 2 + !end if + + ! Ensure that blanket material fractions allow non-zero space for steel + ! CCFE HCPB Model + + if (istell == 0) then + if ((iblanket == 1).or.(iblanket == 3)) then + fsum = breeder_multiplier + vfcblkt + vfpblkt + if (fsum >= 1.0D0) then + idiags(1) = iblanket + fdiags(2) = breeder_multiplier + fdiags(3) = vfcblkt + fdiags(4) = vfpblkt + fdiags(5) = fsum + call report_error(165) + end if + end if + end if + + ! Initialise superconductor cable parameters + if(i_tf_sup==1)then + call initialise_cables() + end if + + ! Check that the temperature margins are not overdetermined + if(tmargmin>0.0001d0)then + ! This limit has been input and will be applied to both TFC and CS + if(tmargmin_tf>0.0001d0)then + write(*,*)'tmargmin_tf and tmargmin should not both be specified in IN.DAT.' + write(*,*)'tmargmin_tf has been ignored.' + end if + if(tmargmin_cs>0.0001d0)then + write(*,*)'tmargmin_cs and tmargmin should not both be specified in IN.DAT.' + write(*,*)'tmargmin_cs has been ignored.' + end if + tmargmin_tf = tmargmin + tmargmin_cs = tmargmin + end if + + if (tauee_in.ge.1.0D-10.and.isc.ne.48) then + ! Report error if confinement time is in the input + ! but the scaling to use it is not selected. + call report_error(220) + end if + + if (aspect.gt.1.7D0.and.isc.eq.46) then + ! NSTX scaling is for A<1.7 + call report_error(221) + end if + + if (i_plasma_current.eq.2.and.isc.eq.42) then + call report_error(222) + end if + + ! Cannot use temperature margin constraint with REBCO TF coils + if(any(icc == 36) .and. ((i_tf_sc_mat == 8).or.(i_tf_sc_mat == 9))) then + call report_error(265) + endif + + ! Cannot use temperature margin constraint with REBCO CS coils + if(any(icc == 60) .and. (isumatoh == 8)) then + call report_error(264) + endif + + ! Cold end of the cryocooler should be colder than the TF + if(tmpcry > tftmp) then + call report_error(273) + endif + + ! Cannot use TF coil strain limit if i_str_wp is off: + if(any(icc == 88) .and. (i_str_wp == 0)) then + call report_error(275) + endif + + errors_on = .false. + + ! Disable error logging only after all checks have been performed. + ! (CPSS #1582: Why is error logging disabled at all?) + errors_on = .false. + + + end subroutine check + end module init_module diff --git a/source/fortran/initial.f90 b/source/fortran/initial.f90 deleted file mode 100755 index 8b75f5f7eb..0000000000 --- a/source/fortran/initial.f90 +++ /dev/null @@ -1,977 +0,0 @@ -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -subroutine initial - - !! Routine to initialise - !! author: P J Knight, CCFE, Culham Science Centre - !! None - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use define_iteration_variables, only: init_itv_1, init_itv_2, init_itv_3, & - init_itv_4, init_itv_5, init_itv_6, init_itv_7, init_itv_8, init_itv_9, & - init_itv_10, init_itv_11, init_itv_12, init_itv_13, init_itv_14, init_itv_15, & - init_itv_16, init_itv_17, init_itv_18, init_itv_19, init_itv_20, init_itv_21, & - init_itv_23, init_itv_25, init_itv_26, init_itv_27, init_itv_28, init_itv_29, & - init_itv_30, init_itv_31, init_itv_32, init_itv_33, init_itv_34, init_itv_35, & - init_itv_36, init_itv_37, init_itv_38, init_itv_39, init_itv_40, init_itv_41, & - init_itv_42, init_itv_44, init_itv_45, init_itv_46, init_itv_47, init_itv_48, & - init_itv_49, init_itv_50, init_itv_51, init_itv_53, init_itv_54, & - init_itv_56, init_itv_57, init_itv_58, init_itv_59, init_itv_60, init_itv_61, & - init_itv_62, init_itv_63, init_itv_64, init_itv_65, init_itv_66, init_itv_67, & - init_itv_68, init_itv_69, init_itv_70, init_itv_71, init_itv_72, init_itv_73, & - init_itv_74, init_itv_75, init_itv_79, init_itv_81, init_itv_82, init_itv_83, & - init_itv_84, init_itv_85, init_itv_86, init_itv_89, init_itv_90, init_itv_91, & - init_itv_92, init_itv_93, init_itv_94, init_itv_95, init_itv_96, init_itv_97, & - init_itv_98, init_itv_103, init_itv_104, init_itv_105, & - init_itv_106, init_itv_107, init_itv_108, init_itv_109, init_itv_110, & - init_itv_111, init_itv_112, init_itv_113, init_itv_114, init_itv_115, & - init_itv_116, init_itv_117, init_itv_118, init_itv_119, init_itv_120, & - init_itv_121, init_itv_122, init_itv_123, init_itv_124, init_itv_125, & - init_itv_126, init_itv_127, init_itv_128, init_itv_129, init_itv_130, & - init_itv_131, init_itv_132, init_itv_133, init_itv_134, init_itv_135, & - init_itv_136, init_itv_137, init_itv_138, init_itv_139, init_itv_140, & - init_itv_141, init_itv_142, init_itv_143, init_itv_144, init_itv_145, & - init_itv_146, init_itv_147, init_itv_148, init_itv_149, & - init_itv_152, init_itv_153, init_itv_154, init_itv_155, & - init_itv_156, init_itv_157, init_itv_158, init_itv_159, init_itv_160, & - init_itv_161, init_itv_162, init_itv_163, init_itv_164, init_itv_165, & - init_itv_166, init_itv_167, init_itv_168, init_itv_169, init_itv_170, & - init_itv_171, init_itv_172, init_itv_173, init_itv_174, init_itv_175 -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - - implicit none - - ! Arguments - - ! Local variables - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - !! boundl(ipnvars) /../ : lower bounds on iteration variables - !! boundu(ipnvars) /../ : upper bounds on iteration variables - - ! Issue #287 The initialization subroutines for the iteration variables are called - call init_itv_1 - call init_itv_2 - call init_itv_3 - call init_itv_4 - call init_itv_5 - call init_itv_6 - call init_itv_7 - call init_itv_8 - call init_itv_9 - call init_itv_10 - call init_itv_11 - call init_itv_12 - call init_itv_13 - call init_itv_14 - call init_itv_15 - call init_itv_16 - call init_itv_17 - call init_itv_18 - call init_itv_19 - call init_itv_20 - call init_itv_21 - - call init_itv_23 - - call init_itv_25 - call init_itv_26 - call init_itv_27 - call init_itv_28 - call init_itv_29 - call init_itv_30 - call init_itv_31 - call init_itv_32 - call init_itv_33 - call init_itv_34 - call init_itv_35 - call init_itv_36 - call init_itv_37 - call init_itv_38 - call init_itv_39 - call init_itv_40 - call init_itv_41 - call init_itv_42 - - call init_itv_44 - call init_itv_45 - call init_itv_46 - call init_itv_47 - call init_itv_48 - call init_itv_49 - call init_itv_50 - call init_itv_51 - - call init_itv_53 - call init_itv_54 - - call init_itv_56 - call init_itv_57 - call init_itv_58 - call init_itv_59 - call init_itv_60 - call init_itv_61 - call init_itv_62 - call init_itv_63 - call init_itv_64 - call init_itv_65 - call init_itv_66 - call init_itv_67 - call init_itv_68 - call init_itv_69 - call init_itv_70 - call init_itv_71 - call init_itv_72 - call init_itv_73 - call init_itv_74 - call init_itv_75 - - call init_itv_79 - - call init_itv_81 - call init_itv_82 - call init_itv_83 - - call init_itv_85 - call init_itv_86 - - call init_itv_89 - call init_itv_90 - call init_itv_91 - call init_itv_92 - call init_itv_93 - call init_itv_94 - call init_itv_95 - call init_itv_96 - call init_itv_97 - call init_itv_98 - !Not used - call init_itv_103 - call init_itv_104 - call init_itv_105 - call init_itv_106 - call init_itv_107 - call init_itv_108 - call init_itv_109 - call init_itv_110 - call init_itv_111 - call init_itv_112 - call init_itv_113 - call init_itv_114 - call init_itv_115 - call init_itv_116 - call init_itv_117 - call init_itv_118 - call init_itv_119 - call init_itv_120 - call init_itv_121 - call init_itv_122 - call init_itv_123 - call init_itv_124 - call init_itv_125 - call init_itv_126 - call init_itv_127 - call init_itv_128 - call init_itv_129 - call init_itv_130 - call init_itv_131 - call init_itv_132 - call init_itv_133 - call init_itv_134 - call init_itv_135 - call init_itv_136 - call init_itv_137 - call init_itv_138 - call init_itv_139 - call init_itv_140 - call init_itv_141 - call init_itv_142 - call init_itv_143 - call init_itv_144 - call init_itv_145 - call init_itv_146 - call init_itv_147 - call init_itv_148 - call init_itv_149 - call init_itv_152 - call init_itv_153 - call init_itv_154 - call init_itv_155 - call init_itv_156 - call init_itv_157 - call init_itv_158 - call init_itv_159 - call init_itv_160 - call init_itv_161 - call init_itv_162 - call init_itv_163 - call init_itv_164 - call init_itv_165 - call init_itv_166 - call init_itv_167 - call init_itv_168 - call init_itv_169 - call init_itv_170 - call init_itv_171 - call init_itv_172 - call init_itv_173 - call init_itv_174 - call init_itv_175 - - -end subroutine initial - -subroutine check - - !! Routine to reset specific variables if certain options are - !! being used - !! author: P J Knight, CCFE, Culham Science Centre - !! None - !! This routine performs a sanity check of the input variables - !! and ensures other dependent variables are given suitable values. - - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use build_variables, only: blnkith, bore, gapoh, ohcth, precomp, iprecomp, & - i_r_cp_top, r_cp_top, vgaptop, vgap_xpoint_divertor, shldtth, shldlth, d_vv_top, d_vv_bot, tf_in_cs - use buildings_variables, only: esbldgm3, triv - use current_drive_variables, only: gamcd, iefrf, irfcd - use error_handling, only: errors_on, idiags, fdiags, report_error - use fwbs_variables, only: breeder_multiplier, iblanket, vfcblkt, vfpblkt, & - iblnkith - use global_variables, only: icase - use heat_transport_variables, only: trithtmw - use ife_variables, only: ife - use impurity_radiation_module, only: nimp, impurity_arr_frac, fimp - use numerics, only: ixc, icc, ioptimz, neqns, nineqns, nvar, boundl, & - boundu - use pfcoil_variables, only: ipfres, ngrp, pfclres, ipfloc, ncls, isumatoh - use physics_variables, only: aspect, f_deuterium, fgwped, f_helium3, & - fgwsep, f_tritium, i_bootstrap_current, i_single_null, i_plasma_current, idivrt, ishape, & - iradloss, isc, ipedestal, ilhthresh, itart, nesep, rhopedn, rhopedt, & - rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar, i_diamagnetic_current - use pulse_variables, only: lpulse - use reinke_variables, only: fzactual, impvardiv - use tfcoil_variables, only: casthi, casthi_is_fraction, casths, i_tf_sup, & - tcoolin, tcpav, tfc_sidewall_is_fraction, tmargmin, tmargmin_cs, & - tmargmin_tf, eff_tf_cryo, eyoung_ins, i_tf_bucking, i_tf_shape, & - n_tf_graded_layers, n_tf_stress_layers, tlegav, i_tf_stress_model, & - i_tf_sc_mat, i_tf_wp_geom, i_tf_turns_integer, tinstf, thwcndut, & - tfinsgap, rcool, dhecoil, thicndut, i_cp_joints, t_turn_tf_is_input, & - t_turn_tf, tftmp, t_cable_tf, t_cable_tf_is_input, tftmp, tmpcry, & - i_tf_cond_eyoung_axial, eyoung_cond_axial, eyoung_cond_trans, & - i_tf_cond_eyoung_trans, i_str_wp - use stellarator_variables, only: istell - use sctfcoil_module, only: initialise_cables - use vacuum_variables, only: vacuum_model - use, intrinsic :: iso_fortran_env, only: dp=>real64 - - implicit none - - ! Local variables - - integer :: i,j,k,imp - real(dp) :: fsum - - real(dp) :: dr_tf_wp_min - !! Minimal WP or conductor layer thickness [m] - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - errors_on = .true. - - ! Check that there are sufficient iteration variables - if (nvar < neqns) then - idiags(1) = nvar ; idiags(2) = neqns - call report_error(137) - end if - - ! Check that sufficient elements of ixc and icc have been specified - if ( any(ixc(1:nvar) == 0) ) then - idiags(1) = nvar - call report_error(139) - end if - - - if ( any(icc(1:neqns+nineqns) == 0) ) then - idiags(1) = neqns ; idiags(2) = nineqns - call report_error(140) - end if - - ! Deprecate constraints 3 and 4 - if ( any(icc(1:neqns+nineqns) == 3) ) then - call report_error(162) - write(*,*) 'PROCESS stopping' - stop 1 - end if - - if ( any(icc(1:neqns+nineqns) == 4) ) then - call report_error(163) - write(*,*) 'PROCESS stopping' - stop 1 - end if - - - ! MDK Report error if constraint 63 is used with old vacuum model - if (any(icc(1:neqns+nineqns) == 63).and.(vacuum_model.ne.'simple') ) then - write(*,*) 'Constraint 63 is requested without the correct vacuum model ("simple").' - write(*,*) 'vacuum_model = ', vacuum_model - write(*,*) 'PROCESS stopping' - stop 1 - end if - - if ( any(icc(1:neqns+nineqns) == 74) ) then - write(*,*)'Constraint 74 (TF coil quench temperature for Croco HTS conductor) is not yet implemented' - write(*,*) 'PROCESS stopping' - stop 1 - end if - - ! Fuel ion fractions must add up to 1.0 - if (abs(1.0D0 - f_deuterium - f_tritium - f_helium3) > 1.0D-6) then - fdiags(1) = f_deuterium; fdiags(2) = f_tritium ; fdiags(3) = f_helium3 - call report_error(36) - end if - - if (f_tritium < 1.0D-3) then ! tritium fraction is negligible - triv = 0.0D0 - trithtmw = 0.0D0 - end if - - if (fimp(2) .ne. 0.1D0) then - write(*,*)'The thermal alpha/electron density ratio should be controlled using ralpne (itv 109) and not fimp(2).' - write(*,*)'fimp(2) should be removed from the input file, or set to the default value 0.1D0.' - stop 1 - end if - - ! Impurity fractions - do imp = 1,nimp - impurity_arr_frac(imp) = fimp(imp) - end do - - ! The 1/R B field dependency constraint variable is being depreciated - ! Stop the run if the constraint 10 is used - if ( any( icc == 10 ) ) then - call report_error(236) - stop 1 - end if - - ! Stop the run if oacdcp is used as an optimisation variable - ! As the current density is now calculated from bt without constraint 10 - if ( any( ixc == 12 ) ) then - call report_error(236) - stop 1 - end if - - ! Warn if ion power balance equation is being used with the new radiation model - if (any(icc == 3)) then - call report_error(138) - end if - - ! Plasma profile consistency checks - if (ife /= 1) then - if (ipedestal == 1) then - - ! Temperature checks - if (teped < tesep) then - fdiags(1) = teped ; fdiags(2) = tesep - call report_error(146) - end if - - if ((abs(rhopedt-1.0D0) <= 1.0D-7).and.((teped-tesep) >= 1.0D-7)) then - fdiags(1) = rhopedt ; fdiags(2) = teped ; fdiags(3) = tesep - call report_error(147) - end if - - ! Core temperature should always be calculated (later) as being - ! higher than the pedestal temperature, if and only if the - ! volume-averaged temperature never drops below the pedestal - ! temperature. Prevent this by adjusting te, and its lower bound - ! (which will only have an effect if this is an optimisation run) - if (te <= teped) then - fdiags(1) = te ; fdiags(2) = teped - te = teped*1.001D0 - call report_error(149) - end if - - if ((ioptimz >= 0).and.(any(ixc == 4)).and.(boundl(4) < teped*1.001D0)) then - call report_error(150) - boundl(4) = teped*1.001D0 - boundu(4) = max(boundu(4), boundl(4)) - end if - - ! Density checks - ! Case where pedestal density is set manually - ! --------------- - if ( (fgwped < 0) .or. (.not.any(ixc==145)) ) then - - ! Issue #589 Pedestal density is set manually using neped but it is less than nesep. - if ( neped < nesep ) then - fdiags(1) = neped ; fdiags(2) = nesep - call report_error(151) - end if - - ! Issue #589 Pedestal density is set manually using neped, - ! but pedestal width = 0. - if ( (abs(rhopedn-1.0D0) <= 1.0D-7).and.((neped-nesep) >= 1.0D-7) ) then - fdiags(1) = rhopedn ; fdiags(2) = neped ; fdiags(3) = nesep - call report_error(152) - end if - end if - - ! Issue #862 : Variable ne0/neped ratio without constraint eq 81 (ne0>neped) - ! -> Potential hollowed density profile - if ( (ioptimz >= 0) .and. (.not.any(icc==81)) ) then - if ( any(ixc == 145 )) call report_error(154) - if ( any(ixc == 6 )) call report_error(155) - end if - end if - end if - ! --------------- - - - ! Cannot use Psep/R and PsepB/qAR limits at the same time - if(any(icc == 68) .and. any(icc == 56)) then - call report_error(178) - endif - - if ((any(ixc==145)) .and. (boundl(145) < fgwsep)) then !if lower bound of fgwped < fgwsep - fdiags(1) = boundl(145); fdiags(2) = fgwsep - call report_error(186) - end if - - if (any(icc == 78)) then - - !If Reinke criterion is used tesep is calculated and cannot be an - !iteration variable - if (any(ixc == 119)) then - call report_error(219) - endif - - !If Reinke criterion is used need to enforce LH-threshold - !using Martin scaling for consistency - if (.not. ilhthresh == 6) then - call report_error(218) - endif - if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then - call report_error(218) - endif - - - endif - - if (any(icc == 78)) then - - !If Reinke criterion is used tesep is calculated and cannot be an - !iteration variable - if (any(ixc == 119)) then - call report_error(219) - endif - - !If Reinke criterion is used need to enforce LH-threshold - !using Martin scaling for consistency - if (.not. ilhthresh == 6) then - call report_error(218) - endif - if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then - call report_error(218) - endif - - - endif - - if (i_single_null == 0) then - idivrt = 2 - vgaptop = vgap_xpoint_divertor - shldtth = shldlth - d_vv_top = d_vv_bot - call report_error(272) - else ! i_single_null == 1 - idivrt = 1 - end if - - - ! Tight aspect ratio options (ST) - ! -------------------------------- - if ( itart == 1 ) then - - icase = 'Tight aspect ratio tokamak model' - - ! Disabled Forcing that no inboard breeding blanket is used - ! Disabled iblnkith = 0 - - ! Check if the choice of plasma current is addapted for ST - ! 2 : Peng Ip scaling (See STAR code documentation) - ! 9 : Fiesta Ip scaling - if (i_plasma_current /= 2 .and. i_plasma_current /= 9) then - idiags(1) = i_plasma_current ; call report_error(37) - end if - - !! If using Peng and Strickler (1986) model (itartpf == 0) - ! Overwrite the location of the TF coils - ! 2 : PF coil on top of TF coil - ! 3 : PF coil outside of TF coil - if (itartpf == 0) then - ipfloc(1) = 2 - ipfloc(2) = 3 - ipfloc(3) = 3 - end if - - ! Water cooled copper magnets initalisation / checks - if ( i_tf_sup == 0 ) then - ! Check if the initial centrepost coolant loop adapted to the magnet technology - ! Ice cannot flow so tcoolin > 273.15 K - if ( tcoolin < 273.15D0 ) call report_error(234) - - ! Temperature of the TF legs cannot be cooled down - if ( abs(tlegav+1.0D0) > epsilon(tlegav) .and. tlegav < 273.15D0 ) call report_error(239) - - ! Check if conductor upper limit is properly set to 50 K or below - if ( any(ixc == 20 ) .and. boundu(20) < 273.15D0 ) call report_error(241) - - ! Call a lvl 3 error if superconductor magnets are used - else if ( i_tf_sup == 1 ) then - call report_error(233) - - ! Aluminium magnets initalisation / checks - ! Initialize the CP conductor temperature to cryogenic temperature for cryo-al magnets (20 K) - else if ( i_tf_sup == 2 ) then - - ! Call a lvl 3 error if the inlet coolant temperature is too large - ! Motivation : ill-defined aluminium resistivity fit for T > 40-50 K - if ( tcoolin > 40.0D0 ) call report_error(235) - - ! Check if the leg average temperature is low enough for the resisitivity fit - if ( tlegav > 50.0D0 ) call report_error(238) - - ! Check if conductor upper limit is properly set to 50 K or below - if ( any(ixc == 20 ) .and. boundu(20) > 50.0D0 ) call report_error(240) - - ! Otherwise intitialise the average conductor temperature at - tcpav = tcoolin - - end if - - ! Check if the boostrap current selection is addapted to ST - if (i_bootstrap_current == 1) call report_error(38) - - ! Check if a single null divertor is used in double null machine - if (i_single_null == 0 .and. (ftar == 1.0 .or. ftar == 0.0)) then - call report_error(39) - end if - - ! Set the TF coil shape to picture frame (if default value) - if ( i_tf_shape == 0 ) i_tf_shape = 2 - - ! Warning stating that the CP fast neutron fluence calculation - ! is not addapted for cryoaluminium calculations yet - if ( i_tf_sup == 2 .and. any( icc == 85 ) .and. itart == 1 ) then - call report_error(260) - end if - - ! Setting the CP joints default options : - ! 0 : No joints for superconducting magents (i_tf_sup = 1) - ! 1 : Sliding joints for resistive magnets (i_tf_sup = 0, 2) - if ( i_cp_joints == -1 ) then - if ( i_tf_sup == 1 ) then - i_cp_joints = 0 - else - i_cp_joints = 1 - end if - end if - - ! Checking the CP TF top radius - if ( ( abs(r_cp_top) > epsilon(r_cp_top) .or. any(ixc(1:nvar) == 174) ) & - .and. i_r_cp_top /= 1 ) then - call report_error(267) - end if - ! -------------------------------- - - - ! Conventionnal aspect ratios specific - ! ------------------------------------ - else - - if (i_plasma_current == 2 .or. i_plasma_current == 9) call report_error(40) - - ! Set the TF coil shape to PROCESS D-shape (if default value) - if ( i_tf_shape == 0 ) i_tf_shape = 1 - - ! Check PF coil configurations - j = 0 ; k = 0 - do i = 1, ngrp - if ((ipfloc(i) /= 2).and.(ncls(i) /= 2)) then - idiags(1) = i ; idiags(2) = ncls(i) - call report_error(41) - end if - - if (ipfloc(i) == 2) then - j = j + 1 - k = k + ncls(i) - end if - end do - - if (k == 1) call report_error(42) - if (k > 2) call report_error(43) - if ((i_single_null == 1).and.(j < 2)) call report_error(44) - - ! Constraint 10 is dedicated to ST designs with demountable joints - if ( any(icc(1:neqns+nineqns) == 10 ) ) call report_error(259) - - end if - ! ------------------------------------ - - ! Pulsed power plant model - if (lpulse == 1) then - icase = 'Pulsed tokamak model' - else - esbldgm3 = 0.0D0 - end if - - ! Ensure minimum cycle time constraint is turned off - ! (not currently available, as routine thrmal has been commented out) - if ( any(icc == 42) ) then - call report_error(164) - end if - - - - ! TF coil - ! ------- - ! TF stress model not defined of r_tf_inboard = 0 - ! Unless i_tf_stress_model == 2 - ! -> If bore + gapoh + ohcth = 0 and fixed and stress constraint is used - ! Generate a lvl 3 error proposing not to use any stress constraints - if ( ( .not. ( any(ixc == 16 ) .or. any(ixc == 29 ) .or. any(ixc == 42 ) ) ) & ! No bore,gapoh, ohcth iteration - .and. ( abs(bore + gapoh + ohcth + precomp) < epsilon(bore) ) & ! bore + gapoh + ohcth = 0 - .and. ( any(icc == 31) .or. any(icc == 32) ) & ! Stress constraints (31 or 32) is used - .and. ( i_tf_stress_model /= 2 ) ) then ! TF stress model can't handle no bore - - call report_error(246) - stop 1 - end if - - ! Make sure that plane stress model is not used for resistive magnets - if ( i_tf_stress_model == 1 .and. i_tf_sup /= 1 ) call report_error(253) - - ! bucking cylinder default option setting - ! - bucking (casing) for SC i_tf_bucking ( i_tf_bucking = 1 ) - ! - No bucking for copper magnets ( i_tf_bucking = 0 ) - ! - Bucking for aluminium magnets ( i_tf_bucking = 1 ) - if ( i_tf_bucking == -1 ) then - if ( i_tf_sup == 0 ) then - i_tf_bucking = 0 - else - i_tf_bucking = 1 - end if - end if - - ! Ensure that the TF isnt placed against the - ! CS which is now outside it - if ( i_tf_bucking >= 2 .and. tf_in_cs == 1 ) then - call report_error(281) - end if - ! Ensure that no pre-compression structure - ! is used for bucked and wedged design - if ( i_tf_bucking >= 2 .and. iprecomp == 1 ) then - call report_error(252) - end if - - ! Number of stress calculation layers - ! +1 to add in the inboard TF coil case on the plasma side, per Issue #1509 - n_tf_stress_layers = i_tf_bucking + n_tf_graded_layers + 1 - - ! If TFC sidewall has not been set by user - if ( casths < 0.1d-10 ) tfc_sidewall_is_fraction = .true. - - ! If inboard TF coil case plasma side thickness has not been set by user - if( casthi < 0.1d-10 ) casthi_is_fraction = .true. - - ! Setting the default cryo-plants efficiencies - !-! - if ( abs(eff_tf_cryo + 1.0D0) < epsilon(eff_tf_cryo) ) then - - ! The ITER cyoplant efficiency is used for SC - if ( i_tf_sup == 1 ) then - eff_tf_cryo = 0.13D0 - - ! Strawbrige plot extrapolation is used for Cryo-Al - else if ( i_tf_sup == 2 ) then - eff_tf_cryo = 0.40D0 - end if - - ! Cryo-plane efficiency must be in [0-1.0] - else if ( eff_tf_cryo > 1.0D0 .or. eff_tf_cryo < 0.0D0 ) then - call report_error(248) - stop 1 - end if - !-! - - ! Integer turns option not yet available for REBCO taped turns - !-! - if ( i_tf_sc_mat == 6 .and. i_tf_turns_integer == 1 ) then - call report_error(254) - stop 1 - end if - !-! - - - ! Setting up insulation layer young modulae default values [Pa] - !-! - if ( abs(eyoung_ins - 1.0D8 ) < epsilon(eyoung_ins) ) then - - ! Copper magnets, no insulation material defined - ! But use the ITER design by default - if ( i_tf_sup == 0 ) then - eyoung_ins = 20.0D9 - - ! SC magnets - ! Value from DDD11-2 v2 2 (2009) - else if ( i_tf_sup == 1 ) then - eyoung_ins = 20.0D9 - - ! Cryo-aluminum magnets (Kapton polymer) - else if ( i_tf_sup == 2 ) then - eyoung_ins = 2.5D9 - end if - end if - !-! - - !-! Setting the default WP geometry - !-! - if ( i_tf_wp_geom == -1 ) then - if ( i_tf_turns_integer == 0 ) i_tf_wp_geom = 1 - if ( i_tf_turns_integer == 1 ) i_tf_wp_geom = 0 - end if - !-! - - !-! Setting the TF coil conductor elastic properties - !-! - if ( i_tf_cond_eyoung_axial == 0 ) then - ! Conductor stiffness is not considered - eyoung_cond_axial = 0 - eyoung_cond_trans = 0 - else if ( i_tf_cond_eyoung_axial == 2 ) then - ! Select sensible defaults from the literature - select case (i_tf_sc_mat) - case (1,4,5) - ! Nb3Sn: Nyilas, A et. al, Superconductor Science and Technology 16, no. 9 (2003): 1036–42. https://doi.org/10.1088/0953-2048/16/9/313. - eyoung_cond_axial = 32D9 - case (2) - ! Bi-2212: Brown, M. et al, IOP Conference Series: Materials Science and Engineering 279 (2017): 012022. https://doi.org/10.1088/1757-899X/279/1/012022. - eyoung_cond_axial = 80D9 - case (3,7) - ! NbTi: Vedrine, P. et. al, IEEE Transactions on Applied Superconductivity 9, no. 2 (1999): 236–39. https://doi.org/10.1109/77.783280. - eyoung_cond_axial = 6.8D9 - case (6,8,9) - ! REBCO: Fujishiro, H. et. al, Physica C: Superconductivity, 426–431 (2005): 699–704. https://doi.org/10.1016/j.physc.2005.01.045. - eyoung_cond_axial = 145D9 - end select - - if ( i_tf_cond_eyoung_trans == 0) then - ! Transverse stiffness is not considered - eyoung_cond_trans = 0 - else - ! Transverse stiffness is significant - eyoung_cond_trans = eyoung_cond_axial - end if - end if - !-! - - ! Check if the WP/conductor radial thickness (dr_tf_wp) is large enough - ! To contains the insulation, cooling and the support structure - ! Rem : Only verified if the WP thickness is used - if ( any(ixc(1:nvar) == 140) ) then - - ! Minimal WP thickness - if ( i_tf_sup == 1 ) then - dr_tf_wp_min = 2.0D0 * ( tinstf + tfinsgap + thicndut + dhecoil ) - - ! Steel conduit thickness (can be an iteration variable) - if ( any(ixc(1:nvar) == 58 ) ) then - dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * boundl(58) - else - dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * thwcndut - end if - - ! Minimal conductor layer thickness - else if ( i_tf_sup == 0 .or. i_tf_sup == 2 ) then - dr_tf_wp_min = 2.0D0 * ( thicndut + tinstf ) + 4.0D0 * rcool - end if - - if ( boundl(140) < dr_tf_wp_min ) then - fdiags(1) = dr_tf_wp_min - call report_error(255) - end if - end if - - ! Setting t_turn_tf_is_input to true if t_turn_tf is an input - if ( abs(t_turn_tf) < epsilon(t_turn_tf) ) then - t_turn_tf_is_input = .false. - else - t_turn_tf_is_input = .true. - end if - - ! Impossible to set the turn size of integer turn option - if ( t_turn_tf_is_input .and. i_tf_turns_integer == 1 ) then - call report_error(269) - end if - - if ( i_tf_wp_geom /= 0 .and. i_tf_turns_integer == 1 ) then - call report_error(283) - end if - - if ( i_bootstrap_current == 5 .and. i_diamagnetic_current /= 0 ) then - call report_error(284) - end if - - ! Setting t_cable_tf_is_input to true if t_cable_tf is an input - if ( abs(t_cable_tf) < epsilon(t_cable_tf) ) then - t_cable_tf_is_input = .false. - else - t_cable_tf_is_input = .true. - end if - - ! Impossible to set the cable size of integer turn option - if ( t_cable_tf_is_input .and. i_tf_turns_integer == 1 ) then - call report_error(269) - end if - - ! Impossible to set both the TF coil turn and the cable dimension - if ( t_turn_tf_is_input .and. t_cable_tf_is_input ) then - call report_error(271) - end if - - ! Checking the SC temperature for LTS - if ( ( i_tf_sc_mat == 1 .or. & - i_tf_sc_mat == 3 .or. & - i_tf_sc_mat == 4 .or. & - i_tf_sc_mat == 5 ) .and. tftmp > 10.0D0 ) then - call report_error(270) - end if - ! ------- - - - - ! PF coil resistivity is zero if superconducting - if (ipfres == 0) pfclres = 0.0D0 - - ! If there is no NBI, then hot beam density should be zero - if (irfcd == 1) then - if ((iefrf /= 5).and.(iefrf /= 8)) rnbeam = 0.0D0 - else - rnbeam = 0.0D0 - end if - - ! Set inboard blanket thickness to zero if no inboard blanket switch - ! used (Issue #732) - if (iblnkith == 0) blnkith = 0.0D0 - - ! Solid breeder assumed if ipowerflow=0 - - !if (ipowerflow == 0) blkttype = 3 - - ! Set coolant fluid type - - !if ((blkttype == 1).or.(blkttype == 2)) then - ! coolwh = 2 ! water - !else - ! coolwh = 1 ! helium - !end if - - ! But... set coolant to water if blktmodel > 0 - ! Although the *blanket* is by definition helium-cooled in this case, - ! the shield etc. are assumed to be water-cooled, and since water is - ! heavier (and the unit cost of pumping it is higher), the calculation - ! for coolmass is better done with coolwh=2 if blktmodel > 0 to give - ! slightly pessimistic results. - - !if (blktmodel > 0) then - ! secondary_cycle = 0 - ! blkttype = 3 ! HCPB - ! coolwh = 2 - !end if - - ! Ensure that blanket material fractions allow non-zero space for steel - ! CCFE HCPB Model - - if (istell == 0) then - if ((iblanket == 1).or.(iblanket == 3)) then - fsum = breeder_multiplier + vfcblkt + vfpblkt - if (fsum >= 1.0D0) then - idiags(1) = iblanket - fdiags(2) = breeder_multiplier - fdiags(3) = vfcblkt - fdiags(4) = vfpblkt - fdiags(5) = fsum - call report_error(165) - end if - end if - end if - - ! Initialise superconductor cable parameters - if(i_tf_sup==1)then - call initialise_cables() - end if - - ! Check that the temperature margins are not overdetermined - if(tmargmin>0.0001d0)then - ! This limit has been input and will be applied to both TFC and CS - if(tmargmin_tf>0.0001d0)then - write(*,*)'tmargmin_tf and tmargmin should not both be specified in IN.DAT.' - write(*,*)'tmargmin_tf has been ignored.' - end if - if(tmargmin_cs>0.0001d0)then - write(*,*)'tmargmin_cs and tmargmin should not both be specified in IN.DAT.' - write(*,*)'tmargmin_cs has been ignored.' - end if - tmargmin_tf = tmargmin - tmargmin_cs = tmargmin - end if - - if (tauee_in.ge.1.0D-10.and.isc.ne.48) then - ! Report error if confinement time is in the input - ! but the scaling to use it is not selected. - call report_error(220) - end if - - if (aspect.gt.1.7D0.and.isc.eq.46) then - ! NSTX scaling is for A<1.7 - call report_error(221) - end if - - if (i_plasma_current.eq.2.and.isc.eq.42) then - call report_error(222) - end if - - ! Cannot use temperature margin constraint with REBCO TF coils - if(any(icc == 36) .and. ((i_tf_sc_mat == 8).or.(i_tf_sc_mat == 9))) then - call report_error(265) - endif - - ! Cannot use temperature margin constraint with REBCO CS coils - if(any(icc == 60) .and. (isumatoh == 8)) then - call report_error(264) - endif - - ! Cold end of the cryocooler should be colder than the TF - if(tmpcry > tftmp) then - call report_error(273) - endif - - ! Cannot use TF coil strain limit if i_str_wp is off: - if(any(icc == 88) .and. (i_str_wp == 0)) then - call report_error(275) - endif - - errors_on = .false. - - ! Disable error logging only after all checks have been performed. - ! (CPSS #1582: Why is error logging disabled at all?) - errors_on = .false. - - -end subroutine check diff --git a/tests/integration/test_vmcon.py b/tests/integration/test_vmcon.py index 00022265a6..9147f40e2f 100644 --- a/tests/integration/test_vmcon.py +++ b/tests/integration/test_vmcon.py @@ -5,14 +5,15 @@ Expected answers for tests 1 to 3 are given in VMCON documentation ANL-80-64 """ -from process.evaluators import Evaluators -from process.fortran import init_module -from process.fortran import error_handling + import pytest import numpy as np import logging from abc import ABC, abstractmethod from process.solver import get_solver +from process.init import init_all_module_vars +from process.evaluators import Evaluators +from process.fortran import error_handling # Debug-level terminal output logging logger = logging.getLogger(__name__) @@ -24,7 +25,7 @@ @pytest.fixture(autouse=True) def reinit(): """Re-initialise Fortran module variables before each test is run.""" - init_module.init_all_module_vars() + init_all_module_vars() class Case: diff --git a/tests/unit/conftest.py b/tests/unit/conftest.py index 3fa07185b4..bd3ff4d842 100644 --- a/tests/unit/conftest.py +++ b/tests/unit/conftest.py @@ -2,9 +2,10 @@ Define fixtures that will be shared across unit test modules. """ + import pytest -from process import fortran from pathlib import Path +from process.init import init_all_module_vars @pytest.fixture(scope="module", autouse=True) @@ -17,7 +18,7 @@ def reinit_fix(): the module variable values. autouse ensures that this fixture is used automatically by any test function in the unit directory. """ - fortran.init_module.init_all_module_vars() + init_all_module_vars() @pytest.fixture() diff --git a/tests/unit/test_availability.py b/tests/unit/test_availability.py index 6c08ae3649..6c7df164a8 100644 --- a/tests/unit/test_availability.py +++ b/tests/unit/test_availability.py @@ -10,6 +10,7 @@ from process.fortran import times_variables as tv from process.fortran import ife_variables as ifev from process.fortran import divertor_variables as dv +from process.init import init_all_module_vars import pytest from pytest import approx @@ -81,7 +82,7 @@ def test_avail_1(monkeypatch, availability): :type availability: tests.unit.test_availability.availability (functional fixture) """ # Initialise fortran variables to keep test isolated from others - fortran.init_module.init_all_module_vars() + init_all_module_vars() # Mock module vars monkeypatch.setattr(cv, "iavail", 1) @@ -104,7 +105,7 @@ def test_avail_1(monkeypatch, availability): assert pytest.approx(cfactr_exp) == cfactr_obs # Initialise fortran variables again to reset for other tests - fortran.init_module.init_all_module_vars() + init_all_module_vars() def test_calc_u_unplanned_hcd(availability): @@ -486,6 +487,7 @@ def test_avail_2(monkeypatch, availability): :param availability: fixture containing an initialised `Availability` object :type availability: tests.unit.test_availability.availability (functional fixture) """ + # Mock return values for for functions called in avail_2 def mock_calc_u_planned(*args, **kwargs): return 0.01 @@ -567,8 +569,7 @@ def test_avail_st(monkeypatch, availability): :type availability: tests.unit.test_availability.availability (functional fixture) """ # Initialise fortran variables to keep test isolated from others - fortran.init_module.init_all_module_vars() - + init_all_module_vars() monkeypatch.setattr(cv, "tmain", 1.0) monkeypatch.setattr(cv, "tlife", 30.0) monkeypatch.setattr(cv, "u_unplanned_cp", 0.05) @@ -588,9 +589,6 @@ def test_avail_st(monkeypatch, availability): assert pytest.approx(cv.cfactr) == 0.27008858 assert pytest.approx(cv.cpfact, abs=1.0e-8) == 0.00015005 - # Initialise fortran variables again to reset for other tests - fortran.init_module.init_all_module_vars() - @pytest.mark.parametrize("i_tf_sup, exp", ((1, 6.337618), (0, 4))) def test_cp_lifetime(monkeypatch, availability, i_tf_sup, exp): diff --git a/tests/unit/test_input.py b/tests/unit/test_input.py index 96ddbd8656..8756a02cca 100644 --- a/tests/unit/test_input.py +++ b/tests/unit/test_input.py @@ -1,6 +1,7 @@ import pytest from process import fortran from process.utilities.f2py_string_patch import string_to_f2py_compatible +import process.init as init def _create_input_file(directory, content: str): @@ -57,6 +58,6 @@ def test_parse_real(epsvmc, expected, tmp_path): fortran.global_variables.fileprefix, _create_input_file(tmp_path, f"epsvmc = {epsvmc}"), ) - fortran.init_module.init() + init.init_process() assert fortran.numerics.epsvmc.item() == expected From 1ffcea9b7e438f5aedae5c35d9fcfd364b658d89 Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Thu, 14 Nov 2024 17:16:30 +0000 Subject: [PATCH 2/3] Add new error classes for PROCESS-raised errors --- process/exceptions.py | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) create mode 100644 process/exceptions.py diff --git a/process/exceptions.py b/process/exceptions.py new file mode 100644 index 0000000000..b7d945d11c --- /dev/null +++ b/process/exceptions.py @@ -0,0 +1,19 @@ +class ProcessException(Exception): + def __init__(self, *args, **kwargs): + super().__init__(*args) + self._diagnostics = kwargs + + def __str__(self): + exception_message = super().__str__() + diagnostics_message = "\n".join( + [f"\t{d}: {repr(v)}" for d, v in self._diagnostics.items()] + ) + + if diagnostics_message: + return f"{exception_message}\n{diagnostics_message}" + + return exception_message + + +class ProcessValueError(ProcessException, ValueError): + pass From 65e198e31852e9d9c4b8c0d1cf4784f55f6724bd Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Thu, 21 Nov 2024 17:08:47 +0000 Subject: [PATCH 3/3] Convert process input checking routines to Python --- process/exceptions.py | 10 + process/init.py | 842 ++++++++++++++++++++++++++++++++- source/fortran/init_module.f90 | 752 ----------------------------- 3 files changed, 851 insertions(+), 753 deletions(-) diff --git a/process/exceptions.py b/process/exceptions.py index b7d945d11c..f087bb129d 100644 --- a/process/exceptions.py +++ b/process/exceptions.py @@ -1,4 +1,6 @@ class ProcessException(Exception): + """A base Exception to derive other PROCESS exceptions from""" + def __init__(self, *args, **kwargs): super().__init__(*args) self._diagnostics = kwargs @@ -15,5 +17,13 @@ def __str__(self): return exception_message +class ProcessValidationError(ProcessException): + """Exception raised when validating PROCESS input. + + E.g. initial values, constraint/variable combinations, switch combinations""" + + pass + + class ProcessValueError(ProcessException, ValueError): pass diff --git a/process/init.py b/process/init.py index f58a51c793..3460828187 100644 --- a/process/init.py +++ b/process/init.py @@ -1,4 +1,6 @@ +from warnings import warn import process.fortran as fortran +from process.exceptions import ProcessValidationError def init_process(): @@ -26,7 +28,7 @@ def init_process(): fortran.stellarator_module.stinit() # Check input data for errors/ambiguities - fortran.init_module.check() + check_process() fortran.main_module.run_summary() @@ -250,3 +252,841 @@ def initialise_iterative_variables(): fortran.define_iteration_variables.init_itv_173() fortran.define_iteration_variables.init_itv_174() fortran.define_iteration_variables.init_itv_175() + + +def check_process(): + """Routine to reset specific variables if certain options are + being used + author: P J Knight, CCFE, Culham Science Centre + None + This routine performs a sanity check of the input variables + and ensures other dependent variables are given suitable values. + """ + # error_handling.errors_on = True + + # Check that there are sufficient iteration variables + if fortran.numerics.nvar < fortran.numerics.neqns: + raise ProcessValidationError( + "Insufficient iteration variables to solve the problem! NVAR < NEQNS", + nvar=fortran.numerics.nvar, + neqns=fortran.numerics.neqns, + ) + + # Check that sufficient elements of ixc and icc have been specified + if (fortran.numerics.ixc[: fortran.numerics.nvar] == 0).any(): + raise ProcessValidationError( + "The number of iteration variables specified is smaller than the number stated in ixc", + nvar=fortran.numerics.nvar, + ) + + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] == 0 + ).any(): + raise ProcessValidationError( + "The number of constraints specified is smaller than the number stated in neqns+nineqns", + neqns=fortran.numerics.neqns, + nineqns=fortran.numerics.nineqns, + ) + + # Deprecate constraints + for depcrecated_constraint in [3, 4, 10, 74, 42]: + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] + == depcrecated_constraint + ).any(): + raise ProcessValidationError( + "Constraint equation is no longer available", icc=depcrecated_constraint + ) + + # MDK Report error if constraint 63 is used with old vacuum model + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] == 63 + ).any() and fortran.vacuum_variables.vacuum_model != "simple": + raise ProcessValidationError( + "Constraint 63 is requested without the correct vacuum model (simple)" + ) + + # Fuel ion fractions must add up to 1.0 + if ( + abs( + 1.0 + - fortran.physics_variables.f_deuterium + - fortran.physics_variables.f_tritium + - fortran.physics_variables.f_helium3 + ) + > 1e-6 + ): + raise ProcessValidationError( + "Fuel ion fractions do not sum to 1.0", + f_deuterium=fortran.physics_variables.f_deuterium, + f_tritium=fortran.physics_variables.f_tritium, + f_helium3=fortran.physics_variables.f_helium3, + ) + + if fortran.physics_variables.f_tritium < 1.0e-3: # tritium fraction is negligible + fortran.buildings_variables.triv = 0.0 + fortran.heat_transport_variables.trithtmw = 0.0 + + if fortran.impurity_radiation_module.fimp[1] != 0.1: + raise ProcessValidationError( + "The thermal alpha/electron density ratio should be controlled using ralpne (itv 109) and not fimp(2)." + "fimp(2) should be removed from the input file, or set to the default value 0.1D0." + ) + + # Impurity fractions + for imp in range(fortran.impurity_radiation_module.nimp): + fortran.impurity_radiation_module.impurity_arr_frac[ + imp + ] = fortran.impurity_radiation_module.fimp[imp] + + # Stop the run if oacdcp is used as an optimisation variable + # As the current density is now calculated from bt without constraint 10 + + if (fortran.numerics.ixc[: fortran.numerics.nvar] == 12).any(): + raise ProcessValidationError( + "The 1/R toroidal B field dependency constraint is being depreciated" + ) + + # Plasma profile consistency checks + if fortran.ife_variables.ife != 1: + if fortran.physics_variables.ipedestal == 1: + + # Temperature checks + if fortran.physics_variables.teped < fortran.physics_variables.tesep: + raise ProcessValidationError( + "Pedestal temperature is lower than separatrix temperature", + teped=fortran.physics_variables.teped, + tesep=fortran.physics_variables.tesep, + ) + + if (abs(fortran.physics_variables.rhopedt - 1.0) <= 1e-7) and ( + (fortran.physics_variables.teped - fortran.physics_variables.tesep) + >= 1e-7 + ): + warn( + f"Temperature pedestal is at plasma edge, but teped " + f"({fortran.physics_variables.teped}) differs from tesep " + f"({fortran.physics_variables.tesep})" + ) + + # Core temperature should always be calculated (later) as being + # higher than the pedestal temperature, if and only if the + # volume-averaged temperature never drops below the pedestal + # temperature. Prevent this by adjusting te, and its lower bound + # (which will only have an effect if this is an optimisation run) + if fortran.physics_variables.te <= fortran.physics_variables.teped: + warn( + f"Volume-averaged temperature ({fortran.physics_variables.te}) has been " + f"forced to exceed input pedestal height ({fortran.physics_variables.teped}). " + "Changing to te = teped*1.001" + ) + fortran.physics_variables.te = fortran.physics_variables.teped * 1.001 + + if ( + fortran.numerics.ioptimz >= 0 + and (fortran.numerics.ixc[: fortran.numerics.nvar] == 4).any() + and fortran.numerics.boundl[3] < fortran.physics_variables.teped * 1.001 + ): + warn( + "Lower limit of volume averaged electron temperature (te) has been raised to ensure te > teped" + ) + fortran.numerics.boundl[3] = fortran.physics_variables.teped * 1.001 + fortran.numerics.boundu[3] = max( + fortran.numerics.boundu[3], fortran.numerics.boundl[3] + ) + + # Density checks + # Case where pedestal density is set manually + if ( + fortran.physics_variables.fgwped < 0 + or not (fortran.numerics.ixc[: fortran.numerics.nvar] == 145).any() + ): + + # Issue #589 Pedestal density is set manually using neped but it is less than nesep. + if fortran.physics_variables.neped < fortran.physics_variables.nesep: + raise ProcessValidationError( + "Density pedestal is lower than separatrix density", + neped=fortran.physics_variables.neped, + nesep=fortran.physics_variables.nesep, + ) + + # Issue #589 Pedestal density is set manually using neped, + # but pedestal width = 0. + if ( + abs(fortran.physics_variables.rhopedn - 1.0) <= 1e-7 + and ( + fortran.physics_variables.neped + - fortran.physics_variables.nesep + ) + >= 1e-7 + ): + warn( + "Density pedestal is at plasma edge " + f"({fortran.physics_variables.rhopedn = }), but neped " + f"({fortran.physics_variables.neped}) differs from " + f"nesep ({fortran.physics_variables.nesep})" + ) + + # Issue #862 : Variable ne0/neped ratio without constraint eq 81 (ne0>neped) + # -> Potential hollowed density profile + if ( + fortran.numerics.ioptimz >= 0 + and not ( + fortran.numerics.icc[ + : fortran.numerics.neqns + fortran.numerics.nineqns + ] + == 81 + ).any() + ): + if (fortran.numerics.ixc[: fortran.numerics.nvar] == 145).any(): + warn("neped set with fgwped without constraint eq 81 (neped 273.15 K + if fortran.tfcoil_variables.tcoolin < 273.15: + raise ProcessValidationError( + "Coolant temperature (tcoolin) cannot be < 0 C (273.15 K) for water cooled copper magents" + ) + + # Temperature of the TF legs cannot be cooled down + if ( + fortran.tfcoil_variables.tlegav > 0 + and fortran.tfcoil_variables.tlegav < 273.15 + ): + raise ProcessValidationError( + "TF legs conductor temperature (tlegav) cannot be < 0 C (273.15 K) for water cooled magents" + ) + + # Check if conductor upper limit is properly set to 50 K or below + if ( + fortran.numerics.ixc[: fortran.numerics.nvar] == 20 + ).any() and fortran.numerics.boundu[19] < 273.15: + raise ProcessValidationError( + "Too low CP conductor temperature (tcpav). Lower limit for copper > 273.15 K" + ) + + # Call a lvl 3 error if superconductor magnets are used + elif fortran.tfcoil_variables.i_tf_sup == 1: + warn( + "Joints res not cal. for SC (itart = 1) TF (fortran.tfcoil_variables.i_tf_sup = 1)" + ) + + # Aluminium magnets initalisation / checks + # Initialize the CP conductor temperature to cryogenic temperature for cryo-al magnets (20 K) + elif fortran.tfcoil_variables.i_tf_sup == 2: + + # Call a lvl 3 error if the inlet coolant temperature is too large + # Motivation : ill-defined aluminium resistivity fit for T > 40-50 K + if fortran.tfcoil_variables.tcoolin > 40.0: + raise ProcessValidationError( + "Coolant temperature (tcoolin) should be < 40 K for the cryo-al resistivity to be defined" + ) + + # Check if the leg average temperature is low enough for the resisitivity fit + if fortran.tfcoil_variables.tlegav > 50.0: + raise ProcessValidationError( + "TF legs conductor temperature (tlegav) should be < 40 K for the cryo-al resistivity to be defined" + ) + + # Check if conductor upper limit is properly set to 50 K or below + if ( + fortran.numerics.ixc[: fortran.numerics.nvar] == 20 + ).any() and fortran.numerics.boundu[19] > 50.0: + raise ProcessValidationError( + "Too large CP conductor temperature (tcpav). Upper limit for cryo-al < 50 K" + ) + + # Otherwise intitialise the average conductor temperature at + fortran.tfcoil_variables.tcpav = fortran.tfcoil_variables.tcoolin + + # Check if the boostrap current selection is addapted to ST + if fortran.physics_variables.i_bootstrap_current == 1: + raise ProcessValidationError( + "Invalid boostrap current law for ST, do not use i_bootstrap_current = 1" + ) + + # Check if a single null divertor is used in double null machine + if fortran.physics_variables.i_single_null == 0 and ( + fortran.physics_variables.ftar == 1.0 + or fortran.physics_variables.ftar == 0.0 + ): + warn("Operating with a single null in a double null machine") + + # Set the TF coil shape to picture frame (if default value) + if fortran.tfcoil_variables.i_tf_shape == 0: + fortran.tfcoil_variables.i_tf_shape = 2 + + # Warning stating that the CP fast neutron fluence calculation + # is not addapted for cryoaluminium calculations yet + if ( + fortran.tfcoil_variables.i_tf_sup == 2 + and ( + fortran.numerics.icc[ + : fortran.numerics.neqns + fortran.numerics.nineqns + ] + == 85 + ).any() + and fortran.physics_variables.itart == 1 + ): + raise ProcessValidationError( + "Al TF coil fluence not calculated properly for Al CP, do not use constraint 85" + ) + + # Setting the CP joints default options : + # 0 : No joints for superconducting magents (fortran.tfcoil_variables.i_tf_sup = 1) + # 1 : Sliding joints for resistive magnets (fortran.tfcoil_variables.i_tf_sup = 0, 2) + if fortran.tfcoil_variables.i_cp_joints == -1: + if fortran.tfcoil_variables.i_tf_sup == 1: + fortran.tfcoil_variables.i_cp_joints = 0 + else: + fortran.tfcoil_variables.i_cp_joints = 1 + + # Checking the CP TF top radius + if ( + abs(fortran.build_variables.r_cp_top) > 0 + or (fortran.numerics.ixc[: fortran.numerics.nvar] == 174).any() + ) and fortran.build_variables.i_r_cp_top != 1: + raise ProcessValidationError( + "To set the TF CP top value, you must use i_r_cp_top = 1" + ) + + # Conventionnal aspect ratios specific + else: + + if ( + fortran.physics_variables.i_plasma_current == 2 + or fortran.physics_variables.i_plasma_current == 9 + ): + raise ProcessValidationError( + "i_plasma_current=2,9 is not a valid option for a non-TART device" + ) + + # Set the TF coil shape to PROCESS D-shape (if default value) + if fortran.tfcoil_variables.i_tf_shape == 0: + fortran.tfcoil_variables.i_tf_shape = 1 + + # Check PF coil configurations + j = 0 + k = 0 + for i in range(fortran.pfcoil_variables.ngrp): + if ( + fortran.pfcoil_variables.ipfloc[i] != 2 + and fortran.pfcoil_variables.ncls[i] != 2 + ): + raise ProcessValidationError( + "ncls(i) .ne. 2 is not a valid option except for (ipfloc = 2)" + ) + + if fortran.pfcoil_variables.ipfloc[i] == 2: + j = j + 1 + k = k + fortran.pfcoil_variables.ncls[i] + + if k == 1: + raise ProcessValidationError( + "Only 1 divertor coil (ipfloc = 2) is not a valid configuration" + ) + if k > 2: + raise ProcessValidationError( + "More than 2 divertor coils (ipfloc = 2) is not a valid configuration" + ) + if fortran.physics_variables.i_single_null == 1 and j < 2: + raise ProcessValidationError( + "If snull=1, use 2 individual divertor coils (ipfloc = 2, 2; ncls = 1, 1)" + ) + + # Constraint 10 is dedicated to ST designs with demountable joints + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] + == 10 + ).any(): + raise ProcessValidationError( + "Constraint equation 10 (CP lifetime) to used with ST desing (itart=1)" + ) + + # Pulsed power plant model + if fortran.pulse_variables.lpulse == 1: + fortran.global_variables.icase = "Pulsed tokamak model" + else: + fortran.buildings_variables.esbldgm3 = 0.0 + + # TF coil + # ------- + # TF stress model not defined of r_tf_inboard = 0 + # Unless i_tf_stress_model == 2 + # -> If bore + gapoh + ohcth = 0 and fixed and stress constraint is used + # Generate a lvl 3 error proposing not to use any stress constraints + if ( + ( + not ( + (fortran.numerics.ixc[: fortran.numerics.nvar] == 16).any() + or (fortran.numerics.ixc[: fortran.numerics.nvar] == 29).any() + or (fortran.numerics.ixc[: fortran.numerics.nvar] == 42).any() + ) + ) # No bore,gapoh, ohcth iteration + and ( + abs( + fortran.build_variables.bore + + fortran.build_variables.gapoh + + fortran.build_variables.ohcth + + fortran.build_variables.precomp + ) + <= 0 + ) # bore + gapoh + ohcth = 0 + and ( + ( + fortran.numerics.icc[ + : fortran.numerics.neqns + fortran.numerics.nineqns + ] + == 31 + ).any() + or ( + fortran.numerics.icc[ + : fortran.numerics.neqns + fortran.numerics.nineqns + ] + == 32 + ).any() + ) # Stress constraints (31 or 32) is used + and ( + fortran.tfcoil_variables.i_tf_stress_model != 2 + ) # TF stress model can't handle no bore + ): + raise ProcessValidationError( + "Invalid stress model if bore + gapoh + ohcth = 0. Don't use constraint 31" + ) + + # Make sure that plane stress model is not used for resistive magnets + if ( + fortran.tfcoil_variables.i_tf_stress_model == 1 + and fortran.tfcoil_variables.i_tf_sup != 1 + ): + raise ProcessValidationError( + "Use generalized plane strain for resistive magnets (i_tf_stress_model = 0 or 2)" + ) + + # bucking cylinder default option setting + # - bucking (casing) for SC i_tf_bucking ( i_tf_bucking = 1 ) + # - No bucking for copper magnets ( i_tf_bucking = 0 ) + # - Bucking for aluminium magnets ( i_tf_bucking = 1 ) + if fortran.tfcoil_variables.i_tf_bucking == -1: + if fortran.tfcoil_variables.i_tf_sup == 0: + fortran.tfcoil_variables.i_tf_bucking = 0 + else: + fortran.tfcoil_variables.i_tf_bucking = 1 + + # Ensure that the TF isnt placed against the + # CS which is now outside it + if ( + fortran.tfcoil_variables.i_tf_bucking >= 2 + and fortran.build_variables.tf_in_cs == 1 + ): + raise ProcessValidationError("Cannot have i_tf_bucking >= 2 when tf_in_cs = 1") + + # Ensure that no pre-compression structure + # is used for bucked and wedged design + if ( + fortran.tfcoil_variables.i_tf_bucking >= 2 + and fortran.build_variables.iprecomp == 1 + ): + raise ProcessValidationError( + "No CS precompression structure for bucked and wedged, use iprecomp = 0" + ) + + # Number of stress calculation layers + # +1 to add in the inboard TF coil case on the plasma side, per Issue #1509 + fortran.tfcoil_variables.n_tf_stress_layers = ( + fortran.tfcoil_variables.i_tf_bucking + + fortran.tfcoil_variables.n_tf_graded_layers + + 1 + ) + + # If TFC sidewall has not been set by user + if fortran.tfcoil_variables.casths < 0.1e-10: + fortran.tfcoil_variables.tfc_sidewall_is_fraction = True + + # If inboard TF coil case plasma side thickness has not been set by user + if fortran.tfcoil_variables.casthi < 0.1e-10: + fortran.tfcoil_variables.casthi_is_fraction = True + + # Setting the default cryo-plants efficiencies + if abs(fortran.tfcoil_variables.eff_tf_cryo + 1) < 1e-6: + + # The ITER cyoplant efficiency is used for SC + if fortran.tfcoil_variables.i_tf_sup == 1: + fortran.tfcoil_variables.eff_tf_cryo = 0.13 + + # Strawbrige plot extrapolation is used for Cryo-Al + elif fortran.tfcoil_variables.i_tf_sup == 2: + fortran.tfcoil_variables.eff_tf_cryo = 0.40 + + # Cryo-plane efficiency must be in [0-1.0] + elif ( + fortran.tfcoil_variables.eff_tf_cryo > 1.0 + or fortran.tfcoil_variables.eff_tf_cryo < 0.0 + ): + raise ProcessValidationError( + "TF cryo-plant efficiency `eff_tf_cryo` must be within [0-1]" + ) + + # Integer turns option not yet available for REBCO taped turns + + if ( + fortran.tfcoil_variables.i_tf_sc_mat == 6 + and fortran.tfcoil_variables.i_tf_turns_integer == 1 + ): + raise ProcessValidationError( + "Integer turns (i_tf_turns_integer = 1) not supported for REBCO (i_tf_sc_mat = 6)" + ) + + # Setting up insulation layer young modulae default values [Pa] + + if fortran.tfcoil_variables.eyoung_ins <= 1.0e8: + + # Copper magnets, no insulation material defined + # But use the ITER design by default + if fortran.tfcoil_variables.i_tf_sup == 0: + fortran.tfcoil_variables.eyoung_ins = 20.0e9 + + # SC magnets + # Value from DDD11-2 v2 2 (2009) + elif fortran.tfcoil_variables.i_tf_sup == 1: + fortran.tfcoil_variables.eyoung_ins = 20.0e9 + + # Cryo-aluminum magnets (Kapton polymer) + elif fortran.tfcoil_variables.i_tf_sup == 2: + fortran.tfcoil_variables.eyoung_ins = 2.5e9 + + # Setting the default WP geometry + + if fortran.tfcoil_variables.i_tf_wp_geom == -1: + if fortran.tfcoil_variables.i_tf_turns_integer == 0: + fortran.tfcoil_variables.i_tf_wp_geom = 1 + if fortran.tfcoil_variables.i_tf_turns_integer == 1: + fortran.tfcoil_variables.i_tf_wp_geom = 0 + + # Setting the TF coil conductor elastic properties + + if fortran.tfcoil_variables.i_tf_cond_eyoung_axial == 0: + # Conductor stiffness is not considered + fortran.tfcoil_variables.eyoung_cond_axial = 0 + fortran.tfcoil_variables.eyoung_cond_trans = 0 + elif fortran.tfcoil_variables.i_tf_cond_eyoung_axial == 2: + # Select sensible defaults from the literature + if fortran.tfcoil_variables.i_tf_sc_mat in [1, 4, 5]: + # Nb3Sn: Nyilas, A et. al, Superconductor Science and Technology 16, no. 9 (2003): 1036–42. https://doi.org/10.1088/0953-2048/16/9/313. + fortran.tfcoil_variables.eyoung_cond_axial = 32e9 + elif fortran.tfcoil_variables.i_tf_sc_mat == 2: + # Bi-2212: Brown, M. et al, IOP Conference Series: Materials Science and Engineering 279 (2017): 012022. https://doi.org/10.1088/1757-899X/279/1/012022. + fortran.tfcoil_variables.eyoung_cond_axial = 80e9 + elif fortran.tfcoil_variables.i_tf_sc_mat in [3, 7]: + # NbTi: Vedrine, P. et. al, IEEE Transactions on Applied Superconductivity 9, no. 2 (1999): 236–39. https://doi.org/10.1109/77.783280. + fortran.tfcoil_variables.eyoung_cond_axial = 6.8e9 + elif fortran.tfcoil_variables.i_tf_sc_mat in [6, 8, 9]: + # REBCO: Fujishiro, H. et. al, Physica C: Superconductivity, 426–431 (2005): 699–704. https://doi.org/10.1016/j.physc.2005.01.045. + fortran.tfcoil_variables.eyoung_cond_axial = 145e9 + + if fortran.tfcoil_variables.i_tf_cond_eyoung_trans == 0: + # Transverse stiffness is not considered + fortran.tfcoil_variables.eyoung_cond_trans = 0 + else: + # Transverse stiffness is significant + fortran.tfcoil_variables.eyoung_cond_trans = ( + fortran.tfcoil_variables.eyoung_cond_axial + ) + + # Check if the WP/conductor radial thickness (dr_tf_wp) is large enough + # To contains the insulation, cooling and the support structure + # Rem : Only verified if the WP thickness is used + if (fortran.numerics.ixc[: fortran.numerics.nvar] == 140).any(): + + # Minimal WP thickness + if fortran.tfcoil_variables.i_tf_sup == 1: + dr_tf_wp_min = 2.0 * ( + fortran.tfcoil_variables.tinstf + + fortran.tfcoil_variables.tfinsgap + + fortran.tfcoil_variables.thicndut + + fortran.tfcoil_variables.dhecoil + ) + + # Steel conduit thickness (can be an iteration variable) + if (fortran.numerics.ixc[: fortran.numerics.nvar] == 58).any(): + dr_tf_wp_min = dr_tf_wp_min + 2.0 * fortran.numerics.boundl[57] + else: + dr_tf_wp_min = dr_tf_wp_min + 2.0 * fortran.tfcoil_variables.thwcndut + + # Minimal conductor layer thickness + elif ( + fortran.tfcoil_variables.i_tf_sup == 0 + or fortran.tfcoil_variables.i_tf_sup == 2 + ): + dr_tf_wp_min = ( + 2.0 + * (fortran.tfcoil_variables.thicndut + fortran.tfcoil_variables.tinstf) + + 4.0 * fortran.tfcoil_variables.rcool + ) + + if fortran.numerics.boundl[139] < dr_tf_wp_min: + raise ProcessValidationError( + "The TF coil WP thickness (dr_tf_wp) must be at least", + dr_tf_wp_min=dr_tf_wp_min, + ) + + # Setting t_turn_tf_is_input to true if t_turn_tf is an input + fortran.tfcoil_variables.t_turn_tf_is_input = ( + abs(fortran.tfcoil_variables.t_turn_tf) > 0 + ) + + # Impossible to set the turn size of integer turn option + if ( + fortran.tfcoil_variables.t_turn_tf_is_input + and fortran.tfcoil_variables.i_tf_turns_integer == 1 + ): + raise ProcessValidationError( + "Impossible to set the TF turn/cable size with the integer turn option (i_tf_turns_integer: 1)" + ) + + if ( + fortran.tfcoil_variables.i_tf_wp_geom != 0 + and fortran.tfcoil_variables.i_tf_turns_integer == 1 + ): + raise ProcessValidationError( + "Can only have i_tf_turns_integer = 1 with i_tf_wp_geom = 0" + ) + + if ( + fortran.physics_variables.i_bootstrap_current == 5 + and fortran.physics_variables.i_diamagnetic_current != 0 + ): + raise ProcessValidationError( + "i_diamagnetic_current = 0 should be used with the Sakai plasma current scaling" + ) + + # Setting t_cable_tf_is_input to true if t_cable_tf is an input + fortran.tfcoil_variables.t_cable_tf_is_input = ( + abs(fortran.tfcoil_variables.t_cable_tf) > 0 + ) + + # Impossible to set the cable size of integer turn option + if ( + fortran.tfcoil_variables.t_cable_tf_is_input + and fortran.tfcoil_variables.i_tf_turns_integer == 1 + ): + raise ProcessValidationError( + "Impossible to set the TF turn/cable size with the integer turn option (i_tf_turns_integer: 1)" + ) + + # Impossible to set both the TF coil turn and the cable dimension + if ( + fortran.tfcoil_variables.t_turn_tf_is_input + and fortran.tfcoil_variables.t_cable_tf_is_input + ): + raise ProcessValidationError( + "Impossible to set the TF coil turn and cable size simultaneously" + ) + + # Checking the SC temperature for LTS + if ( + fortran.tfcoil_variables.i_tf_sc_mat in [1, 3, 4, 5] + and fortran.tfcoil_variables.tftmp > 10.0 + ): + raise ProcessValidationError( + "The LTS conductor temperature (tftmp) has to be lower than 10" + ) + + # PF coil resistivity is zero if superconducting + if fortran.pfcoil_variables.ipfres == 0: + fortran.pfcoil_variables.pfclres = 0.0 + + # If there is no NBI, then hot beam density should be zero + if fortran.current_drive_variables.irfcd == 1: + if ( + fortran.current_drive_variables.iefrf != 5 + and fortran.current_drive_variables.iefrf != 8 + ): + fortran.physics_variables.rnbeam = 0.0 + else: + fortran.physics_variables.rnbeam = 0.0 + + # Set inboard blanket thickness to zero if no inboard blanket switch + # used (Issue #732) + if fortran.fwbs_variables.iblnkith == 0: + fortran.build_variables.blnkith = 0.0 + + # Ensure that blanket material fractions allow non-zero space for steel + # CCFE HCPB Model + + if fortran.stellarator_variables.istell == 0: + if fortran.fwbs_variables.iblanket == 1 or fortran.fwbs_variables.iblanket == 3: + fsum = ( + fortran.fwbs_variables.breeder_multiplier + + fortran.fwbs_variables.vfcblkt + + fortran.fwbs_variables.vfpblkt + ) + if fsum >= 1.0: + raise ProcessValidationError( + "Blanket material fractions do not sum to 1.0", + iblanket=fortran.fwbs_variables.iblanket, + breeder_multiplier=fortran.fwbs_variables.breeder_multiplier, + vfcblkt=fortran.fwbs_variables.vfcblkt, + vfpblkt=fortran.fwbs_variables.vfpblkt, + fsum=fsum, + ) + + # Initialise superconductor cable parameters + if fortran.tfcoil_variables.i_tf_sup == 1: + fortran.sctfcoil_module.initialise_cables() + + # Check that the temperature margins are not overdetermined + if fortran.tfcoil_variables.tmargmin > 0.0001: + # This limit has been input and will be applied to both TFC and CS + if fortran.tfcoil_variables.tmargmin_tf > 0.0001: + warn( + "tmargmin_tf and tmargmin should not both be specified in IN.DAT " + "tmargmin_tf has been ignored" + ) + if fortran.tfcoil_variables.tmargmin_cs > 0.0001: + warn( + "tmargmin_cs and tmargmin should not both be specified in IN.DAT " + "tmargmin_cs has been ignored" + ) + + fortran.tfcoil_variables.tmargmin_tf = fortran.tfcoil_variables.tmargmin + fortran.tfcoil_variables.tmargmin_cs = fortran.tfcoil_variables.tmargmin + + if ( + fortran.physics_variables.tauee_in > 1e-10 + and fortran.physics_variables.isc != 48 + ): + # Report error if confinement time is in the input + # but the scaling to use it is not selected. + warn("tauee_in is for use with isc=48 only") + + if fortran.physics_variables.aspect > 1.7 and fortran.physics_variables.isc == 46: + # NSTX scaling is for A<1.7 + warn("NSTX scaling is for A<1.7") + + if ( + fortran.physics_variables.i_plasma_current == 2 + and fortran.physics_variables.isc == 42 + ): + raise ProcessValidationError( + "Lang 2012 confinement scaling cannot be used for i_plasma_current=2 due to wrong q" + ) + + # Cannot use temperature margin constraint with REBCO TF coils + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] == 36 + ).any() and ( + fortran.tfcoil_variables.i_tf_sc_mat == 8 + or fortran.tfcoil_variables.i_tf_sc_mat == 9 + ): + raise ProcessValidationError( + "turn off TF temperature margin constraint icc = 36 when using REBCO" + ) + + # Cannot use temperature margin constraint with REBCO CS coils + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] == 60 + ).any() and fortran.pfcoil_variables.isumatoh == 8: + raise ProcessValidationError( + "turn off CS temperature margin constraint icc = 60 when using REBCO" + ) + + # Cold end of the cryocooler should be colder than the TF + if fortran.tfcoil_variables.tmpcry > fortran.tfcoil_variables.tftmp: + raise ProcessValidationError("tmpcry should be lower than tftmp") + + # Cannot use TF coil strain limit if i_str_wp is off: + if ( + fortran.numerics.icc[: fortran.numerics.neqns + fortran.numerics.nineqns] == 88 + ).any() and fortran.tfcoil_variables.i_str_wp == 0: + raise ProcessValidationError("Can't use constraint 88 if i_strain_tf == 0") diff --git a/source/fortran/init_module.f90 b/source/fortran/init_module.f90 index ee4d42462a..d2b00d683f 100644 --- a/source/fortran/init_module.f90 +++ b/source/fortran/init_module.f90 @@ -96,756 +96,4 @@ subroutine finish close(unit = opt_file) if (verbose == 1) close(unit = vfile) end subroutine finish - - subroutine check - - !! Routine to reset specific variables if certain options are - !! being used - !! author: P J Knight, CCFE, Culham Science Centre - !! None - !! This routine performs a sanity check of the input variables - !! and ensures other dependent variables are given suitable values. - - !! ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use build_variables, only: blnkith, bore, gapoh, ohcth, precomp, iprecomp, & - i_r_cp_top, r_cp_top, vgaptop, vgap_xpoint_divertor, shldtth, shldlth, d_vv_top, d_vv_bot, tf_in_cs - use buildings_variables, only: esbldgm3, triv - use current_drive_variables, only: gamcd, iefrf, irfcd - use error_handling, only: errors_on, idiags, fdiags, report_error - use fwbs_variables, only: breeder_multiplier, iblanket, vfcblkt, vfpblkt, & - iblnkith - use global_variables, only: icase - use heat_transport_variables, only: trithtmw - use ife_variables, only: ife - use impurity_radiation_module, only: nimp, impurity_arr_frac, fimp - use numerics, only: ixc, icc, ioptimz, neqns, nineqns, nvar, boundl, & - boundu - use pfcoil_variables, only: ipfres, ngrp, pfclres, ipfloc, ncls, isumatoh - use physics_variables, only: aspect, f_deuterium, fgwped, f_helium3, & - fgwsep, f_tritium, i_bootstrap_current, i_single_null, i_plasma_current, idivrt, ishape, & - iradloss, isc, ipedestal, ilhthresh, itart, nesep, rhopedn, rhopedt, & - rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar, i_diamagnetic_current - use pulse_variables, only: lpulse - use reinke_variables, only: fzactual, impvardiv - use tfcoil_variables, only: casthi, casthi_is_fraction, casths, i_tf_sup, & - tcoolin, tcpav, tfc_sidewall_is_fraction, tmargmin, tmargmin_cs, & - tmargmin_tf, eff_tf_cryo, eyoung_ins, i_tf_bucking, i_tf_shape, & - n_tf_graded_layers, n_tf_stress_layers, tlegav, i_tf_stress_model, & - i_tf_sc_mat, i_tf_wp_geom, i_tf_turns_integer, tinstf, thwcndut, & - tfinsgap, rcool, dhecoil, thicndut, i_cp_joints, t_turn_tf_is_input, & - t_turn_tf, tftmp, t_cable_tf, t_cable_tf_is_input, tftmp, tmpcry, & - i_tf_cond_eyoung_axial, eyoung_cond_axial, eyoung_cond_trans, & - i_tf_cond_eyoung_trans, i_str_wp - use stellarator_variables, only: istell - use sctfcoil_module, only: initialise_cables - use vacuum_variables, only: vacuum_model - use, intrinsic :: iso_fortran_env, only: dp=>real64 - - implicit none - - ! Local variables - - integer :: i,j,k,imp - real(dp) :: fsum - - real(dp) :: dr_tf_wp_min - !! Minimal WP or conductor layer thickness [m] - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - errors_on = .true. - - ! Check that there are sufficient iteration variables - if (nvar < neqns) then - idiags(1) = nvar ; idiags(2) = neqns - call report_error(137) - end if - - ! Check that sufficient elements of ixc and icc have been specified - if ( any(ixc(1:nvar) == 0) ) then - idiags(1) = nvar - call report_error(139) - end if - - - if ( any(icc(1:neqns+nineqns) == 0) ) then - idiags(1) = neqns ; idiags(2) = nineqns - call report_error(140) - end if - - ! Deprecate constraints 3 and 4 - if ( any(icc(1:neqns+nineqns) == 3) ) then - call report_error(162) - write(*,*) 'PROCESS stopping' - stop 1 - end if - - if ( any(icc(1:neqns+nineqns) == 4) ) then - call report_error(163) - write(*,*) 'PROCESS stopping' - stop 1 - end if - - - ! MDK Report error if constraint 63 is used with old vacuum model - if (any(icc(1:neqns+nineqns) == 63).and.(vacuum_model.ne.'simple') ) then - write(*,*) 'Constraint 63 is requested without the correct vacuum model ("simple").' - write(*,*) 'vacuum_model = ', vacuum_model - write(*,*) 'PROCESS stopping' - stop 1 - end if - - if ( any(icc(1:neqns+nineqns) == 74) ) then - write(*,*)'Constraint 74 (TF coil quench temperature for Croco HTS conductor) is not yet implemented' - write(*,*) 'PROCESS stopping' - stop 1 - end if - - ! Fuel ion fractions must add up to 1.0 - if (abs(1.0D0 - f_deuterium - f_tritium - f_helium3) > 1.0D-6) then - fdiags(1) = f_deuterium; fdiags(2) = f_tritium ; fdiags(3) = f_helium3 - call report_error(36) - end if - - if (f_tritium < 1.0D-3) then ! tritium fraction is negligible - triv = 0.0D0 - trithtmw = 0.0D0 - end if - - if (fimp(2) .ne. 0.1D0) then - write(*,*)'The thermal alpha/electron density ratio should be controlled using ralpne (itv 109) and not fimp(2).' - write(*,*)'fimp(2) should be removed from the input file, or set to the default value 0.1D0.' - stop 1 - end if - - ! Impurity fractions - do imp = 1,nimp - impurity_arr_frac(imp) = fimp(imp) - end do - - ! The 1/R B field dependency constraint variable is being depreciated - ! Stop the run if the constraint 10 is used - if ( any( icc == 10 ) ) then - call report_error(236) - stop 1 - end if - - ! Stop the run if oacdcp is used as an optimisation variable - ! As the current density is now calculated from bt without constraint 10 - if ( any( ixc == 12 ) ) then - call report_error(236) - stop 1 - end if - - ! Warn if ion power balance equation is being used with the new radiation model - if (any(icc == 3)) then - call report_error(138) - end if - - ! Plasma profile consistency checks - if (ife /= 1) then - if (ipedestal == 1) then - - ! Temperature checks - if (teped < tesep) then - fdiags(1) = teped ; fdiags(2) = tesep - call report_error(146) - end if - - if ((abs(rhopedt-1.0D0) <= 1.0D-7).and.((teped-tesep) >= 1.0D-7)) then - fdiags(1) = rhopedt ; fdiags(2) = teped ; fdiags(3) = tesep - call report_error(147) - end if - - ! Core temperature should always be calculated (later) as being - ! higher than the pedestal temperature, if and only if the - ! volume-averaged temperature never drops below the pedestal - ! temperature. Prevent this by adjusting te, and its lower bound - ! (which will only have an effect if this is an optimisation run) - if (te <= teped) then - fdiags(1) = te ; fdiags(2) = teped - te = teped*1.001D0 - call report_error(149) - end if - - if ((ioptimz >= 0).and.(any(ixc == 4)).and.(boundl(4) < teped*1.001D0)) then - call report_error(150) - boundl(4) = teped*1.001D0 - boundu(4) = max(boundu(4), boundl(4)) - end if - - ! Density checks - ! Case where pedestal density is set manually - ! --------------- - if ( (fgwped < 0) .or. (.not.any(ixc==145)) ) then - - ! Issue #589 Pedestal density is set manually using neped but it is less than nesep. - if ( neped < nesep ) then - fdiags(1) = neped ; fdiags(2) = nesep - call report_error(151) - end if - - ! Issue #589 Pedestal density is set manually using neped, - ! but pedestal width = 0. - if ( (abs(rhopedn-1.0D0) <= 1.0D-7).and.((neped-nesep) >= 1.0D-7) ) then - fdiags(1) = rhopedn ; fdiags(2) = neped ; fdiags(3) = nesep - call report_error(152) - end if - end if - - ! Issue #862 : Variable ne0/neped ratio without constraint eq 81 (ne0>neped) - ! -> Potential hollowed density profile - if ( (ioptimz >= 0) .and. (.not.any(icc==81)) ) then - if ( any(ixc == 145 )) call report_error(154) - if ( any(ixc == 6 )) call report_error(155) - end if - end if - end if - ! --------------- - - - ! Cannot use Psep/R and PsepB/qAR limits at the same time - if(any(icc == 68) .and. any(icc == 56)) then - call report_error(178) - endif - - if ((any(ixc==145)) .and. (boundl(145) < fgwsep)) then !if lower bound of fgwped < fgwsep - fdiags(1) = boundl(145); fdiags(2) = fgwsep - call report_error(186) - end if - - if (any(icc == 78)) then - - !If Reinke criterion is used tesep is calculated and cannot be an - !iteration variable - if (any(ixc == 119)) then - call report_error(219) - endif - - !If Reinke criterion is used need to enforce LH-threshold - !using Martin scaling for consistency - if (.not. ilhthresh == 6) then - call report_error(218) - endif - if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then - call report_error(218) - endif - - - endif - - if (any(icc == 78)) then - - !If Reinke criterion is used tesep is calculated and cannot be an - !iteration variable - if (any(ixc == 119)) then - call report_error(219) - endif - - !If Reinke criterion is used need to enforce LH-threshold - !using Martin scaling for consistency - if (.not. ilhthresh == 6) then - call report_error(218) - endif - if (.not. any(icc==15) .and. (ipedestal .ne. 3)) then - call report_error(218) - endif - - - endif - - if (i_single_null == 0) then - idivrt = 2 - vgaptop = vgap_xpoint_divertor - shldtth = shldlth - d_vv_top = d_vv_bot - call report_error(272) - else ! i_single_null == 1 - idivrt = 1 - end if - - - ! Tight aspect ratio options (ST) - ! -------------------------------- - if ( itart == 1 ) then - - icase = 'Tight aspect ratio tokamak model' - - ! Disabled Forcing that no inboard breeding blanket is used - ! Disabled iblnkith = 0 - - ! Check if the choice of plasma current is addapted for ST - ! 2 : Peng Ip scaling (See STAR code documentation) - ! 9 : Fiesta Ip scaling - if (i_plasma_current /= 2 .and. i_plasma_current /= 9) then - idiags(1) = i_plasma_current ; call report_error(37) - end if - - !! If using Peng and Strickler (1986) model (itartpf == 0) - ! Overwrite the location of the TF coils - ! 2 : PF coil on top of TF coil - ! 3 : PF coil outside of TF coil - if (itartpf == 0) then - ipfloc(1) = 2 - ipfloc(2) = 3 - ipfloc(3) = 3 - end if - - ! Water cooled copper magnets initalisation / checks - if ( i_tf_sup == 0 ) then - ! Check if the initial centrepost coolant loop adapted to the magnet technology - ! Ice cannot flow so tcoolin > 273.15 K - if ( tcoolin < 273.15D0 ) call report_error(234) - - ! Temperature of the TF legs cannot be cooled down - if ( abs(tlegav+1.0D0) > epsilon(tlegav) .and. tlegav < 273.15D0 ) call report_error(239) - - ! Check if conductor upper limit is properly set to 50 K or below - if ( any(ixc == 20 ) .and. boundu(20) < 273.15D0 ) call report_error(241) - - ! Call a lvl 3 error if superconductor magnets are used - else if ( i_tf_sup == 1 ) then - call report_error(233) - - ! Aluminium magnets initalisation / checks - ! Initialize the CP conductor temperature to cryogenic temperature for cryo-al magnets (20 K) - else if ( i_tf_sup == 2 ) then - - ! Call a lvl 3 error if the inlet coolant temperature is too large - ! Motivation : ill-defined aluminium resistivity fit for T > 40-50 K - if ( tcoolin > 40.0D0 ) call report_error(235) - - ! Check if the leg average temperature is low enough for the resisitivity fit - if ( tlegav > 50.0D0 ) call report_error(238) - - ! Check if conductor upper limit is properly set to 50 K or below - if ( any(ixc == 20 ) .and. boundu(20) > 50.0D0 ) call report_error(240) - - ! Otherwise intitialise the average conductor temperature at - tcpav = tcoolin - - end if - - ! Check if the boostrap current selection is addapted to ST - if (i_bootstrap_current == 1) call report_error(38) - - ! Check if a single null divertor is used in double null machine - if (i_single_null == 0 .and. (ftar == 1.0 .or. ftar == 0.0)) then - call report_error(39) - end if - - ! Set the TF coil shape to picture frame (if default value) - if ( i_tf_shape == 0 ) i_tf_shape = 2 - - ! Warning stating that the CP fast neutron fluence calculation - ! is not addapted for cryoaluminium calculations yet - if ( i_tf_sup == 2 .and. any( icc == 85 ) .and. itart == 1 ) then - call report_error(260) - end if - - ! Setting the CP joints default options : - ! 0 : No joints for superconducting magents (i_tf_sup = 1) - ! 1 : Sliding joints for resistive magnets (i_tf_sup = 0, 2) - if ( i_cp_joints == -1 ) then - if ( i_tf_sup == 1 ) then - i_cp_joints = 0 - else - i_cp_joints = 1 - end if - end if - - ! Checking the CP TF top radius - if ( ( abs(r_cp_top) > epsilon(r_cp_top) .or. any(ixc(1:nvar) == 174) ) & - .and. i_r_cp_top /= 1 ) then - call report_error(267) - end if - ! -------------------------------- - - - ! Conventionnal aspect ratios specific - ! ------------------------------------ - else - - if (i_plasma_current == 2 .or. i_plasma_current == 9) call report_error(40) - - ! Set the TF coil shape to PROCESS D-shape (if default value) - if ( i_tf_shape == 0 ) i_tf_shape = 1 - - ! Check PF coil configurations - j = 0 ; k = 0 - do i = 1, ngrp - if ((ipfloc(i) /= 2).and.(ncls(i) /= 2)) then - idiags(1) = i ; idiags(2) = ncls(i) - call report_error(41) - end if - - if (ipfloc(i) == 2) then - j = j + 1 - k = k + ncls(i) - end if - end do - - if (k == 1) call report_error(42) - if (k > 2) call report_error(43) - if ((i_single_null == 1).and.(j < 2)) call report_error(44) - - ! Constraint 10 is dedicated to ST designs with demountable joints - if ( any(icc(1:neqns+nineqns) == 10 ) ) call report_error(259) - - end if - ! ------------------------------------ - - ! Pulsed power plant model - if (lpulse == 1) then - icase = 'Pulsed tokamak model' - else - esbldgm3 = 0.0D0 - end if - - ! Ensure minimum cycle time constraint is turned off - ! (not currently available, as routine thrmal has been commented out) - if ( any(icc == 42) ) then - call report_error(164) - end if - - - - ! TF coil - ! ------- - ! TF stress model not defined of r_tf_inboard = 0 - ! Unless i_tf_stress_model == 2 - ! -> If bore + gapoh + ohcth = 0 and fixed and stress constraint is used - ! Generate a lvl 3 error proposing not to use any stress constraints - if ( ( .not. ( any(ixc == 16 ) .or. any(ixc == 29 ) .or. any(ixc == 42 ) ) ) & ! No bore,gapoh, ohcth iteration - .and. ( abs(bore + gapoh + ohcth + precomp) < epsilon(bore) ) & ! bore + gapoh + ohcth = 0 - .and. ( any(icc == 31) .or. any(icc == 32) ) & ! Stress constraints (31 or 32) is used - .and. ( i_tf_stress_model /= 2 ) ) then ! TF stress model can't handle no bore - - call report_error(246) - stop 1 - end if - - ! Make sure that plane stress model is not used for resistive magnets - if ( i_tf_stress_model == 1 .and. i_tf_sup /= 1 ) call report_error(253) - - ! bucking cylinder default option setting - ! - bucking (casing) for SC i_tf_bucking ( i_tf_bucking = 1 ) - ! - No bucking for copper magnets ( i_tf_bucking = 0 ) - ! - Bucking for aluminium magnets ( i_tf_bucking = 1 ) - if ( i_tf_bucking == -1 ) then - if ( i_tf_sup == 0 ) then - i_tf_bucking = 0 - else - i_tf_bucking = 1 - end if - end if - - ! Ensure that the TF isnt placed against the - ! CS which is now outside it - if ( i_tf_bucking >= 2 .and. tf_in_cs == 1 ) then - call report_error(281) - end if - ! Ensure that no pre-compression structure - ! is used for bucked and wedged design - if ( i_tf_bucking >= 2 .and. iprecomp == 1 ) then - call report_error(252) - end if - - ! Number of stress calculation layers - ! +1 to add in the inboard TF coil case on the plasma side, per Issue #1509 - n_tf_stress_layers = i_tf_bucking + n_tf_graded_layers + 1 - - ! If TFC sidewall has not been set by user - if ( casths < 0.1d-10 ) tfc_sidewall_is_fraction = .true. - - ! If inboard TF coil case plasma side thickness has not been set by user - if( casthi < 0.1d-10 ) casthi_is_fraction = .true. - - ! Setting the default cryo-plants efficiencies - !-! - if ( abs(eff_tf_cryo + 1.0D0) < epsilon(eff_tf_cryo) ) then - - ! The ITER cyoplant efficiency is used for SC - if ( i_tf_sup == 1 ) then - eff_tf_cryo = 0.13D0 - - ! Strawbrige plot extrapolation is used for Cryo-Al - else if ( i_tf_sup == 2 ) then - eff_tf_cryo = 0.40D0 - end if - - ! Cryo-plane efficiency must be in [0-1.0] - else if ( eff_tf_cryo > 1.0D0 .or. eff_tf_cryo < 0.0D0 ) then - call report_error(248) - stop 1 - end if - !-! - - ! Integer turns option not yet available for REBCO taped turns - !-! - if ( i_tf_sc_mat == 6 .and. i_tf_turns_integer == 1 ) then - call report_error(254) - stop 1 - end if - !-! - - - ! Setting up insulation layer young modulae default values [Pa] - !-! - if ( abs(eyoung_ins - 1.0D8 ) < epsilon(eyoung_ins) ) then - - ! Copper magnets, no insulation material defined - ! But use the ITER design by default - if ( i_tf_sup == 0 ) then - eyoung_ins = 20.0D9 - - ! SC magnets - ! Value from DDD11-2 v2 2 (2009) - else if ( i_tf_sup == 1 ) then - eyoung_ins = 20.0D9 - - ! Cryo-aluminum magnets (Kapton polymer) - else if ( i_tf_sup == 2 ) then - eyoung_ins = 2.5D9 - end if - end if - !-! - - !-! Setting the default WP geometry - !-! - if ( i_tf_wp_geom == -1 ) then - if ( i_tf_turns_integer == 0 ) i_tf_wp_geom = 1 - if ( i_tf_turns_integer == 1 ) i_tf_wp_geom = 0 - end if - !-! - - !-! Setting the TF coil conductor elastic properties - !-! - if ( i_tf_cond_eyoung_axial == 0 ) then - ! Conductor stiffness is not considered - eyoung_cond_axial = 0 - eyoung_cond_trans = 0 - else if ( i_tf_cond_eyoung_axial == 2 ) then - ! Select sensible defaults from the literature - select case (i_tf_sc_mat) - case (1,4,5) - ! Nb3Sn: Nyilas, A et. al, Superconductor Science and Technology 16, no. 9 (2003): 1036–42. https://doi.org/10.1088/0953-2048/16/9/313. - eyoung_cond_axial = 32D9 - case (2) - ! Bi-2212: Brown, M. et al, IOP Conference Series: Materials Science and Engineering 279 (2017): 012022. https://doi.org/10.1088/1757-899X/279/1/012022. - eyoung_cond_axial = 80D9 - case (3,7) - ! NbTi: Vedrine, P. et. al, IEEE Transactions on Applied Superconductivity 9, no. 2 (1999): 236–39. https://doi.org/10.1109/77.783280. - eyoung_cond_axial = 6.8D9 - case (6,8,9) - ! REBCO: Fujishiro, H. et. al, Physica C: Superconductivity, 426–431 (2005): 699–704. https://doi.org/10.1016/j.physc.2005.01.045. - eyoung_cond_axial = 145D9 - end select - - if ( i_tf_cond_eyoung_trans == 0) then - ! Transverse stiffness is not considered - eyoung_cond_trans = 0 - else - ! Transverse stiffness is significant - eyoung_cond_trans = eyoung_cond_axial - end if - end if - !-! - - ! Check if the WP/conductor radial thickness (dr_tf_wp) is large enough - ! To contains the insulation, cooling and the support structure - ! Rem : Only verified if the WP thickness is used - if ( any(ixc(1:nvar) == 140) ) then - - ! Minimal WP thickness - if ( i_tf_sup == 1 ) then - dr_tf_wp_min = 2.0D0 * ( tinstf + tfinsgap + thicndut + dhecoil ) - - ! Steel conduit thickness (can be an iteration variable) - if ( any(ixc(1:nvar) == 58 ) ) then - dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * boundl(58) - else - dr_tf_wp_min = dr_tf_wp_min + 2.0D0 * thwcndut - end if - - ! Minimal conductor layer thickness - else if ( i_tf_sup == 0 .or. i_tf_sup == 2 ) then - dr_tf_wp_min = 2.0D0 * ( thicndut + tinstf ) + 4.0D0 * rcool - end if - - if ( boundl(140) < dr_tf_wp_min ) then - fdiags(1) = dr_tf_wp_min - call report_error(255) - end if - end if - - ! Setting t_turn_tf_is_input to true if t_turn_tf is an input - if ( abs(t_turn_tf) < epsilon(t_turn_tf) ) then - t_turn_tf_is_input = .false. - else - t_turn_tf_is_input = .true. - end if - - ! Impossible to set the turn size of integer turn option - if ( t_turn_tf_is_input .and. i_tf_turns_integer == 1 ) then - call report_error(269) - end if - - if ( i_tf_wp_geom /= 0 .and. i_tf_turns_integer == 1 ) then - call report_error(283) - end if - - if ( i_bootstrap_current == 5 .and. i_diamagnetic_current /= 0 ) then - call report_error(284) - end if - - ! Setting t_cable_tf_is_input to true if t_cable_tf is an input - if ( abs(t_cable_tf) < epsilon(t_cable_tf) ) then - t_cable_tf_is_input = .false. - else - t_cable_tf_is_input = .true. - end if - - ! Impossible to set the cable size of integer turn option - if ( t_cable_tf_is_input .and. i_tf_turns_integer == 1 ) then - call report_error(269) - end if - - ! Impossible to set both the TF coil turn and the cable dimension - if ( t_turn_tf_is_input .and. t_cable_tf_is_input ) then - call report_error(271) - end if - - ! Checking the SC temperature for LTS - if ( ( i_tf_sc_mat == 1 .or. & - i_tf_sc_mat == 3 .or. & - i_tf_sc_mat == 4 .or. & - i_tf_sc_mat == 5 ) .and. tftmp > 10.0D0 ) then - call report_error(270) - end if - ! ------- - - - - ! PF coil resistivity is zero if superconducting - if (ipfres == 0) pfclres = 0.0D0 - - ! If there is no NBI, then hot beam density should be zero - if (irfcd == 1) then - if ((iefrf /= 5).and.(iefrf /= 8)) rnbeam = 0.0D0 - else - rnbeam = 0.0D0 - end if - - ! Set inboard blanket thickness to zero if no inboard blanket switch - ! used (Issue #732) - if (iblnkith == 0) blnkith = 0.0D0 - - ! Solid breeder assumed if ipowerflow=0 - - !if (ipowerflow == 0) blkttype = 3 - - ! Set coolant fluid type - - !if ((blkttype == 1).or.(blkttype == 2)) then - ! coolwh = 2 ! water - !else - ! coolwh = 1 ! helium - !end if - - ! But... set coolant to water if blktmodel > 0 - ! Although the *blanket* is by definition helium-cooled in this case, - ! the shield etc. are assumed to be water-cooled, and since water is - ! heavier (and the unit cost of pumping it is higher), the calculation - ! for coolmass is better done with coolwh=2 if blktmodel > 0 to give - ! slightly pessimistic results. - - !if (blktmodel > 0) then - ! secondary_cycle = 0 - ! blkttype = 3 ! HCPB - ! coolwh = 2 - !end if - - ! Ensure that blanket material fractions allow non-zero space for steel - ! CCFE HCPB Model - - if (istell == 0) then - if ((iblanket == 1).or.(iblanket == 3)) then - fsum = breeder_multiplier + vfcblkt + vfpblkt - if (fsum >= 1.0D0) then - idiags(1) = iblanket - fdiags(2) = breeder_multiplier - fdiags(3) = vfcblkt - fdiags(4) = vfpblkt - fdiags(5) = fsum - call report_error(165) - end if - end if - end if - - ! Initialise superconductor cable parameters - if(i_tf_sup==1)then - call initialise_cables() - end if - - ! Check that the temperature margins are not overdetermined - if(tmargmin>0.0001d0)then - ! This limit has been input and will be applied to both TFC and CS - if(tmargmin_tf>0.0001d0)then - write(*,*)'tmargmin_tf and tmargmin should not both be specified in IN.DAT.' - write(*,*)'tmargmin_tf has been ignored.' - end if - if(tmargmin_cs>0.0001d0)then - write(*,*)'tmargmin_cs and tmargmin should not both be specified in IN.DAT.' - write(*,*)'tmargmin_cs has been ignored.' - end if - tmargmin_tf = tmargmin - tmargmin_cs = tmargmin - end if - - if (tauee_in.ge.1.0D-10.and.isc.ne.48) then - ! Report error if confinement time is in the input - ! but the scaling to use it is not selected. - call report_error(220) - end if - - if (aspect.gt.1.7D0.and.isc.eq.46) then - ! NSTX scaling is for A<1.7 - call report_error(221) - end if - - if (i_plasma_current.eq.2.and.isc.eq.42) then - call report_error(222) - end if - - ! Cannot use temperature margin constraint with REBCO TF coils - if(any(icc == 36) .and. ((i_tf_sc_mat == 8).or.(i_tf_sc_mat == 9))) then - call report_error(265) - endif - - ! Cannot use temperature margin constraint with REBCO CS coils - if(any(icc == 60) .and. (isumatoh == 8)) then - call report_error(264) - endif - - ! Cold end of the cryocooler should be colder than the TF - if(tmpcry > tftmp) then - call report_error(273) - endif - - ! Cannot use TF coil strain limit if i_str_wp is off: - if(any(icc == 88) .and. (i_str_wp == 0)) then - call report_error(275) - endif - - errors_on = .false. - - ! Disable error logging only after all checks have been performed. - ! (CPSS #1582: Why is error logging disabled at all?) - errors_on = .false. - - - end subroutine check - end module init_module