! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! module scan_module !! Module containing routines to perform a parameter scan !! author: P J Knight, CCFE, Culham Science Centre !! None !! This module contains routines to perform a parameter scan !! over a range of values of a particular scanning variable. !! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! #ifndef dp use, intrinsic :: iso_fortran_env, only: dp=>real64 #endif implicit none public integer, parameter :: ipnscns = 1000 !! Maximum number of scan points integer, parameter :: ipnscnv = 81 !! Number of available scan variables integer, parameter :: noutvars = 84 integer, parameter :: width = 110 integer :: scan_dim !! 1-D or 2-D scan switch (1=1D, 2=2D) integer :: isweep !! Number of scan points to calculate integer :: isweep_2 !! Number of 2D scan points to calculate integer :: nsweep !! Switch denoting quantity to scan:<UL> !! <LI> 1 aspect !! <LI> 2 hldivlim !! <LI> 3 pnetelin !! <LI> 4 hfact !! <LI> 5 oacdcp !! <LI> 6 walalw !! <LI> 7 beamfus0 !! <LI> 8 fqval !! <LI> 9 te !! <LI> 10 boundu(15: fvs) !! <LI> 11 dnbeta !! <LI> 12 bootstrap_current_fraction_max !! <LI> 13 boundu(10: hfact) !! <LI> 14 fiooic !! <LI> 15 fjprot !! <LI> 16 rmajor !! <LI> 17 bmxlim !! <LI> 18 gammax !! <LI> 19 boundl(16: ohcth) !! <LI> 20 t_burn_min !! <LI> 21 not used !! <LI> 22 cfactr (N.B. requires iavail=0) !! <LI> 23 boundu(72: fipir) !! <LI> 24 powfmax !! <LI> 25 kappa !! <LI> 26 triang !! <LI> 27 tbrmin (for blktmodel > 0 only) !! <LI> 28 bt !! <LI> 29 coreradius !! <LI> 30 fimpvar # OBSOLETE !! <LI> 31 taulimit !! <LI> 32 epsvmc !! <LI> 33 ttarget !! <LI> 34 qtargettotal !! <LI> 35 lambda_q_omp !! <LI> 36 lambda_target !! <LI> 37 lcon_factor !! <LI> 38 Neon upper limit !! <LI> 39 Argon upper limit !! <LI> 40 Xenon upper limit !! <LI> 41 blnkoth !! <LI> 42 Argon fraction fimp(9) !! <LI> 43 normalised minor radius at which electron cyclotron current drive is maximum !! <LI> 44 Allowable maximum shear stress (Tresca) in tf coil structural material !! <LI> 45 Minimum allowable temperature margin ; tf coils !! <LI> 46 boundu(150) fgwsep !! <LI> 47 impurity_enrichment(9) Argon impurity enrichment !! <LI> 48 TF coil - n_pancake (integer turn winding pack) !! <LI> 49 TF coil - n_layer (integer turn winding pack) !! <LI> 50 Xenon fraction fimp(13) !! <LI> 51 Power fraction to lower DN Divertor ftar !! <LI> 52 SoL radiation fraction !! <LI> 54 GL_nbti upper critical field at 0 Kelvin !! <LI> 55 `shldith` : Inboard neutron shield thickness !! <LI> 56 crypmw_max: Maximum cryogenic power (ixx=164, ixc=87) !! <LI> 57 `bt` lower boundary !! <LI> 58 `scrapli` : Inboard plasma-first wall gap !! <LI> 59 `scraplo` : Outboard plasma-first wall gap !! <LI> 60 sig_tf_wp_max: Allowable stress in TF Coil conduit (Tresca) !! <LI> 61 copperaoh_m2_max : CS coil current / copper area !! <LI> 62 coheof : CS coil current density at EOF !! <LI> 63 ohcth : CS thickness (m) !! <LI> 64 ohhghf : CS height (m) !! <LI> 65 n_cycle_min : Minimum cycles for CS stress model constraint 90 !! <LI> 66 oh_steel_frac: Steel fraction in CS coil !! <LI> 67 t_crack_vertical: Initial crack vertical dimension (m) </UL> !! <LI> 68 `inlet_temp_liq' : Inlet temperature of blanket liquid metal coolant/breeder (K) !! <LI> 69 `outlet_temp_liq' : Outlet temperature of blanket liquid metal coolant/breeder (K) !! <LI> 70 `blpressure_liq' : Blanket liquid metal breeder/coolant pressure (Pa) !! <LI> 71 `n_liq_recirc' : Selected number of liquid metal breeder recirculations per day !! <LI> 72 `bz_channel_conduct_liq' : Conductance of liquid metal breeder duct walls (A V-1 m-1) !! <LI> 73 `pnuc_fw_ratio_dcll' : Ratio of FW nuclear power as fraction of total (FW+BB) !! <LI> 74 `f_nuc_pow_bz_struct' : Fraction of BZ power cooled by primary coolant for dual-coolant balnket !! <LI> 75 pitch : pitch of first wall cooling channels (m) !! <LI> 76 etath : Thermal conversion eff. !! <LI> 77 startupratio : Gyrotron redundancy !! <LI> 78 fkind : Multiplier for Nth of a kind costs !! <LI> 79 etaech : ECH wall plug to injector efficiency integer :: nsweep_2 !! nsweep_2 /3/ : switch denoting quantity to scan for 2D scan: real(dp), dimension(ipnscns) :: sweep !! sweep(ipnscns) /../: actual values to use in scan real(dp), dimension(ipnscns) :: sweep_2 !! sweep_2(ipnscns) /../: actual values to use in 2D scan ! Vars in subroutines scan_1d and scan_2d requiring re-initialising before ! each new run logical :: first_call_1d logical :: first_call_2d contains subroutine init_scan_module !! Initialise module variables implicit none scan_dim = 1 isweep = 0 isweep_2 = 0 nsweep = 1 nsweep_2 = 3 sweep = 0.0D0 sweep_2 = 0.0D0 first_call_1d = .true. first_call_2d = .true. end subroutine init_scan_module subroutine scan_1d_write_point_header(iscan) use global_variables, only: iscan_global, xlabel, vlabel use constants, only: mfile, nout use process_output, only: ovarin, ostars, oblnkl implicit none integer, intent(in) :: iscan !! Scan point number ! Makes iscan available globally (read-only) iscan_global = iscan call scan_select(nsweep, sweep, iscan, vlabel, xlabel) ! Write banner to output file call oblnkl(nout) call ostars(nout,width) write(nout,10) ' ***** Scan point ', iscan,' of ',isweep,': & ',trim(xlabel),', ',trim(vlabel),' = ',sweep(iscan),' *****' 10 format(a,i2,a,i2,5a,1pe10.3,a) call ostars(nout,width) ! Write additional information to mfile call oblnkl(mfile) call ovarin(mfile,'Scan point number','(iscan)',iscan) ! Call the optimization routine VMCON at this scan point write(*,20)'Starting scan point ',iscan,' of ',isweep,': ', trim(xlabel),', & ',trim(vlabel),' = ',sweep(iscan) 20 format(a,i2,a,i2,a,4a,1pe10.3) end subroutine scan_1d_write_point_header subroutine scan_1d_store_output(iscan, ifail, noutvars_, ipnscns_, outvar) use constraint_variables, only: taulimit use cost_variables, only: cdirt, coe, coeoam, coefuelt, c222, ireactor, & capcost, coecap, c221 use current_drive_variables, only: pheat, pinjmw, bootstrap_current_fraction, beam_energy, bigq use divertor_variables, only: hldiv use error_handling, only: errors_on use heat_transport_variables, only: pgrossmw, pinjwp, pnetelmw use impurity_radiation_module, only: fimp use pfcoil_variables, only: whtpf use pf_power_variables, only: srcktpm use process_output, only: oblnkl use numerics, only: sqsumsq use tfcoil_variables, only: tfareain, wwp2, sig_tf_wp, tfcmw, tcpmax, oacdcp, & tfcpmw, fcutfsu, acond, fcoolcp, rcool, whttf, ppump, vcool, wwp1, n_tf, & dr_tf_wp, b_crit_upper_nbti use fwbs_variables, only: tpeak use physics_variables, only: q, aspect, pradmw, dene, fusion_power, btot, tesep, & pdivt, ralpne, ten, betap, hfac, teped, alpha_power_beams, qlim, rmajor, wallmw, & beta, betalim, bt, plasma_current use global_variables, only: verbose, maxcal, runtitle, run_tests use constants, only: nout implicit none integer, intent(in) :: iscan integer, intent(in) :: ifail ! outvar integer, intent(in) :: noutvars_, ipnscns_ real(dp), dimension(noutvars_,ipnscns_), intent(out) :: outvar ! Turn off error reporting (until next output) errors_on = .false. ! Store values for MFILE.DAT output outvar( 1,iscan) = dble(ifail) outvar( 2,iscan) = sqsumsq outvar( 3,iscan) = coe outvar( 4,iscan) = coecap outvar( 5,iscan) = coefuelt outvar( 6,iscan) = coeoam outvar( 7,iscan) = capcost outvar( 8,iscan) = c221 + c222 outvar( 9,iscan) = cdirt / 1.0D3 outvar(10,iscan) = rmajor outvar(11,iscan) = aspect outvar(12,iscan) = 1.0D-6 * plasma_current outvar(13,iscan) = bt outvar(14,iscan) = btot outvar(15,iscan) = q outvar(16,iscan) = qlim outvar(17,iscan) = beta outvar(18,iscan) = betalim outvar(19,iscan) = betap / aspect outvar(20,iscan) = ten/10.0D0 outvar(21,iscan) = dene/1.0D20 outvar(22,iscan) = hfac(6) outvar(23,iscan) = hfac(7) outvar(24,iscan) = fusion_power outvar(25,iscan) = alpha_power_beams * 5.0D0 outvar(26,iscan) = wallmw outvar(27,iscan) = pinjmw outvar(28,iscan) = pinjwp outvar(29,iscan) = pheat outvar(30,iscan) = pinjmw - pheat outvar(31,iscan) = bigq outvar(32,iscan) = bootstrap_current_fraction outvar(33,iscan) = beam_energy/1.0D3 outvar(34,iscan) = hldiv outvar(35,iscan) = tfcmw outvar(36,iscan) = whttf outvar(37,iscan) = sig_tf_wp outvar(38,iscan) = oacdcp/1.0D6 outvar(39,iscan) = tcpmax outvar(40,iscan) = tfcpmw outvar(41,iscan) = fcoolcp outvar(42,iscan) = rcool outvar(43,iscan) = vcool outvar(44,iscan) = ppump/1.0D6 outvar(45,iscan) = 1.0D-3 * srcktpm outvar(46,iscan) = whtpf outvar(47,iscan) = pgrossmw outvar(48,iscan) = pnetelmw if (ireactor == 1) then outvar(49,iscan) = (pgrossmw-pnetelmw) / pgrossmw else outvar(49,iscan) = 0.0D0 end if outvar(50,iscan) = pdivt/rmajor !outvar(51,iscan) = fimpvar #OBSOLETE outvar(51,iscan) = 0.0d0 outvar(52,iscan) = pradmw outvar(53,iscan) = tpeak outvar(54,iscan) = fcutfsu outvar(55,iscan) = (wwp1+wwp2)*dr_tf_wp outvar(56,iscan) = acond outvar(57,iscan) = tfareain/n_tf outvar(58,iscan) = taulimit outvar(66,iscan) = ralpne outvar(69,iscan) = fimp(1) outvar(70,iscan) = fimp(2) outvar(71,iscan) = fimp(3) outvar(72,iscan) = fimp(4) outvar(73,iscan) = fimp(5) outvar(74,iscan) = fimp(6) outvar(75,iscan) = fimp(7) outvar(76,iscan) = fimp(8) outvar(77,iscan) = fimp(9) outvar(78,iscan) = fimp(10) outvar(79,iscan) = fimp(11) outvar(80,iscan) = fimp(12) outvar(81,iscan) = fimp(13) outvar(82,iscan) = fimp(14) outvar(83,iscan) = teped end subroutine scan_1d_store_output subroutine scan_1d_write_plot(iscan, outvar) use global_variables, only: icase, xlabel use constants, only: nplot, mfile use process_output, only: ovarin implicit none integer, intent(inout) :: iscan real(dp), dimension(:,:), intent(in) :: outvar character(len=48) :: tlabel integer :: ivar character(len=25), dimension(noutvars), save :: plabel tlabel = icase ! Set up labels for plotting output ! Use underscores instead of spaces if (first_call_1d) then plabel( 1) = 'Ifail____________________' plabel( 2) = 'Sqsumsq__________________' plabel( 3) = 'Electric_cost_(mil/kwh)__' plabel( 4) = 'Capital_cost_(mil/kwh)___' plabel( 5) = 'Fuel_cost_(mil/kwh)______' plabel( 6) = 'Operations_cost_(mil/kwh)' plabel( 7) = 'Capital_cost_(millions)__' plabel( 8) = 'Core_costs_(millions)____' plabel( 9) = 'Direct_cost_(billions)___' plabel(10) = 'Major_Radius_(m)_________' plabel(11) = 'Aspect_Ratio_____________' plabel(12) = 'Plasma_Current_(MA)______' plabel(13) = 'B_Toroidal_Axis_(T)______' plabel(14) = 'B_total_on_axis_(T)______' plabel(15) = 'Safety_Factor____________' plabel(16) = 'qlim_(zero_if_ishape=0)__' plabel(17) = 'Beta_____________________' plabel(18) = 'Beta_Limit_______________' plabel(19) = 'Epsilon_Beta_Poloidal____' plabel(20) = 'Dens.weight_Te_(10keV)___' plabel(21) = 'Average_Dens_(10^20/m^3)_' plabel(22) = 'H-fact_Iter_Power________' plabel(23) = 'H-fact_Iter_Offset_______' plabel(24) = 'Fusion_Power_(MW)________' plabel(25) = 'nb_Fusion_Power_(MW)_____' plabel(26) = 'Wall_Load_(MW/m^2)_______' plabel(27) = 'Injection_Power_(MW)_____' plabel(28) = 'Inject_Pwr_Wall_Plug_(MW)' plabel(29) = 'Heating_Power_(MW)_______' plabel(30) = 'Current_Drive_(MW)_______' plabel(31) = 'Big_Q____________________' plabel(32) = 'Bootstrap_Fraction_______' plabel(33) = 'Neutral_Beam_Energy_(MeV)' plabel(34) = 'Divertor_Heat_(MW/m^2)___' plabel(35) = 'TF_coil_Power_(MW)_______' plabel(36) = 'TF_coil_weight_(kg)______' plabel(37) = 'vM_stress_in_TF_case_(Pa)' plabel(38) = 'J_TF_inboard_leg_(MA/m^2)' plabel(39) = 'Centrepost_max_T_(TART)__' plabel(40) = 'Res_TF_inbrd_leg_Pwr_(MW)' plabel(41) = 'Coolant_Fraction_Ctr.____' plabel(42) = 'C/P_coolant_radius_(m)___' plabel(43) = 'C/P_coolant_velocity(m/s)' plabel(44) = 'C/P_pump_power_(MW)______' plabel(45) = 'PF_coil_Power_(MW)_______' plabel(46) = 'PF_coil_weight_(kg)______' plabel(47) = 'Gross_Elect_Pwr_(MW)_____' plabel(48) = 'Net_electric_Pwr_(MW)____' plabel(49) = 'Recirculating_Fraction___' plabel(50) = 'Psep/R___________________' plabel(51) = '' !OBSOLETE plabel(52) = 'Tot._radiation_power_(MW)' plabel(53) = 'First_wall_peak_temp_(K)_' plabel(54) = 'Cu_frac_TFC_conductor____' plabel(55) = 'Winding_pack_area_TFC(m2)' plabel(56) = 'Conductor_area_TFC_(m2)__' plabel(57) = 'Area_TF_inboard_leg_(m2)_' plabel(58) = 'Taup/taueff_lower_limit__' plabel(59) = 'Plasma_temp_at_sep_[keV]_' plabel(60) = 'SOL_density_at_OMP_______' plabel(61) = 'Power_through__separatrix' plabel(62) = 'neomp/nesep______________' plabel(63) = 'qtargettotal_____________' plabel(64) = 'Total_pressure_at_target_' plabel(65) = 'Temperature_at_target____' plabel(66) = 'Helium_fraction__________' plabel(67) = 'Momentum_loss_factor_____' plabel(68) = 'totalpowerlost_[W]_______' plabel(69) = 'H__concentration_________' plabel(70) = 'He_concentration_________' plabel(71) = 'Be_concentration_________' plabel(72) = 'C__concentration_________' plabel(73) = 'N__concentration_________' plabel(74) = 'O__concentration_________' plabel(75) = 'Ne_concentration_________' plabel(76) = 'Si_concentration_________' plabel(77) = 'Ar_concentration_________' plabel(78) = 'Fe_concentration_________' plabel(79) = 'Ni_concentration_________' plabel(80) = 'Kr_concentration_________' plabel(81) = 'Xe_concentration_________' plabel(82) = 'W__concentration_________' plabel(83) = 'teped____________________' plabel(84) = 'Max_field_on_TF_coil_____' call ovarin(mfile,'Number of scan points','(isweep)',isweep) call ovarin(mfile,'Scanning variable number','(nsweep)',nsweep) first_call_1d = .false. end if end subroutine scan_1d_write_plot subroutine scan_2d_init !! Routine to call 2-D scan !! author: J Morris, UKAEA, Culham Science Centre ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use constants, only: mfile use process_output, only: ovarin implicit none ! Set up labels for plotting output ! Use underscores instead of spaces call ovarin(mfile,'Number of first variable scan points','(isweep)',isweep) call ovarin(mfile,'Number of second variable scan points','(isweep_2)',isweep_2) call ovarin(mfile,'Scanning first variable number','(nsweep)',nsweep) call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) call ovarin(mfile,'Scanning second variable number','(nsweep_2)',nsweep_2) end subroutine scan_2d_init subroutine scan_2d_write_point_header(iscan, iscan_1, iscan_2, iscan_R) use process_output, only: oblnkl, ostars, ovarin use global_variables, only: vlabel, vlabel_2, xlabel, xlabel_2, iscan_global use constants, only: nout, mfile implicit none integer, intent(in) :: iscan integer, intent(in) :: iscan_1 integer, intent(in) :: iscan_2 integer, intent(out) :: iscan_R integer :: ifail if (mod(iscan_1,2)==0) then iscan_R = isweep_2 - iscan_2 + 1 else iscan_R = iscan_2 end if ! Makes iscan available globally (read-only) iscan_global = iscan call scan_select(nsweep, sweep, iscan_1, vlabel, xlabel) call scan_select(nsweep_2, sweep_2, iscan_R, vlabel_2, xlabel_2) ! Write banner to output file call oblnkl(nout) call ostars(nout,width) write(nout,10) iscan, isweep*isweep_2, trim(vlabel), & sweep(iscan_1), trim(vlabel_2), sweep_2(iscan_R) ! 10 format(a, i2, a, i2, 5a, 1pe10.3, a) 10 format(' ***** 2D scan point ', i3, ' of ', i3, ' : ', a, ' = ', & 1pe10.3, ' and ', a, ' = ', 1pe10.3, ' *****') call ostars(nout,width) ! Write additional information to mfile call oblnkl(mfile) call ovarin(mfile,'Scan point number','(iscan)',iscan) ! Call the optimization routine VMCON at this scan point write(*,20) iscan, trim(xlabel), trim(vlabel), sweep(iscan_1), & trim(xlabel_2), trim(vlabel_2), sweep_2(iscan_R) ! 20 format(a,i2,a,4a,1pe10.3) 20 format('Starting scan point ', i3, ': ', a, ', ', a, ' = ', & 1pe10.3, ' and ', a, ', ', a, ' = ', 1pe10.3) end subroutine scan_2d_write_point_header subroutine scan_2d_store_output(ifail, iscan_1, iscan_R, iscan, noutvars_, ipnscns_, outvar, & sweep_1_vals, sweep_2_vals) implicit none integer, intent(in) :: ifail integer, intent(in) :: iscan_1 integer, intent(in) :: iscan_R integer, intent(in) :: iscan integer, intent(in) :: noutvars_, ipnscns_ ! Required for shape of intent(out) arrays real(dp), dimension(noutvars_,ipnscns_), intent(out) :: outvar real(dp), dimension(ipnscns_), intent(out) :: sweep_1_vals, sweep_2_vals call scan_1d_store_output(iscan, ifail, noutvars_, ipnscns_, outvar) sweep_1_vals(iscan) = sweep(iscan_1) sweep_2_vals(iscan) = sweep_2(iscan_R) end subroutine scan_2d_store_output subroutine scan_2d_write_plot(iscan, outvar, sweep_1_vals, sweep_2_vals) use constants, only: nplot use global_variables, only: icase, xlabel, xlabel_2 implicit none integer, intent(inout) :: iscan real(dp), dimension(:,:), intent(in) :: outvar real(dp), dimension(:), intent(in) :: sweep_1_vals, sweep_2_vals integer :: ivar character(len=48) :: tlabel character(len=25), dimension(noutvars), save :: plabel plabel( 1) = 'Ifail____________________' plabel( 2) = 'Sqsumsq__________________' plabel( 3) = 'Electric_cost_(mil/kwh)__' plabel( 4) = 'Capital_cost_(mil/kwh)___' plabel( 5) = 'Fuel_cost_(mil/kwh)______' plabel( 6) = 'Operations_cost_(mil/kwh)' plabel( 7) = 'Capital_cost_(millions)__' plabel( 8) = 'Core_costs_(millions)____' plabel( 9) = 'Direct_cost_(billions)___' plabel(10) = 'Major_Radius_(m)_________' plabel(11) = 'Aspect_Ratio_____________' plabel(12) = 'Plasma_Current_(MA)______' plabel(13) = 'B_Toroidal_Axis_(T)______' plabel(14) = 'B_total_on_axis_(T)______' plabel(15) = 'Safety_Factor____________' plabel(16) = 'qlim_(zero_if_ishape=0)__' plabel(17) = 'Beta_____________________' plabel(18) = 'Beta_Limit_______________' plabel(19) = 'Epsilon_Beta_Poloidal____' plabel(20) = 'Dens.weight_Te_(10keV)___' plabel(21) = 'Average_Dens_(10^20/m^3)_' plabel(22) = 'H-fact_Iter_Power________' plabel(23) = 'H-fact_Iter_Offset_______' plabel(24) = 'Fusion_Power_(MW)________' plabel(25) = 'nb_Fusion_Power_(MW)_____' plabel(26) = 'Wall_Load_(MW/m^2)_______' plabel(27) = 'Injection_Power_(MW)_____' plabel(28) = 'Inject_Pwr_Wall_Plug_(MW)' plabel(29) = 'Heating_Power_(MW)_______' plabel(30) = 'Current_Drive_(MW)_______' plabel(31) = 'Big_Q____________________' plabel(32) = 'Bootstrap_Fraction_______' plabel(33) = 'Neutral_Beam_Energy_(MeV)' plabel(34) = 'Divertor_Heat_(MW/m^2)___' plabel(35) = 'TF_coil_Power_(MW)_______' plabel(36) = 'TF_coil_weight_(kg)______' plabel(37) = 'vM_stress_in_TF_cond_(Pa)' plabel(38) = 'J_TF_inboard_leg_(MA/m^2)' plabel(39) = 'Centrepost_max_T_(TART)__' plabel(40) = 'Res_TF_inbrd_leg_Pwr_(MW)' plabel(41) = 'Coolant_Fraction_Ctr.____' plabel(42) = 'C/P_coolant_radius_(m)___' plabel(43) = 'C/P_coolant_velocity(m/s)' plabel(44) = 'C/P_pump_power_(MW)______' plabel(45) = 'PF_coil_Power_(MW)_______' plabel(46) = 'PF_coil_weight_(kg)______' plabel(47) = 'Gross_Elect_Pwr_(MW)_____' plabel(48) = 'Net_electric_Pwr_(MW)____' plabel(49) = 'Recirculating_Fraction___' plabel(50) = 'Psep/R___________________' plabel(51) = '' !OBSOLETE plabel(52) = 'Tot._radiation_power_(MW)' plabel(53) = 'First_wall_peak_temp_(K)_' plabel(54) = 'Cu_frac_TFC_conductor____' plabel(55) = 'Winding_pack_area_TFC(m2)' plabel(56) = 'Conductor_area_TFC_(m2)__' plabel(57) = 'Area_TF_inboard_leg_(m2)_' plabel(58) = 'Taup/taueff_lower_limit__' plabel(59) = 'Plasma_temp_at_sep_[keV]_' plabel(60) = 'SOL_density_at_OMP_______' plabel(61) = 'Power_through__separatrix' plabel(62) = 'neomp/nesep______________' plabel(63) = 'qtargettotal_____________' plabel(64) = 'Total_pressure_at_target_' plabel(65) = 'Temperature_at_target____' plabel(66) = 'Helium_fraction__________' plabel(67) = 'Momentum_loss_factor_____' plabel(68) = 'totalpowerlost_[W]_______' plabel(69) = 'H__concentration_________' plabel(70) = 'He_concentration_________' plabel(71) = 'Be_concentration_________' plabel(72) = 'C__concentration_________' plabel(73) = 'N__concentration_________' plabel(74) = 'O__concentration_________' plabel(75) = 'Ne_concentration_________' plabel(76) = 'Si_concentration_________' plabel(77) = 'Ar_concentration_________' plabel(78) = 'Fe_concentration_________' plabel(79) = 'Ni_concentration_________' plabel(80) = 'Kr_concentration_________' plabel(81) = 'Xe_concentration_________' plabel(82) = 'W__concentration_________' plabel(83) = 'teped____________________' plabel(84) = 'Max_field_on_TF_coil_____' tlabel = icase end subroutine scan_2d_write_plot subroutine scan_select(nwp, swp, iscn, vlab, xlab) !! Routine to select first scan case !! author: J Morris, UKAEA, Culham Science Centre ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use build_variables, only: blnkoth, shldith, scrapli, scraplo, ohcth use constraint_variables, only: fiooic, walalw, bmxlim, fqval, taulimit, & gammax, t_burn_min, tbrmin, fjprot, pnetelin, powfmax use cost_variables, only: cfactr, iavail, fkind, startupratio use current_drive_variables, only: bootstrap_current_fraction_max, etaech use divertor_variables, only: hldivlim use error_handling, only: idiags, report_error use fwbs_variables, only: inlet_temp_liq, outlet_temp_liq, blpressure_liq, & n_liq_recirc, bz_channel_conduct_liq, pnuc_fw_ratio_dcll, f_nuc_pow_bz_struct, pitch use impurity_radiation_module, only: fimp, coreradius, impurity_arr_frac use physics_variables, only: kappa, dnbeta, te, aspect, ftar, bt, & rad_fraction_sol, triang, rmajor, beamfus0, hfact use numerics, only: epsvmc, boundu, boundl use tfcoil_variables, only: tmargmin_tf, sig_tf_case_max, n_pancake, oacdcp, & n_layer, b_crit_upper_nbti, sig_tf_wp_max, fcoolcp, n_tf_turn use heat_transport_variables, only: crypmw_max, etath use rebco_variables, only: copperaoh_m2_max use pfcoil_variables, only: coheof, ohhghf, oh_steel_frac use CS_fatigue_variables, only: n_cycle_min, t_crack_vertical implicit none ! Arguments integer, intent(in) :: nwp, iscn real(dp), intent(in), dimension(:) :: swp character(len=25), intent(out) :: vlab, xlab select case (nwp) ! Use underscores instead of spaces in xlabel ! MDK Remove the "=" from vlabel, to make it easier to compare with ! list of iteration variables case (1) aspect = swp(iscn) vlab = 'aspect' ; xlab = 'Aspect_ratio' case (2) hldivlim = swp(iscn) vlab = 'hldivlim' ; xlab = 'Div_heat_limit_(MW/m2)' case (3) pnetelin = swp(iscn) vlab = 'pnetelin' ; xlab = 'Net_electric_power_(MW)' case (4) hfact = swp(iscn) vlab = 'hfact' ; xlab = 'Confinement_H_factor' case (5) oacdcp = swp(iscn) vlab = 'oacdcp' ; xlab = 'TF_inboard_leg_J_(MA/m2)' case (6) walalw = swp(iscn) vlab = 'walalw' ; xlab = 'Allow._wall_load_(MW/m2)' case (7) beamfus0 = swp(iscn) vlab = 'beamfus0' ; xlab = 'Beam_bkgrd_multiplier' case (8) fqval = swp(iscn) vlab = 'fqval' ; xlab = 'Big_Q_f-value' case (9) te = swp(iscn) vlab = 'te' ; xlab = 'Electron_temperature_keV' case (10) boundu(15) = swp(iscn) vlab = 'boundu(15)' ; xlab = 'Volt-second_upper_bound' case (11) dnbeta = swp(iscn) vlab = 'dnbeta' ; xlab = 'Beta_coefficient' case (12) bootstrap_current_fraction_max = swp(iscn) vlab = 'bootstrap_current_fraction_max' ; xlab = 'Bootstrap_fraction' case (13) boundu(10) = swp(iscn) vlab = 'boundu(10)' ; xlab = 'H_factor_upper_bound' case (14) fiooic = swp(iscn) vlab = 'fiooic' ; xlab = 'TFC_Iop_/_Icrit_f-value' case (15) fjprot = swp(iscn) vlab = 'fjprot' ; xlab = 'TFC_Jprot_limit_f-value' case (16) rmajor = swp(iscn) vlab = 'rmajor' ; xlab = 'Plasma_major_radius_(m)' case (17) bmxlim = swp(iscn) vlab = 'bmxlim' ; xlab = 'Max_toroidal_field_(T)' case (18) gammax = swp(iscn) vlab = 'gammax' ; xlab = 'Maximum_CD_gamma' case (19) boundl(16) = swp(iscn) vlab = 'boundl(16)' ; xlab = 'CS_thickness_lower_bound' case (20) t_burn_min = swp(iscn) vlab = 't_burn_min' ; xlab = 'Minimum_burn_time_(s)' case (21) ! sigpfalw = swp(iscn) vlab = 'obsolete' ; xlab = 'obsolete' case (22) if (iavail == 1) call report_error(95) cfactr = swp(iscn) vlab = 'cfactr' ; xlab = 'Plant_availability_factor' case (23) boundu(72) = swp(iscn) vlab = 'boundu(72)' ; xlab = 'Ip/Irod_upper_bound' case (24) powfmax = swp(iscn) vlab = 'powfmax' ; xlab = 'Fusion_power_limit_(MW)' case (25) kappa = swp(iscn) vlab = 'kappa' ; xlab = 'Plasma_elongation' case (26) triang = swp(iscn) vlab = 'triang' ; xlab = 'Plasma_triangularity' case (27) tbrmin = swp(iscn) vlab = 'tbrmin' ; xlab = 'Min_tritium_breed._ratio' case (28) bt = swp(iscn) vlab = 'bt' ; xlab = 'Tor._field_on_axis_(T)' case (29) coreradius = swp(iscn) vlab = 'coreradius' ; xlab = 'Core_radius' case (30) !fimpvar = swp(iscn) vlab = 'OBSOLETE' ; xlab = 'OBSOLETE' case (31) taulimit = swp(iscn) vlab = 'taulimit' ; xlab = 'Taup/taueff_lower_limit' case (32) epsvmc = swp(iscn) vlab = 'epsvmc' ; xlab = 'VMCON error tolerance' case (33, 34, 35, 36, 37, 47) write(*,*) 'Kallenbach model has been removed, remove the kallenbach scan variables' stop 1 case (38) boundu(129) = swp(iscn) vlab = 'boundu(129)' ; xlab = ' Neon upper limit' case (39) boundu(131) = swp(iscn) vlab = 'boundu(131)' ; xlab = ' Argon upper limit' case (40) boundu(135) = swp(iscn) vlab = 'boundu(135)' ; xlab = ' Xenon upper limit' case (41) blnkoth = swp(iscn) vlab = 'blnkoth' ; xlab = 'Outboard blanket thick.' case (42) fimp(9) = swp(iscn) impurity_arr_frac(9) = fimp(9) vlab = 'fimp(9)' ; xlab = 'Argon fraction' case (43) ! rho_ecrh = swp(iscn) vlab = 'obsolete' ; xlab = 'obsolete' case (44) sig_tf_case_max = swp(iscn) vlab = 'sig_tf_case_max' ; xlab = 'Allowable_stress_in_tf_coil_case_Tresca_(pa)' case (45) tmargmin_tf = swp(iscn) vlab = 'tmargmin_tf' ; xlab = 'Minimum_allowable_temperature_margin' case (46) boundu(152) = swp(iscn) vlab = 'boundu(152)' ; xlab = 'Max allowable fgwsep' case (48) n_pancake = int(swp(iscn)) vlab = 'n_pancake' ; xlab = 'TF Coil - n_pancake' case (49) n_layer = int(swp(iscn)) vlab = 'n_layer' ; xlab = 'TF Coil - n_layer' case (50) fimp(13) = swp(iscn) impurity_arr_frac(13) = fimp(13) vlab = 'fimp(13)' ; xlab = 'Xenon fraction' case (51) ftar = swp(iscn) vlab = 'ftar' ; xlab = 'lower_divertor_power_fraction' case (52) rad_fraction_sol = swp(iscn) vlab = 'rad_fraction_sol' ; xlab = 'SoL radiation fraction' case (53) boundu(157) = swp(iscn) vlab = 'boundu(157)' ; xlab = 'Max allowable fvssu' case (54) b_crit_upper_nbti = swp(iscn) vlab = 'Bc2(0K)' ; xlab = 'GL_NbTi Bc2(0K)' case(55) shldith = swp(iscn) vlab = 'shldith' ; xlab = 'Inboard neutronic shield' case(56) crypmw_max = swp(iscn) vlab = 'crypmw_max' ; xlab = 'max allowable crypmw' case(57) boundl(2) = swp(iscn) vlab = 'boundl(2)' ; xlab = 'bt minimum' case(58) scrapli = swp(iscn) vlab = 'scrapli' ; xlab = 'Inboard FW-plasma sep gap' case(59) scraplo = swp(iscn) vlab = 'scraplo' ; xlab = 'Outboard FW-plasma sep gap' case (60) sig_tf_wp_max = swp(iscn) vlab = 'sig_tf_wp_max' ; xlab = 'Allowable_stress_in_tf_coil_conduit_Tresca_(pa)' case (61) copperaoh_m2_max = swp(iscn) vlab = 'copperaoh_m2_max' ; xlab = 'Max CS coil current / copper area' case (62) coheof = swp(iscn) vlab = 'coheof' ; xlab = 'CS coil current density at EOF (A/m2)' case (63) ohcth = swp(iscn) vlab = 'ohcth' ; xlab = 'CS coil thickness (m)' case (64) ohhghf = swp(iscn) vlab = 'ohhghf' ; xlab = 'CS height (m)' case (65) n_cycle_min = swp(iscn) vlab = 'n_cycle_min' ; xlab = 'CS stress cycles min' case (66) oh_steel_frac = swp(iscn) vlab = 'oh_steel_frac' ; xlab = 'CS steel fraction' case (67) t_crack_vertical = swp(iscn) vlab = 't_crack_vertical' ; xlab = 'Initial crack vertical size (m)' case (68) inlet_temp_liq = swp(iscn) vlab = 'inlet_temp_liq' ; xlab = 'Inlet Temperature Liquid Metal Breeder/Coolant (K)' case (69) outlet_temp_liq = swp(iscn) vlab = 'outlet_temp_liq' ; xlab = 'Outlet Temperature Liquid Metal Breeder/Coolant (K)' case(70) blpressure_liq = swp(iscn) vlab = 'blpressure_liq' ; xlab = 'Blanket liquid metal breeder/coolant pressure (Pa)' case(71) n_liq_recirc = swp(iscn) vlab = 'n_liq_recirc' ; xlab = 'Selected number of liquid metal breeder recirculations per day' case(72) bz_channel_conduct_liq = swp(iscn) vlab = 'bz_channel_conduct_liq' ; xlab = 'Conductance of liquid metal breeder duct walls (A V-1 m-1)' case(73) pnuc_fw_ratio_dcll = swp(iscn) vlab = 'pnuc_fw_ratio_dcll' ; xlab = 'Ratio of FW nuclear power as fraction of total (FW+BB)' case(74) f_nuc_pow_bz_struct = swp(iscn) vlab = 'f_nuc_pow_bz_struct' ; xlab = 'Fraction of BZ power cooled by primary coolant for dual-coolant balnket' case(75) pitch = swp(iscn) vlab = 'pitch' ; xlab = 'pitch of first wall cooling channels (m)' case (76) etath = swp(iscn) vlab = 'etath' ; xlab = 'Thermal conversion eff.' case (77) startupratio = swp(iscn) vlab = 'startupratio' ; xlab = 'Gyrotron redundancy' case (78) fkind = swp(iscn) vlab = 'fkind' ; xlab = 'Multiplier for Nth of a kind costs' case (79) etaech = swp(iscn) vlab = 'etaech' ; xlab = 'ECH wall plug to injector efficiency' case (80) fcoolcp = swp(iscn) vlab = 'fcoolcp' ; xlab = 'Coolant fraction of TF' case (81) n_tf_turn = swp(iscn) vlab = 'n_tf_turn' ; xlab = 'Number of turns in TF' case default idiags(1) = nwp ; call report_error(96) end select end subroutine scan_select subroutine post_optimise(ifail) !! Called after calling the optimising equation solver from Python. !! author: P J Knight, CCFE, Culham Science Centre !! ifail : input integer : error flag !! use constraints use error_handling use numerics use process_output use utilities, only:upper_case use main_module, only:verror ! for ipedestal = 2 option use global_variables, only: convergence_parameter use constants, only: iotty, nout, mfile use physics_variables, only: ipedestal use define_iteration_variables, only: boundxc, loadxc implicit none ! Arguments integer, intent(in) :: ifail ! Local variables integer :: ii,inn,iflag real(dp) :: summ,xcval,xmaxx,xminn,f,xnorm real(dp), dimension(ipeqns) :: con1, con2, err character(len=1), dimension(ipeqns) :: sym character(len=10), dimension(ipeqns) :: lab character(len=60) :: string1, string2 ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Check on accuracy of solution by summing the ! squares of the residuals of the equality constraints summ = 0.0D0 do ii = 1,neqns summ = summ + rcm(ii)*rcm(ii) end do sqsumsq = sqrt(summ) ! Turn on error reporting errors_on = .true. ! Print out information on solution call oheadr(nout,'Numerics') call ocmmnt(nout,'PROCESS has performed a VMCON (optimisation) run.') if (ifail /= 1) then !call ocmmnt(nout,'but could not find a feasible set of parameters.') ! call oheadr(nout,'PROCESS COULD NOT FIND A FEASIBLE SOLUTION') ! call ovarin(iotty,'VMCON error flag (ifail)','',ifail) call ovarin(nout,'VMCON error flag','(ifail)',ifail) call oheadr(iotty,'PROCESS COULD NOT FIND A FEASIBLE SOLUTION') call oblnkl(iotty) idiags(1) = ifail ; call report_error(132) else call ocmmnt(nout,'and found a feasible set of parameters.') call oblnkl(nout) call ovarin(nout,'VMCON error flag','(ifail)',ifail) call oheadr(iotty,'PROCESS found a feasible solution') end if !call oblnkl(nout) ! If necessary, write out a relevant error message if (ifail /= 1) then call verror(ifail) call oblnkl(nout) call oblnkl(iotty) else ! Show a warning if the constraints appear high even if allegedly converged if (sqsumsq >= 1.0D-2) then call oblnkl(nout) call ocmmnt(nout,'WARNING: Constraint residues are HIGH; consider re-running') call ocmmnt(nout,' with lower values of EPSVMC to confirm convergence...') call ocmmnt(nout,' (should be able to get down to about 1.0E-8 okay)') call oblnkl(nout) call ocmmnt(iotty,'WARNING: Constraint residues are HIGH; consider re-running') call ocmmnt(iotty,' with lower values of EPSVMC to confirm convergence...') call ocmmnt(iotty,' (should be able to get down to about 1.0E-8 okay)') call oblnkl(iotty) fdiags(1) = sqsumsq ; call report_error(134) end if end if call ovarin(nout,'Number of iteration variables','(nvar)',nvar) call ovarin(nout,'Number of constraints (total)','(neqns+nineqns)',neqns+nineqns) call ovarin(nout,'Optimisation switch','(ioptimz)',ioptimz) call ovarin(nout,'Figure of merit switch','(minmax)',minmax) ! if (ifail /= 1) then ! call ovarin(nout,'VMCON error flag','(ifail)',ifail) ! end if objf_name = '"'//trim(lablmm(abs(minmax)))//'"' ! Quotes required for string parsing in MFILE call ovarst(nout,'Objective function name','(objf_name)',objf_name) call ovarre(nout,'Normalised objective function','(norm_objf)',norm_objf, 'OP ') call ovarre(nout,'Square root of the sum of squares of the constraint residuals','(sqsumsq)',sqsumsq, 'OP ') call ovarre(nout,'VMCON convergence parameter','(convergence_parameter)',convergence_parameter, 'OP ') call ovarin(nout,'Number of VMCON iterations','(nviter)',nviter, 'OP ') call oblnkl(nout) if (ifail == 1) then string1 = 'PROCESS has successfully optimised the iteration variables' else string1 = 'PROCESS has tried to optimise the iteration variables' end if if (minmax > 0) then string2 = ' to minimise the figure of merit: ' else string2 = ' to maximise the figure of merit: ' end if call upper_case(objf_name) write(nout,10) trim(string1) // trim(string2), trim(objf_name) 10 format(a90, t92, a22) call oblnkl(nout) ! Check which variables are at bounds iflag = 0 do ii = 1,nvar xminn = 1.01D0*bondl(ii) xmaxx = 0.99D0*bondu(ii) if (xcm(ii) < xminn) then if (iflag == 0) then call ocmmnt(nout, & 'Certain operating limits have been reached,') call ocmmnt(nout, & 'as shown by the following iteration variables that are') call ocmmnt(nout, & 'at or near to the edge of their prescribed range :') call oblnkl(nout) iflag = 1 end if xcval = xcm(ii)*scafc(ii) !write(nout,30) ii,lablxc(ixc(ii)),xcval,bondl(ii)*scafc(ii) write(nout,30) lablxc(ixc(ii)),xcval,bondl(ii)*scafc(ii) end if if (xcm(ii) > xmaxx) then if (iflag == 0) then call ocmmnt(nout, & 'Certain operating limits have been reached,') call ocmmnt(nout, & 'as shown by the following iteration variables that are') call ocmmnt(nout, & 'at or near to the edge of their prescribed range :') call oblnkl(nout) iflag = 1 end if xcval = xcm(ii)*scafc(ii) write(nout,40) lablxc(ixc(ii)),xcval,bondu(ii)*scafc(ii) end if end do !30 format(t4,'Variable ',i3,' (',a9, & ! ',',1pe12.4,') is at or below its lower bound:',1pe12.4) 30 format(t4, a30, '=',1pe12.4,' is at or below its lower bound:',1pe12.4) 40 format(t4, a30, '=',1pe12.4,' is at or above its upper bound:',1pe12.4) !40 format(t4,'Variable ',i3,' (',a9, & ! ',',1pe12.4,') is at or above its upper bound:',1pe12.4) ! Print out information on numerics call osubhd(nout,'The solution vector is comprised as follows :') ! write(nout,50) ! Remove Lagrange multipliers as no-one understands them. ! MFILE not changed !50 format(t47,'lower',t59,'upper') write(nout,60) !60 format(t23,'final',t33,'fractional',t46,'Lagrange',t58,'Lagrange') 60 format(t43,'final',t55,'final /') write(nout,70) !70 format(t5,'i',t23,'value',t35,'change',t45,'multiplier', & ! t57,'multiplier') 70 format(t5,'i',t43,'value',t55,'initial') call oblnkl(nout) do inn = 1,nvar xcs(inn) = xcm(inn)*scafc(inn) ! write(nout,80) inn,lablxc(ixc(inn)),xcs(inn),xcm(inn), & ! vlam(neqns+nineqns+inn), vlam(neqns+nineqns+1+inn+nvar) write(nout,80) inn,lablxc(ixc(inn)),xcs(inn),xcm(inn) !80 format(t2,i4,t8,a9,t19,4(1pe12.4)) !80 format(t2,i4,t8,a30,t39,2(1pe12.4)) 80 format(t2,i4,t8,a30,t39,1pe12.4, t52, 0pf10.4) ! MDK The 0p is needed because of a bizarre "feature"/bug in fortran: ! the 1p in the previous format continues until changed. call ovarre(mfile,lablxc(ixc(inn)),'(itvar'//int_to_string3(inn)//')',xcs(inn)) ! 'Range-normalised' iteration variable values for MFILE: ! 0.0 (at lower bound) to 1.0 (at upper bound) if (bondl(inn) == bondu(inn)) then xnorm = 1.0D0 else xnorm = (xcm(inn) - bondl(inn)) / (bondu(inn) - bondl(inn)) xnorm = max(xnorm, 0.0D0) xnorm = min(xnorm, 1.0D0) end if ! Added ratio final/initial to MFILE call ovarre(mfile,trim(lablxc(ixc(inn)))//' (final value/initial value)', & '(xcm'//int_to_string3(inn)//')',xcm(inn)) call ovarre(mfile,trim(lablxc(ixc(inn)))//' (range normalised)', & '(nitvar'//int_to_string3(inn)//')',xnorm) end do call osubhd(nout, & 'The following equality constraint residues should be close to zero :') call constraint_eqns(neqns+nineqns,-1,con1,con2,err,sym,lab) write(nout,90) 90 format(t48,'physical',t73,'constraint',t100,'normalised') write(nout,100) 100 format(t47,'constraint',t74,'residue',t101,'residue') call oblnkl(nout) do inn = 1,neqns write(nout,110) inn,lablcc(icc(inn)),sym(inn),con2(inn), & lab(inn),err(inn),lab(inn),con1(inn) call ovarre(mfile,lablcc(icc(inn))//' normalised residue', & '(eq_con'//int_to_string3(icc(inn))//')',con1(inn)) end do 110 format(t2,i4,t8,a33,t46,a1,t47,1pe12.4,t60,a10,t71,1pe12.4,t84,a10,t98,1pe12.4) if (nineqns > 0) then call osubhd(nout, & 'The following inequality constraint residues should be greater than or approximately equal to zero :') do inn = neqns+1,neqns+nineqns !write(nout,120) inn,lablcc(icc(inn)),rcm(inn),vlam(inn) write(nout,110) inn,lablcc(icc(inn)),sym(inn),con2(inn), & lab(inn), err(inn), lab(inn) call ovarre(mfile,lablcc(icc(inn)),'(ineq_con'//int_to_string3(icc(inn))//')',rcm(inn)) end do end if ! 120 format(t2,i4,t8,a33,t45,1pe12.4,1pe12.4) end subroutine post_optimise ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! end module scan_module