🚧 3273 calculated bootstrap current fraction is overridden

documentation/proc-pages/physics-models/plasma_current.md

@@ -62,6 +62,10 @@ methods, as summarised here. Note that methods `ibss = 1-3` do not take into acc
 existence of pedestals, whereas the Sauter et al. scaling 
 (`ibss = 4`) allows general profiles to be used. 
 
+!!! Warning "Bootstrap Current Enforcement"
+    The constraint equation `icc = 92` is responsible for making sure the calculated plasma current is whithin the user's desired bounds. `bscfmax` sets the maximum desired value of the bootstrap fraction. `fboot_max` is the f-value to scale this input for the constraint equation. 
+    It is not uncommon for the intial conditions of an input file to give a bootstrap current fraction greater than one for some of the models. It is recommended to then set initially wide bounds of the associated variables. If this fails, the bootstrap current can be set directly by setting `bscfmax` to the negative value of the required fraction. This will cause the bootstrap fraction models to be ignored in the rest of the optimisation.
+
 | `ibss` | Description |
 | :-: | - |
 | 1 | ITER scaling -- To use the ITER scaling method for the bootstrap current fraction.  Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). This method is valid at high aspect ratio only.

process/physics.py

@@ -1023,13 +1023,6 @@ def physics(self):
                 error_handling.idiags[0] = physics_variables.ibss
                 error_handling.report_error(75)
 
-            physics_module.err242 = 0
-            if current_drive_variables.bootipf > current_drive_variables.bscfmax:
-                current_drive_variables.bootipf = min(
-                    current_drive_variables.bootipf, current_drive_variables.bscfmax
-                )
-                physics_module.err242 = 1
-
             if physics_variables.idia == 1:
                 current_drive_variables.diaipf = current_drive_variables.diacf_hender
             elif physics_variables.idia == 2:
@@ -4110,9 +4103,6 @@ def outplas(self):
                 current_drive_variables.pscf_scene,
                 "OP ",
             )
-            # Error to catch if bootstap fraction limit has been enforced
-            if physics_module.err242 == 1:
-                error_handling.report_error(242)
 
             # Error to catch if self-driven current fraction limit has been enforced
             if physics_module.err243 == 1:

source/fortran/constraint_equations.f90

