3548 implement variable dt ratio

documentation/proc-pages/physics-models/fusion_reactions/plasma_reactions.md

@@ -11,7 +11,7 @@ $$
 
 Roughly 20% of the energy produced is given to the alpha particles (\(^4\)He). The remaining 80% is carried
 away by the neutrons, which deposit their energy within the blanket and shield and other reactor components.
-The fraction of the alpha energy deposited in the plasma is [`f_alpha_plasma`](#coupled-alpha-particle-power). 
+The fraction of the alpha energy deposited in the plasma is [`f_alpha_plasma`](#coupled-alpha-particle-power).
 **`PROCESS` only assumes the alpha power produced is coupled to and self heats the plasma, other charged particles do not.**
 
 PROCESS can also model D-\(^3\)He power plants, which utilise the following
@@ -38,11 +38,11 @@ $$\begin{aligned}
 \mathrm{D + D}  & \Longrightarrow \mathrm{T + p + 4.03 \,MeV}
 \end{aligned}$$
 
-Also, as tritium is produced by the latter reaction, D-T fusion also occurs. 
+Also, as tritium is produced by the latter reaction, D-T fusion also occurs.
 As a result, there is still a small amount of neutron power
 extracted from the plasma.
 
-Pure D-\(^3\)He tokamak power plants do not include breeding blankets, because 
+Pure D-\(^3\)He tokamak power plants do not include breeding blankets, because
 no tritium needs to be produced for fuel.
 
 The contributions from all four of the above fusion reactions are included in
@@ -56,7 +56,9 @@ profiles
 
 The fractional composition of the 'fuel' ions ($\text{D}$, $\text{T}$ and $^3\text{He}$) is
 controlled using the three variables `f_deuterium`, `f_tritium` and `f_helium3`, respectively.
-More information about setting seeded impurities and simulating first wall sputtering can be found in the [composition and impurities section](../plasma_composition.md)
+More information about setting seeded impurities and simulating first wall sputtering can be found in the [composition and impurities section](../plasma_composition.md).
+
+It is also possible to optimise on the deuterium-tritium fuel mixture ratio. For this, `f_tritium` **must** be set as an iteration variable with `ixc = 173`. More info can be found [here](#fuel-ions-mixture-self-consistency).
 
 !!! note "Reactions not calculated"
 
@@ -232,5 +234,13 @@ This constraint can be activated by stating `icc = 28` in the input file.
 
 The value of `bigqmin` can be set to the minimum desired $Q_{\text{plasma}}$ value. The scaling value `fqval` can be varied also.
 
+-------------------------
+
+### Fuel ions mixture self-consistency
+
+This constraint can be activated by stating `icc = 92` in the input file.
+
+This constraint should be activated if the user wishes to allow the ratio of the fuel ions to be iterated upon. For this, `f_tritium` has to be set as an iteration variable with `ixc = 173`. This constraint ensure that the summation of the fuel fractions always sums to 1.0 and will modify the other compositions depending on the value of `f_tritium`.
+
 [^1]: H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and thermal reactivities,” Nuclear Fusion, vol. 32, no. 4, pp. 611–631, Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07.
-[^2]: I. P. E. G. on E. Drive and I. P. B. Editors, “Chapter 5: Physics of energetic ions,” Nuclear Fusion, vol. 39, no. 12, pp. 2471–2495, Dec. 1999, doi: https://doi.org/10.1088/0029-5515/39/12/305.
+[^2]: I. P. E. G. on E. Drive and I. P. B. Editors, “Chapter 5: Physics of energetic ions,” Nuclear Fusion, vol. 39, no. 12, pp. 2471–2495, Dec. 1999, doi: https://doi.org/10.1088/0029-5515/39/12/305.

source/fortran/constraint_equations.f90

@@ -286,6 +286,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 D/T ratio
+        case (92); call constraint_eqn_092(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units)
        case default
 
         idiags(1) = icc(i)
@@ -3387,6 +3389,32 @@ 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 checking is D/T ratio is consistent, and sums to 1.
+      !! author: G Turkington, UKAEA
+      !! args : output structure : residual error; constraint value;
+      !! residual error in physical units; output string; units string
+      !! f_deuterium : input : fraction of deuterium ions
+      !! f_tritium  : input : fraction of tritium ions
+      !! f_helium3  : input : fraction of helium-3 ions
+      use physics_variables, only: f_deuterium, f_tritium, f_helium3
+      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
+
+
+      ! Iterate over f_tritium and calculate f_deuterium
+      f_deuterium = 1.0D0 - (f_tritium + f_helium3)
+      tmp_cc = 1.0D0 - (f_deuterium + f_tritium + f_helium3)
+      tmp_con = 1.0D0
+      tmp_err = tmp_con * tmp_cc
+      tmp_symbol = '='
+      tmp_units = 'fraction'
+
+   end subroutine constraint_eqn_092
 
 
 end module constraints

source/fortran/iteration_variables.f90

@@ -3850,31 +3850,28 @@ 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) f_tritium : Tritium fraction in fuel
     use numerics, only: lablxc, boundl, boundu
     implicit none
-    lablxc(173) = 'DUMMY         '
-    boundl(173) = 1.0d-99
-    boundu(173) = 1.0d99
+    lablxc(173) = 'f_tritium   '
+    boundl(173) = 0.000
+    boundu(173) = 1.000
   end subroutine init_itv_173
 
   real(kind(1.d0)) function itv_173()
-    implicit none
-    itv_173 = DUMMY
+    use physics_variables, only: f_tritium
+    itv_173 = f_tritium
   end function itv_173
 
   subroutine set_itv_173(ratio)
-    implicit none
+    use physics_variables, only: f_tritium
     real(kind(1.d0)) :: ratio
-    DUMMY = ratio
+    f_tritium = ratio
   end subroutine set_itv_173
 
+  !---------------------------------
+  ! DUMMY variables below here
   !---------------------------------
 
   subroutine init_itv_174
@@ -4117,8 +4114,8 @@ subroutine loadxc
            case (170);  xcm(i) = itv_170()
            case (171);  xcm(i) = itv_171()
            case (172);  xcm(i) = itv_172()
-            ! DUMMY Cases
            case (173);  xcm(i) = itv_173()
+            ! DUMMY Cases
            case (174);  xcm(i) = itv_174()
            case (175);  xcm(i) = itv_175()
 

source/fortran/numerics.f90

@@ -19,7 +19,7 @@ module numerics
   integer, parameter :: ipnvars = 175
   !!  ipnvars FIX : total number of variables available for iteration
 
-  integer, parameter :: ipeqns = 91
+  integer, parameter :: ipeqns = 92
   !!  ipeqns  FIX : number of constraint equations available
 
   integer, parameter :: ipnfoms = 19
@@ -193,7 +193,9 @@ module numerics
   !!  <LI> (88) Constraint for TF coil strain absolute value
   !!  <LI> (89) Constraint for CS coil quench protection
   !!  <LI> (90) Lower Limit on number of stress load cycles for CS (itr 167 fncycle)
-  !!  <LI> (91) Checking if the design point is ECRH ignitable (itv 168 fecrh_ignition) </UL>
+  !!  <LI> (91) Checking if the design point is ECRH ignitable (itv 168 fecrh_ignition)
+  !!  <LI> (92) D/T/He3 ratio in fuel sums to 1 </UL>
+
 
   integer, dimension(ipnvars) :: ixc
   !!  ixc(ipnvars) /0/ :
@@ -374,7 +376,7 @@ module numerics
   !! <LI> (170) beta_div : field line angle wrt divertor target plate (degrees)
   !! <LI> (171) casths_fraction : TF side case thickness as fraction of toridal case thickness
   !! <LI> (172) casths : TF side case thickness [m]
-  !! <LI> (173) EMPTY : Description
+  !! <LI> (173) f_deuterium : Deuterium fraction in fuel
   !! <LI> (174) EMPTY : Description
   !! <LI> (175) EMPTY : Description
   ! Issue 287 iteration variables are now defined in module define_iteration_variables in iteration variables.f90
@@ -545,7 +547,9 @@ subroutine init_numerics()
       'TF coil strain absolute value    ', &
       'CS current/copper area < Max     ', &
       'CS stress load cycles            ', &
-      'ECRH ignitability                '  &
+      'ECRH ignitability                ', &
+      'Fuel composition consistency     ' &
+
       /)
 
     ! Please note: All strings between '...' above must be exactly 33 chars long

tests/unit/test_power.py

@@ -952,6 +952,7 @@ class PfpwrParam(NamedTuple):
                 False,
                 False,
                 False,
+                False,
             ),
             ioptimz=1,
             tim=np.array(
@@ -1694,6 +1695,7 @@ class PfpwrParam(NamedTuple):
                 False,
                 False,
                 False,
+                False,
             ),
             ioptimz=1,
             tim=np.array(

tests/unit/test_pulse.py

@@ -647,6 +647,7 @@ class BurnParam(NamedTuple):
                 False,
                 False,
                 False,
+                False,
             ),
             i_pulsed_plant=1,
             outfile=11,
@@ -1204,6 +1205,7 @@ class BurnParam(NamedTuple):
                 False,
                 False,
                 False,
+                False,
             ),
             i_pulsed_plant=1,
             outfile=11,