diff --git a/documentation/proc-pages/physics-models/plasma_current.md b/documentation/proc-pages/physics-models/plasma_current.md index a1395a0edd..1a2a3f3029 100644 --- a/documentation/proc-pages/physics-models/plasma_current.md +++ b/documentation/proc-pages/physics-models/plasma_current.md @@ -68,6 +68,7 @@ existence of pedestals, whereas the Sauter et al. scaling | 2 | General scaling -- To use a more general scaling method, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). | 3 | Numerically fitted scaling [^8] -- To use a numerically fitted scaling method, valid for all aspect ratios, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). | 4 | Sauter, Angioni and Lin-Liu scaling [^9] [^10] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). +| 5 | Sakai, Fujita and Okamoto scaling [^11] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). The model includes the toroidal diamagnetic current in the calculation due to the dataset, so `idia = 0` can only be used with it !!! Note "Fixed Bootstrap Current" Direct input -- To input the bootstrap current fraction directly, set `bscfmax` @@ -114,4 +115,5 @@ Task meeting EU-JA #16, Fusion for Energy, Garching, 24--25 June 2014 [^7]: Menard et al. (2016), Nuclear Fusion, 56, 106023 [^8]: H.R. Wilson, Nuclear Fusion **32** (1992) 257 [^9]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **6** (1999) 2834 -[^10]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140 \ No newline at end of file +[^10]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140 +[^11]: Ryosuke Sakai, Takaaki Fujita, Atsushi Okamoto, Derivation of bootstrap current fraction scaling formula for 0-D system code analysis, Fusion Engineering and Design, Volume 149, 2019, 111322, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2019.111322. diff --git a/process/physics.py b/process/physics.py index ac895a08cc..e616abff60 100644 --- a/process/physics.py +++ b/process/physics.py @@ -990,6 +990,19 @@ def physics(self): current_drive_variables.cboot * self.bootstrap_fraction_sauter() ) + current_drive_variables.bscf_sakai = ( + current_drive_variables.cboot + * self.bootstrap_fraction_sakai( + betap=physics_variables.betap, + q95=physics_variables.q95, + q0=physics_variables.q0, + alphan=physics_variables.alphan, + alphat=physics_variables.alphat, + eps=physics_variables.eps, + rli=physics_variables.rli, + ) + ) + if current_drive_variables.bscfmax < 0.0e0: current_drive_variables.bootipf = abs(current_drive_variables.bscfmax) current_drive_variables.plasipf = current_drive_variables.bootipf @@ -1002,6 +1015,10 @@ def physics(self): current_drive_variables.bootipf = current_drive_variables.bscf_wilson elif physics_variables.ibss == 4: current_drive_variables.bootipf = current_drive_variables.bscf_sauter + elif physics_variables.ibss == 5: + # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current + # So the diamagnetic current calculation should be turned off when using, (idia = 0). + current_drive_variables.bootipf = current_drive_variables.bscf_sakai else: error_handling.idiags[0] = physics_variables.ibss error_handling.report_error(75) @@ -4064,6 +4081,14 @@ def outplas(self): current_drive_variables.bscf_wilson, "OP ", ) + + po.ovarrf( + self.outfile, + "Bootstrap fraction (Sakai)", + "(bscf_sakai)", + current_drive_variables.bscf_sakai, + "OP ", + ) po.ovarrf( self.outfile, "Diamagnetic fraction (Hender)", @@ -4115,6 +4140,11 @@ def outplas(self): self.outfile, " (Sauter et al bootstrap current fraction model used)", ) + elif physics_variables.ibss == 5: + po.ocmmnt( + self.outfile, + " (Sakai et al bootstrap current fraction model used)", + ) if physics_variables.idia == 0: po.ocmmnt( @@ -4662,6 +4692,55 @@ def bootstrap_fraction_sauter(): return np.sum(da * jboot, axis=0) / physics_variables.plascur + @staticmethod + def bootstrap_fraction_sakai( + betap: float, + q95: float, + q0: float, + alphan: float, + alphat: float, + eps: float, + rli: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Sakai formula. + + Parameters: + betap (float): Plasma poloidal beta. + q95 (float): Safety factor at 95% of the plasma radius. + q0 (float): Safety factor at the magnetic axis. + alphan (float): Density profile index + alphat (float): Temperature profile index + eps (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + The profile assumed for the alphan anf alpat indexes is only a prabolic profile without a pedestal (L-mode). + The Root Mean Squared Error for the fitting database of this formula was 0.025 + Concentrating on the positive shear plasmas using the ACCOME code equilibria with the fully non-inductively driven + conditions with neutral beam (NB) injection only are calculated. + The electron temperature and the ion temperature were assumed to be equal + This can be used for all apsect ratios. + The diamagnetic fraction is included in this formula. + + References: + Ryosuke Sakai, Takaaki Fujita, Atsushi Okamoto, Derivation of bootstrap current fraction scaling formula for 0-D system code analysis, + Fusion Engineering and Design, Volume 149, 2019, 111322, ISSN 0920-3796, + https://doi.org/10.1016/j.fusengdes.2019.111322. + """ + # Sakai states that the ACCOME dataset used has the toridal diamagnetic current included in the bootstrap current + # So the diamganetic current should not be calculated with this. idia = 0 + return ( + 10 ** (0.951 * eps - 0.948) + * betap ** (1.226 * eps + 1.584) + * rli ** (-0.184 - 0.282) + * (q95 / q0) ** (-0.042 * eps - 0.02) + * alphan ** (0.13 * eps + 0.05) + * alphat ** (0.502 * eps - 0.273) + ) + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/process/utilities/errorlist.json b/process/utilities/errorlist.json index dafabd320e..69611229a3 100644 --- a/process/utilities/errorlist.json +++ b/process/utilities/errorlist.json @@ -8,7 +8,7 @@ "comment2": [ "Increment n_errortypes if an error is added to this list" ], - "n_errortypes": 283, + "n_errortypes": 284, "errors": [ { "no": 1, @@ -1424,6 +1424,11 @@ "no": 283, "level": 3, "message": "[tfcoil]: Can only have i_tf_turns_integer = 1 with i_tf_wp_geom = 0" + }, + { + "no": 284, + "level": 3, + "message": "CHECK: idia = 0 should be used with the Sakai plasma current scaling." } ] } diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index 45b357e961..d8db4048a5 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -42,6 +42,9 @@ module current_drive_variables real(dp) :: bscf_wilson !! bootstrap current fraction, Wilson et al model + real(dp) :: bscf_sakai + !! Bootstrap current fraction, Sakai et al model + real(dp) :: cboot !! bootstrap current fraction multiplier (`ibss=1`) diff --git a/source/fortran/initial.f90 b/source/fortran/initial.f90 index 9a1512d740..97f0b0b45d 100755 --- a/source/fortran/initial.f90 +++ b/source/fortran/initial.f90 @@ -258,7 +258,7 @@ subroutine check use physics_variables, only: aspect, fdeut, fgwped, fhe3, & fgwsep, ftrit, ibss, i_single_null, icurr, idivrt, ishape, & iradloss, isc, ipedestal, ilhthresh, itart, nesep, rhopedn, rhopedt, & - rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar + rnbeam, neped, te, tauee_in, tesep, teped, itartpf, ftar, idia use pulse_variables, only: lpulse use reinke_variables, only: fzactual, impvardiv use tfcoil_variables, only: casthi, casthi_is_fraction, casths, i_tf_sup, & @@ -828,6 +828,10 @@ subroutine check call report_error(283) end if + if ( ibss == 5 .and. idia /= 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. diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 9654d2e262..f5dd77b16b 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -619,7 +619,7 @@ subroutine parse_input_file(in_file,out_file,show_changes) call parse_real_variable('taumax', taumax, 0.1D0, 100.0D0, & 'Maximum allowed energy confinement time (s)') case ('ibss') - call parse_int_variable('ibss', ibss, 1, 4, & + call parse_int_variable('ibss', ibss, 1, 5, & 'Switch for bootstrap scaling') case ('iculbl') call parse_int_variable('iculbl', iculbl, 0, 3, & diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index 995df178f9..7999fc5b7d 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -249,6 +249,7 @@ module physics_variables !! - =2 for Nevins et al general scaling !! - =3 for Wilson et al numerical scaling !! - =4 for Sauter et al scaling + !! - =5 for Sakai et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 125569741a..32bbac8ce6 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -468,6 +468,90 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs) +class BootstrapFractionSakaiParam(NamedTuple): + betap: Any = None + + q95: Any = None + + q0: Any = None + + alphan: Any = None + + alphat: Any = None + + eps: Any = None + + rli: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionsakaiparam", + ( + BootstrapFractionSakaiParam( + betap=1.3184383457774960, + q95=3.5151046634673557, + q0=1.0, + alphan=1.0, + alphat=1.45, + eps=1 / 3, + rli=1.2098126022585098, + expected_bfs=0.34204201506155418, + ), + BootstrapFractionSakaiParam( + betap=1.1701245502231756, + q95=5.1746754543339177, + q0=2.0, + alphan=0.9, + alphat=1.3999999999999999, + eps=1 / 1.8, + rli=0.3, + expected_bfs=0.90349498124262029, + ), + ), +) +def test_bootstrap_fraction_sakai(bootstrapfractionsakaiparam, monkeypatch, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_sakai. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + and tests/regression/input_files/st_regression.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + + :param monkeypatch: pytest fixture used to mock module/class variables + :type monkeypatch: _pytest.monkeypatch.monkeypatch + """ + + monkeypatch.setattr(physics_variables, "betap", bootstrapfractionsakaiparam.betap) + + monkeypatch.setattr(physics_variables, "q95", bootstrapfractionsakaiparam.q95) + + monkeypatch.setattr(physics_variables, "q0", bootstrapfractionsakaiparam.q0) + + monkeypatch.setattr(physics_variables, "alphan", bootstrapfractionsakaiparam.alphan) + + monkeypatch.setattr(physics_variables, "alphat", bootstrapfractionsakaiparam.alphat) + + monkeypatch.setattr(physics_variables, "eps", bootstrapfractionsakaiparam.eps) + + monkeypatch.setattr(physics_variables, "rli", bootstrapfractionsakaiparam.rli) + + bfs = physics.bootstrap_fraction_sakai( + betap=bootstrapfractionsakaiparam.betap, + q95=bootstrapfractionsakaiparam.q95, + q0=bootstrapfractionsakaiparam.q0, + alphan=bootstrapfractionsakaiparam.alphan, + alphat=bootstrapfractionsakaiparam.alphat, + eps=bootstrapfractionsakaiparam.eps, + rli=bootstrapfractionsakaiparam.rli, + ) + + assert bfs == pytest.approx(bootstrapfractionsakaiparam.expected_bfs) + + class CulcurParam(NamedTuple): normalised_total_beta: Any = None