@@ -287,6 +287,8 @@ subroutine constraint_eqns(m,ieqn,cc,con,err,symbol,units)
         case (90); call constraint_eqn_090(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
          ! Constraint for indication of ECRH ignitability
         case (91); call constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+         ! Constraint for bootstrap current fraction
+        case (92); call constraint_eqn_092(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) 
        case default
 
         idiags(1) = icc(i)
@@ -3390,6 +3392,29 @@ subroutine constraint_eqn_091(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
       tmp_units = 'MW'
    end subroutine constraint_eqn_091
 
-
+   subroutine constraint_eqn_092(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
+      !! Equation for bootstrap current fraction upper limit
+      !! author: C Ashe, UKAEA, Culham Campus
+      !! args : output structure : residual error; constraint value;
+      !! residual error in physical units; output string; units string
+      !! Logic change during pre-factoring: err, symbol, units will be assigned only if present.
+      !! bootipf : input real : Bootstrap current fraction
+      !! fboot_max : input real : Desired maximum bootstrap current fraction
+      !! bscfmax : input real : f-value for constraint bootipf < fboot_max
+      use current_drive_variables, only: bootipf, bscfmax
+      use constraint_variables, only: fboot_max
+      implicit none
+            real(dp), intent(out) :: tmp_cc
+      real(dp), intent(out) :: tmp_con
+      real(dp), intent(out) :: tmp_err
+      character(len=1), intent(out) :: tmp_symbol
+      character(len=10), intent(out) :: tmp_units
+
+      tmp_cc =  1.0D0 - fboot_max*bscfmax/bootipf
+      tmp_con = bscfmax
+      tmp_err = bscfmax - bootipf / fboot_max
+      tmp_symbol = '<'
+      tmp_units = ''
+   end subroutine constraint_eqn_092
 
 end module constraints

source/fortran/constraint_variables.f90

@@ -54,6 +54,9 @@ module constraint_variables
   !! f-value for density limit (`constraint equation 5`, `iteration variable 9`)
   !! (invalid if `ipedestal=3`)
 
+  real(dp) :: fboot_max
+  !! f-value for bootstrap current fraction limit (`constraint equation 92`, `iteration variable 173`)
+
   real(dp) :: fdivcol
   !! f-value for divertor collisionality (`constraint equation 22`, `iteration variable 34`)
 

source/fortran/input.f90

@@ -230,7 +230,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
       workshop_l, workshop_w, workshop_h
     use constraint_variables, only: flhthresh, fpeakb, fpsep, fdivcol, ftcycl, &
       betpmx, fpsepbqar, ftmargtf, fradwall, fptfnuc, fnesep, fportsz, tbrmin, &
-      maxradwallload, pseprmax, fdene, fniterpump, fpinj, pnetelin, powfmax, &
+      maxradwallload, pseprmax, fdene, fboot_max, fniterpump, fpinj, pnetelin, powfmax, &
       fgamcd, ftbr, mvalim, taulimit, walalw, fmva, fradpwr, nflutfmax, fipir, &
       fauxmn, fiooic, fcwr, fjohc0, frminor, psepbqarmax, ftpeak, bigqmin, &
       fstrcond, fptemp, ftmargoh, fvs, fbetatry, vvhealw, fpnetel, ftburn, &
@@ -808,6 +808,9 @@ subroutine parse_input_file(in_file,out_file,show_changes)
        case ('fdene')
           call parse_real_variable('fdene', fdene, 0.001D0, 10.0D0, &
                'F-value for density limit')
+       case ('fboot_max')
+          call parse_real_variable('fboot_max', fboot_max, 0.001D0, 10.0D0, &
+               'F-value for bootstrap current fraction limit')       
        case ('fdivcol')
           call parse_real_variable('fdivcol', fdivcol, 0.001D0, 10.0D0, &
                'F-value for divertor collisionality')

source/fortran/iteration_variables.f90

@@ -3850,31 +3850,34 @@ subroutine set_itv_172(ratio)
     casths = ratio
   end subroutine set_itv_172
 
-  !---------------------------------
-  ! DUMMY variables below here
-  !---------------------------------
+
   !---------------------------------
 
   subroutine init_itv_173
-    !! <LI> (173) DUMMY : Description
+    !! <LI> (173) fboot_max : F valie for icc = 92
     use numerics, only: lablxc, boundl, boundu
     implicit none
-    lablxc(173) = 'DUMMY         '
-    boundl(173) = 1.0d-99
-    boundu(173) = 1.0d99
+    lablxc(173) = 'fboot_max         '
+    boundl(173) = 0.001D0
+    boundu(173) = 10.000D0
   end subroutine init_itv_173
 
   real(kind(1.d0)) function itv_173()
+    use constraint_variables, only: fboot_max
     implicit none
-    itv_173 = DUMMY
+    itv_173 = fboot_max
   end function itv_173
 
   subroutine set_itv_173(ratio)
+    use constraint_variables, only: fboot_max
     implicit none
     real(kind(1.d0)) :: ratio
-    DUMMY = ratio
+    fboot_max = ratio
   end subroutine set_itv_173
-
+
+  !---------------------------------
+  ! DUMMY variables below here
+  !---------------------------------
   !---------------------------------
 
   subroutine init_itv_174
@@ -4116,10 +4119,10 @@ subroutine loadxc
            case (168);  xcm(i) = itv_168()
            case (169);  xcm(i) = itv_169()
            case (170);  xcm(i) = itv_170()
-           case (171);  xcm(i) = itv_171()
-           case (172);  xcm(i) = itv_172()
+           case (171);  xcm(i) = itv_171() 
+           case (172);  xcm(i) = itv_172() 
+           case (173);  xcm(i) = itv_173() 
             ! DUMMY Cases
-           case (173);  xcm(i) = itv_173()
            case (174);  xcm(i) = itv_174()
            case (175);  xcm(i) = itv_175()
 
@@ -4384,9 +4387,9 @@ subroutine convxc(xc,nn)
            case (169);  call set_itv_169(ratio)
            case (170);  call set_itv_170(ratio)
            case (171);  call set_itv_171(ratio)
-           case (172);  call set_itv_172(ratio)
-            ! DUMMY Cases
+           case (172);  call set_itv_172(ratio) 
            case (173);  call set_itv_173(ratio)
+            ! DUMMY Cases
            case (174);  call set_itv_174(ratio)
            case (175);  call set_itv_175(ratio)
 

source/fortran/physics.f90

@@ -17,7 +17,7 @@ module physics_module
   !  Module-level variables
 
   integer :: iscz
-  integer :: err242, err243
+  integer :: err243
   real(dp) :: rad_fraction_LCFS
   real(dp) :: total_plasma_internal_energy  ! [J]
   real(dp) :: total_loss_power        ! [W]
@@ -41,7 +41,6 @@ subroutine init_physics_module
 
     first_call = 1
     iscz = 0
-    err242 = 0
     err243 = 0
     rad_fraction_LCFS = 0.0D0
     total_plasma_internal_energy = 0.0D0

tests/integration/ref_dicts.json

@@ -1758,7 +1758,6 @@
         "epsfcn": 0.001,
         "epsilon0": 8.85418781e-12,
         "epsvmc": 1e-06,
-        "err242": 0.0,
         "err243": 0.0,
         "error": ".False.",
         "error_head": "> null()",
@@ -9389,7 +9388,6 @@
         "epsfcn": "epsfcn /1.0e-3/ : finite difference step length for HYBRD/VMCON derivatives",
         "epsilon0": "permittivity of free space [Farad/m]",
         "epsvmc": "epsvmc /1.0e-6/ : error tolerance for VMCON",
-        "err242": "",
         "err243": "",
         "error": "",
         "error_head": "",
@@ -19005,7 +19003,6 @@
         ],
         "physics_module": [
             "iscz",
-            "err242",
             "err243",
             "rad_fraction_LCFS",
             "total_plasma_internal_energy",

tests/unit/test_physics.py

@@ -890,8 +890,6 @@ class PlasmaCompositionParam(NamedTuple):
 
     iscz: Any = None
 
-    err242: Any = None
-
     err243: Any = None
 
     ptarmw: Any = None
@@ -1016,7 +1014,6 @@ class PlasmaCompositionParam(NamedTuple):
             dene=7.5e19,
             dnprot=0,
             iscz=0,
-            err242=0,
             err243=0,
             ptarmw=0,
             lambdaio=0,
@@ -1138,7 +1135,6 @@ class PlasmaCompositionParam(NamedTuple):
             dene=7.5e19,
             dnprot=7500000000000000,
             iscz=0,
-            err242=0,
             err243=0,
             ptarmw=33.990985729118783,
             lambdaio=0.00157,
@@ -1290,8 +1286,6 @@ def test_plasma_composition(plasmacompositionparam, monkeypatch, physics):
 
     monkeypatch.setattr(physics_module, "iscz", plasmacompositionparam.iscz)
 
-    monkeypatch.setattr(physics_module, "err242", plasmacompositionparam.err242)
-
     monkeypatch.setattr(physics_module, "err243", plasmacompositionparam.err243)
 
     monkeypatch.setattr(physics_module, "ptarmw", plasmacompositionparam.ptarmw)