From 3fdb328829ff933e86c9980292e4cef374d45559 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 22 Oct 2024 14:14:16 +0100 Subject: [PATCH 01/10] First implementation of ARIES bootstrap model without tests --- .../plasma_current/bootstrap_current.md | 29 +++++++-- process/physics.py | 62 ++++++++++++++++++- source/fortran/current_drive_variables.f90 | 3 + source/fortran/input.f90 | 2 +- source/fortran/physics_variables.f90 | 1 + 5 files changed, 91 insertions(+), 6 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index e3c44722..a88508fb 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -7,7 +7,7 @@ Some more info can be found [here](https://wiki.fusion.ciemat.es/wiki/Bootstrap_ The fraction of the plasma current provided by the bootstrap effect can be either input into the code directly, or calculated using one of five -methods, as summarised here. Note that methods `i_bootstrap_current = 1-3 & 5` do not take into account the +methods, as summarised here. Note that methods `i_bootstrap_current = 1-3 & 5-6` do not take into account the existence of pedestals, whereas the Sauter et al. scaling (`i_bootstrap_current = 4`) allows general profiles to be used. @@ -460,7 +460,7 @@ The $-1$ subscript in this case refers to the value of the variable in the previ ### Sakai Scaling | `bootstrap_fraction_sakai()` -Is selected by setting `i_bootstrap_current = 5`[^5] +Is selected by setting `i_bootstrap_current = 5`[^6] $$ f_{\text{BS}} = 10^{0.951 \epsilon - 0.948} \cdot \beta_p^{1.226 \epsilon + 1.584} \cdot l_i^{-0.184\epsilon - 0.282} \cdot \left(\frac{q_{95}}{q_0}\right)^{-0.042 \epsilon - 0.02} \\ @@ -469,6 +469,28 @@ $$ The model includes the toroidal diamagnetic current in the calculation due to the dataset, so `i_diamagnetic_current = 0` can only be used with it +------------------- + +### ARIES Scaling | `bootstrap_fraction_aries()` + +Is selected by setting `i_bootstrap_current = 6`[^8] + +$$ +a_1 = 1.10-1.165l_{\text{i}}+0.47l_{\text{i}}^2 +$$ + +$$ +b_1 = 0.806-0.885l_{\text{i}}+0.297l_{\text{i}}^2 +$$ + +$$ +C_{\text{BS}} = a_1+b_1\left(\frac{n(0)}{\langle n \rangle}\right) +$$ + +$$ +f_{\text{BS}} = C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}} +$$ + --------------------- ## Setting of maximum desirable bootstrap current fraction @@ -500,5 +522,4 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN [^5]: O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: “Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime” [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 [^6]: 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. [^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - +[^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. \ No newline at end of file diff --git a/process/physics.py b/process/physics.py index fb794f35..72a8e48f 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1718,6 +1718,17 @@ def physics(self): ) ) + current_drive_variables.bscf_aries = ( + current_drive_variables.cboot + * self.bootstrap_fraction_aries( + betap=physics_variables.betap, + rli=physics_variables.rli, + core_density=physics_variables.ne0, + average_desnity=physics_variables.dene, + inverse_aspect=physics_variables.eps, + ) + ) + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( current_drive_variables.bootstrap_current_fraction_max @@ -1748,6 +1759,10 @@ def physics(self): current_drive_variables.bootstrap_current_fraction = ( current_drive_variables.bscf_sakai ) + elif physics_variables.i_bootstrap_current == 6: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_sauter + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -4837,7 +4852,6 @@ def outplas(self): current_drive_variables.bscf_wilson, "OP ", ) - po.ovarrf( self.outfile, "Bootstrap fraction (Sakai)", @@ -4845,6 +4859,14 @@ def outplas(self): current_drive_variables.bscf_sakai, "OP ", ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (ARIES)", + "(bscf_aries)", + current_drive_variables.bscf_aries, + "OP ", + ) + po.ovarrf( self.outfile, "Diamagnetic fraction (Hender)", @@ -5536,6 +5558,7 @@ def bootstrap_fraction_sakai( alphan (float): Density profile index alphat (float): Temperature profile index eps (float): Inverse aspect ratio. + rli (float): Internal Inductance Returns: float: The calculated bootstrap fraction. @@ -5565,6 +5588,43 @@ def bootstrap_fraction_sakai( * alphat ** (0.502 * eps - 0.273) ) + @staticmethod + def bootstrap_fraction_aries( + betap: float, + rli: float, + core_density: float, + average_desnity: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the ARIES formula. + + Parameters: + betap (float): Plasma poloidal beta. + rli (float): Plasma normalized internal inductance. + core_density (float): Core plasma density. + average_density (float): Average plasma density. + inverse_aspect (float): Inverse aspect ratio.o. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + + References: + - Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants,” + Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, + doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. + + """ + # Using the standard variable naming from the ARIES paper + a_1 = 1.10-1.165*rli+0.47*rli**2 + b_1 = 0.806 - 0.885*rli + 0.297*rli**2 + + c_bs = a_1+b_1*(core_density/average_desnity) + + return c_bs * np.sqrt(inverse_aspect) * betap + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index d68e097d..7524ac1a 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -45,6 +45,9 @@ module current_drive_variables real(dp) :: bscf_sakai !! Bootstrap current fraction, Sakai et al model + real(dp) :: bscf_aries + !! Bootstrap current fraction, ARIES model + real(dp) :: cboot !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 7f884432..94d732bd 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,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 ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 5, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 6, & '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 fdc06965..ca496ea7 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -250,6 +250,7 @@ module physics_variables !! - =3 for Wilson et al numerical scaling !! - =4 for Sauter et al scaling !! - =5 for Sakai et al scaling + !! - =6 for ARIES scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) From 3ebb9872e85516691661c494c50bb5387a6e910b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 22 Oct 2024 15:27:17 +0100 Subject: [PATCH 02/10] Update plasma_profiles.md and physics.py - Add calculation of volume averaged pressure in plasma_profiles.md. This also add the core pressure and volume averaged pressure to the output. Assuming a parabolic profile - Implement bootstrap fraction calculation using Andrade et al formula in physics.py. This is currently implemented and exported with no tests and notes of applicability regimes --- .../profiles/plasma_profiles.md | 7 ++ process/physics.py | 79 ++++++++++++++++++- process/plasma_profiles.py | 4 + source/fortran/current_drive_variables.f90 | 3 + source/fortran/input.f90 | 2 +- source/fortran/physics_variables.f90 | 4 + 6 files changed, 94 insertions(+), 5 deletions(-) diff --git a/documentation/proc-pages/physics-models/profiles/plasma_profiles.md b/documentation/proc-pages/physics-models/profiles/plasma_profiles.md index 7c262472..e211f3d2 100644 --- a/documentation/proc-pages/physics-models/profiles/plasma_profiles.md +++ b/documentation/proc-pages/physics-models/profiles/plasma_profiles.md @@ -407,10 +407,17 @@ $$ With the coefficients used to turn the temperature from $\text{keV}$ back to Joules. Since we multiply the density and temperature profiles together to get the pressure we can find the pressure profile factor by adding the profile factors for temperature and density. + $$ \alpha_p = \alpha_n + \alpha_T $$ +The volume averaged pressure can then be set if we assume the pressure also has a parabolic profile. Using the standard relation used for both density and temeprature we can set the volume averaged pressure as: + +$$ +\langle p \rangle = \frac{p_0}{\alpha_p+1} +$$ + ???+ note "Pressure profile factor" The calculation of $\alpha_p$ is only valid assuming a parabolic profile case. The calculatio of $p_0$ is still true as the core values are calculated independantly for each profile type. diff --git a/process/physics.py b/process/physics.py index 72a8e48f..6fdbbdd4 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1729,6 +1729,16 @@ def physics(self): ) ) + current_drive_variables.bscf_andrade = ( + current_drive_variables.cboot + * self.bootstrap_fraction_andrade( + betap=physics_variables.betap, + core_pressure=physics_variables.p0, + average_pressure=physics_variables.vol_avg_pressure, + inverse_aspect=physics_variables.eps, + ) + ) + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( current_drive_variables.bootstrap_current_fraction_max @@ -1761,8 +1771,12 @@ def physics(self): ) elif physics_variables.i_bootstrap_current == 6: current_drive_variables.bootstrap_current_fraction = ( - current_drive_variables.bscf_sauter - ) + current_drive_variables.bscf_aries + ) + elif physics_variables.i_bootstrap_current == 7: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_andrade + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -3596,6 +3610,20 @@ def outplas(self): physics_variables.dnla, "OP ", ) + po.ovarre( + self.outfile, + "Plasma pressure on axis (Pa)", + "(p0)", + physics_variables.p0, + "OP ", + ) + po.ovarre( + self.outfile, + "Volume averaged plasma pressure (Pa)", + "(vol_avg_pressure)", + physics_variables.vol_avg_pressure, + "OP ", + ) if stellarator_variables.istell == 0: po.ovarre( @@ -4866,7 +4894,14 @@ def outplas(self): current_drive_variables.bscf_aries, "OP ", ) - + po.ovarrf( + self.outfile, + "Bootstrap fraction (Andrade)", + "(bscf_andrade)", + current_drive_variables.bscf_andrade, + "OP ", + ) + po.ovarrf( self.outfile, "Diamagnetic fraction (Hender)", @@ -5604,7 +5639,7 @@ def bootstrap_fraction_aries( rli (float): Plasma normalized internal inductance. core_density (float): Core plasma density. average_density (float): Average plasma density. - inverse_aspect (float): Inverse aspect ratio.o. + inverse_aspect (float): Inverse aspect ratio. Returns: float: The calculated bootstrap fraction. @@ -5625,6 +5660,42 @@ def bootstrap_fraction_aries( return c_bs * np.sqrt(inverse_aspect) * betap + @staticmethod + def bootstrap_fraction_andrade( + betap: float, + core_pressure: float, + average_pressure: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Andrade et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + core_pressure (float): Core plasma pressure. + average_pressure (float): Average plasma pressure. + inverse_aspect (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + + References: + - M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” + Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, + doi: https://doi.org/10.1088/0741-3335/50/6/065001. + + """ + + # Using the standard variable naming from the Andrade et.al. paper + c_p = core_pressure / average_pressure + + # Error +- 0.0007 + c_bs = 0.2340 + + return c_bs * np.sqrt(inverse_aspect) * betap * c_p**0.8 + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index ea804c24..787afec6 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -264,6 +264,10 @@ def calculate_profile_factors(self) -> None: physics_variables.alphap = physics_variables.alphan + physics_variables.alphat + # Shall assume that the pressure profile is parabolic. Can this find volume average from + # profile index and core value the same as for density and temperature + physics_variables.vol_avg_pressure = physics_variables.p0 / (physics_variables.alphap+1) + @staticmethod def calculate_parabolic_profile_factors() -> None: """ diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index 7524ac1a..89609139 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -48,6 +48,9 @@ module current_drive_variables real(dp) :: bscf_aries !! Bootstrap current fraction, ARIES model + real(dp) :: bscf_andrade + !! Bootstrap current fraction, Andrade et al model + real(dp) :: cboot !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 94d732bd..99c6d485 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,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 ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 6, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 7, & '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 ca496ea7..7886853f 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -251,6 +251,7 @@ module physics_variables !! - =4 for Sauter et al scaling !! - =5 for Sakai et al scaling !! - =6 for ARIES scaling + !! - =7 for Andrade et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) @@ -551,6 +552,9 @@ module physics_variables real(dp) :: p0 !! central total plasma pressure (Pa) + real(dp) :: vol_avg_pressure + !! Volume average plasma pressure (Pa) + real(dp) :: palppv !! alpha power per volume (MW/m3) From e79d052bf9b9ffcd8e93a2d5ff08643fd732ea25 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 22 Oct 2024 15:27:17 +0100 Subject: [PATCH 03/10] Update plasma_profiles.md and physics.py - Add calculation of volume averaged pressure in plasma_profiles.md. This also add the core pressure and volume averaged pressure to the output. Assuming a parabolic profile - Implement bootstrap fraction calculation using Andrade et al formula in physics.py. This is currently implemented and exported with no tests and notes of applicability regimes --- .../plasma_current/bootstrap_current.md | 21 ++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index a88508fb..ea17c015 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -493,6 +493,24 @@ $$ --------------------- +### Andrade Scaling | `bootstrap_fraction_andrade()` + +Is selected by setting `i_bootstrap_current = 7`[^9] + +$$ +C_{\text{BS}} = 0.2340 \pm 0.0007 +$$ + +$$ +c_p = \frac{p_0}{\langle p \rangle} +$$ + +$$ +f_{\text{BS}} = C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}}c_p^{0.8} +$$ + +--------------------- + ## Setting of maximum desirable bootstrap current fraction The variable `bootstrap_current_fraction_max` can be set to the value of maximum desirable bootstrap current fraction for a specific design. When optimising if the value of the calculated `bootstrap_current_fraction` for the model selected with `i_bootstrap_current` exceeds this value, then `bootstrap_current_fraction` is set to the value of `bootstrap_current_fraction_max`. @@ -522,4 +540,5 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN [^5]: O. Sauter, C. Angioni, Y. R. Lin-Liu; Erratum: “Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime” [Phys. Plasmas 6, 2834 (1999)]. Phys. Plasmas 1 December 2002; 9 (12): 5140. https://doi.org/10.1063/1.1517052 [^6]: 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. [^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 -[^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. \ No newline at end of file +[^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. +[^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001. \ No newline at end of file From 374345c3049bbfc28b9a506e11e1d9bc22e7d807 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 23 Oct 2024 12:01:18 +0100 Subject: [PATCH 04/10] Add new bootstrap scaling distribution plot to plot_proc.py --- process/io/plot_proc.py | 103 +++++++++++++++++++++++++++++++++------- 1 file changed, 86 insertions(+), 17 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index d04a6d85..192f85a3 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -1876,14 +1876,14 @@ def plot_tf_wp(axis, mfile_data, scan: int) -> None: ) ) - plt.minorticks_on() - plt.xlim(0.0, r_tf_inboard_out * 1.1) - plt.ylim((y14[-1] * 1.25), (-y14[-1] * 1.25)) + axis.minorticks_on() + axis.set_xlim(0.0, r_tf_inboard_out * 1.1) + axis.set_ylim((y14[-1] * 1.25), (-y14[-1] * 1.25)) - plt.title("Top-down view of inboard TF coil at midplane") - plt.xlabel("Radial distance [m]") - plt.ylabel("Toroidal distance [m]") - plt.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left") + axis.set_title("Top-down view of inboard TF coil at midplane") + axis.set_xlabel("Radial distance [m]") + axis.set_ylabel("Toroidal distance [m]") + axis.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left") def plot_tf_turn(axis, mfile_data, scan: int) -> None: @@ -1974,8 +1974,8 @@ def plot_tf_turn(axis, mfile_data, scan: int) -> None: edgecolor="black", ), ) - plt.xlim(-turn_width * 0.05, turn_width * 1.05) - plt.ylim(-turn_width * 0.05, turn_width * 1.05) + axis.set_xlim(-turn_width * 0.05, turn_width * 1.05) + axis.set_ylim(-turn_width * 0.05, turn_width * 1.05) # Non square turns elif integer_turns == 1: @@ -2026,14 +2026,14 @@ def plot_tf_turn(axis, mfile_data, scan: int) -> None: ), ) - plt.xlim(-turn_width * 0.05, turn_width * 1.05) - plt.ylim(-turn_height * 0.05, turn_height * 1.05) + axis.set_xlim(-turn_width * 0.05, turn_width * 1.05) + axis.set_ylim(-turn_height * 0.05, turn_height * 1.05) - plt.minorticks_on() - plt.title("WP Turn Structure") - plt.xlabel("X [mm]") - plt.ylabel("Y [mm]") - plt.legend(bbox_to_anchor=(0.0, -0.25), loc="upper left") + axis.minorticks_on() + axis.set_title("WP Turn Structure") + axis.set_xlabel("X [mm]") + axis.set_ylabel("Y [mm]") + axis.legend(loc="upper right", bbox_to_anchor=(1.0, -0.25)) def plot_pf_coils(axis, mfile_data, scan, colour_scheme): @@ -2898,10 +2898,72 @@ def plot_current_drive_info(axis, mfile_data, scan): plot_info(axis, data, mfile_data, scan) +def plot_bootstrap_comparison(axis, mfile_data, scan): + """Function to plot a scatter box plot of bootstrap current fractions. + + Arguments: + axis --> axis object to plot to + mfile_data --> MFILE data object + scan --> scan number to use + """ + + boot_ipdg = mfile_data.data["bscf_iter89"].get_scan(scan) + boot_sauter = mfile_data.data["bscf_sauter"].get_scan(scan) + boot_nenins = mfile_data.data["bscf_nevins"].get_scan(scan) + boot_wilson = mfile_data.data["bscf_wilson"].get_scan(scan) + boot_sakai = mfile_data.data["bscf_sakai"].get_scan(scan) + boot_aries = mfile_data.data["bscf_aries"].get_scan(scan) + boot_andrade = mfile_data.data["bscf_andrade"].get_scan(scan) + + # Data for the box plot + data = [ + boot_ipdg, + boot_sauter, + boot_nenins, + boot_wilson, + boot_sakai, + boot_aries, + boot_andrade, + ] + labels = [ + "IPDG", + "Sauter", + "Nevins", + "Wilson", + "Sakai", + "ARIES", + "Andrade", + ] + + x = np.ones(len(data)) + + # Labels for the box plot + plt.scatter(x, data, color="black") + # Create the box plot + axis.boxplot(data) + # Calculate average and standard deviation + avg_bootstrap = np.mean(data) + std_bootstrap = np.std(data) + + # Plot average and standard deviation as text + axis.text( + 1.1, 0.9, f"Average: {avg_bootstrap:.2f}", transform=axis.transAxes, fontsize=10 + ) + axis.text( + 1.1, 0.8, f"Std Dev: {std_bootstrap:.2f}", transform=axis.transAxes, fontsize=10 + ) + axis.set_xticks([]) + for i, value in enumerate(data): + axis.text(x[i] + 0.1, value, labels[i], fontsize=8, verticalalignment="center") + axis.set_title("Bootstrap Current Fraction Comparison") + axis.set_ylabel("Bootstrap Current Fraction") + + def main_plot( fig1, fig2, fig3, + fig4, m_file_data, scan, imp="../data/lz_non_corona_14_elements/", @@ -2988,7 +3050,7 @@ def main_plot( plot_current_drive_info(plot_6, m_file_data, scan) fig1.subplots_adjust(wspace=0.25, hspace=0.25) - # Can only plot WP and turn sturcutre if superconducting coil at the moment + # Can only plot WP and turn structure if superconducting coil at the moment if m_file_data.data["i_tf_sup"].get_scan(scan) == 1: # TF coil with WP plot_7 = fig3.add_subplot(321) @@ -2998,6 +3060,9 @@ def main_plot( plot_8 = fig3.add_subplot(322, aspect="equal") plot_tf_turn(plot_8, m_file_data, scan) + plot_9 = fig4.add_subplot(221) + plot_bootstrap_comparison(plot_9, m_file_data, scan) + def main(args=None): # TODO The use of globals here isn't ideal, but is required to get main() @@ -3252,12 +3317,14 @@ def main(args=None): page1 = plt.figure(figsize=(12, 9), dpi=80) page2 = plt.figure(figsize=(12, 9), dpi=80) page3 = plt.figure(figsize=(12, 9), dpi=80) + page4 = plt.figure(figsize=(12, 9), dpi=80) # run main_plot main_plot( page1, page2, page3, + page4, m_file, scan=scan, demo_ranges=demo_ranges, @@ -3269,6 +3336,7 @@ def main(args=None): pdf.savefig(page1) pdf.savefig(page2) pdf.savefig(page3) + pdf.savefig(page4) # show fig if option used if args.show: @@ -3277,6 +3345,7 @@ def main(args=None): plt.close(page1) plt.close(page2) plt.close(page3) + plt.close(page4) if __name__ == "__main__": From 19faafed47c9a09e3641e02e0c5442e78942e0c2 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 23 Oct 2024 15:27:37 +0100 Subject: [PATCH 05/10] Add basics of Hoang bootstrap model and add to plot_proc output, no test added --- .../plasma_current/bootstrap_current.md | 17 ++++- process/io/plot_proc.py | 3 + process/physics.py | 76 ++++++++++++++++++- process/plasma_profiles.py | 4 +- source/fortran/current_drive_variables.f90 | 3 + source/fortran/input.f90 | 2 +- source/fortran/physics_variables.f90 | 1 + 7 files changed, 100 insertions(+), 6 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index ea17c015..ed61cd1e 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -511,6 +511,20 @@ $$ --------------------- +### Hoang Scaling | `bootstrap_fraction_hoang()` + +Is selected by setting `i_bootstrap_current = 8`[^10] + +$$ +C_{\text{BS}} = \sqrt{\frac{\alpha_{\text{p}}}{\alpha_{\text{j}}}} +$$ + +$$ +f_{\text{BS}} = 0.4C_{\text{BS}} \sqrt{\epsilon}\beta_{\text{p}}^{0.9} +$$ + +--------------------- + ## Setting of maximum desirable bootstrap current fraction The variable `bootstrap_current_fraction_max` can be set to the value of maximum desirable bootstrap current fraction for a specific design. When optimising if the value of the calculated `bootstrap_current_fraction` for the model selected with `i_bootstrap_current` exceeds this value, then `bootstrap_current_fraction` is set to the value of `bootstrap_current_fraction_max`. @@ -541,4 +555,5 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN [^6]: 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. [^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 [^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. -[^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001. \ No newline at end of file +[^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001. +[^10]: G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, Jan. 1997, doi: https://doi.org/10.1063/1.53414. \ No newline at end of file diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 192f85a3..373b0294 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2914,6 +2914,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_sakai = mfile_data.data["bscf_sakai"].get_scan(scan) boot_aries = mfile_data.data["bscf_aries"].get_scan(scan) boot_andrade = mfile_data.data["bscf_andrade"].get_scan(scan) + boot_hoang = mfile_data.data["bscf_hoang"].get_scan(scan) # Data for the box plot data = [ @@ -2924,6 +2925,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_sakai, boot_aries, boot_andrade, + boot_hoang, ] labels = [ "IPDG", @@ -2933,6 +2935,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): "Sakai", "ARIES", "Andrade", + "Hoang", ] x = np.ones(len(data)) diff --git a/process/physics.py b/process/physics.py index 6fdbbdd4..35287b44 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1739,6 +1739,16 @@ def physics(self): ) ) + current_drive_variables.bscf_hoang = ( + current_drive_variables.cboot + * self.bootstrap_fraction_hoang( + betap=physics_variables.betap, + pressure_index=physics_variables.alphap, + current_index=physics_variables.alphaj, + inverse_aspect=physics_variables.eps, + ) + ) + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( current_drive_variables.bootstrap_current_fraction_max @@ -1777,6 +1787,10 @@ def physics(self): current_drive_variables.bootstrap_current_fraction = ( current_drive_variables.bscf_andrade ) + elif physics_variables.i_bootstrap_current == 8: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_hoang + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -4901,6 +4915,13 @@ def outplas(self): current_drive_variables.bscf_andrade, "OP ", ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Hoang)", + "(bscf_hoang)", + current_drive_variables.bscf_hoang, + "OP ", + ) po.ovarrf( self.outfile, @@ -4958,6 +4979,21 @@ def outplas(self): self.outfile, " (Sakai et al bootstrap current fraction model used)", ) + elif physics_variables.i_bootstrap_current == 6: + po.ocmmnt( + self.outfile, + " (ARIES bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 7: + po.ocmmnt( + self.outfile, + " (Andrade et al bootstrap current fraction model used)", + ) + elif physics_variables.i_bootstrap_current == 8: + po.ocmmnt( + self.outfile, + " (Hoang et al bootstrap current fraction model used)", + ) if physics_variables.i_diamagnetic_current == 0: po.ocmmnt( @@ -5653,10 +5689,10 @@ def bootstrap_fraction_aries( """ # Using the standard variable naming from the ARIES paper - a_1 = 1.10-1.165*rli+0.47*rli**2 - b_1 = 0.806 - 0.885*rli + 0.297*rli**2 + a_1 = 1.10 - 1.165 * rli + 0.47 * rli**2 + b_1 = 0.806 - 0.885 * rli + 0.297 * rli**2 - c_bs = a_1+b_1*(core_density/average_desnity) + c_bs = a_1 + b_1 * (core_density / average_desnity) return c_bs * np.sqrt(inverse_aspect) * betap @@ -5696,6 +5732,40 @@ def bootstrap_fraction_andrade( return c_bs * np.sqrt(inverse_aspect) * betap * c_p**0.8 + @staticmethod + def bootstrap_fraction_hoang( + betap: float, + pressure_index: float, + current_index: float, + inverse_aspect: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Hoang et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + current_index (float): Current profile index. + inverse_aspect (float): Inverse aspect ratio. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + + References: + - G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, + Jan. 1997, doi: https://doi.org/10.1063/1.53414. + ‌ + """ + + # Using the standard variable naming from the Hoang et.al. paper + # These terms do not equal the profile indexes, though are closely linked + + c_bs = np.sqrt(pressure_index / current_index) + + return 0.4 * np.sqrt(inverse_aspect) * betap**0.9 * c_bs + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index 787afec6..6f5c3d2f 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -266,7 +266,9 @@ def calculate_profile_factors(self) -> None: # Shall assume that the pressure profile is parabolic. Can this find volume average from # profile index and core value the same as for density and temperature - physics_variables.vol_avg_pressure = physics_variables.p0 / (physics_variables.alphap+1) + physics_variables.vol_avg_pressure = physics_variables.p0 / ( + physics_variables.alphap + 1 + ) @staticmethod def calculate_parabolic_profile_factors() -> None: diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index 89609139..b96ad80e 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -51,6 +51,9 @@ module current_drive_variables real(dp) :: bscf_andrade !! Bootstrap current fraction, Andrade et al model + real(dp) :: bscf_hoang + !! Bootstrap current fraction, Hoang et al model + real(dp) :: cboot !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 99c6d485..dbdb1e9f 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,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 ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 7, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 8, & '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 7886853f..efaafd30 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -252,6 +252,7 @@ module physics_variables !! - =5 for Sakai et al scaling !! - =6 for ARIES scaling !! - =7 for Andrade et al scaling + !! - =8 for Hoang et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) From 29e5ed8965415bb1cb12e0ff66249346f4c8bf44 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 23 Oct 2024 17:32:50 +0100 Subject: [PATCH 06/10] Implement the Wong bootstrap model. No test added, Details added to docs --- .../plasma_current/bootstrap_current.md | 29 +++++- process/physics.py | 89 ++++++++++++++++--- source/fortran/current_drive_variables.f90 | 3 + source/fortran/input.f90 | 2 +- source/fortran/physics_variables.f90 | 1 + 5 files changed, 109 insertions(+), 15 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index ed61cd1e..a9f154cd 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -525,6 +525,32 @@ $$ --------------------- +### Wong Scaling | `bootstrap_fraction_wong()` + +Is selected by setting `i_bootstrap_current = 9`[^11] + +$$ +C_{\text{BS}} = 0.773+0.019\kappa +$$ + +$$ +f_{\text{peak}} = \left(\int_0^1 \left(1-\rho^2 \right)^{\alpha_{\text{T}}} \left(1-\rho^2 \right)^{\alpha_{\text{n}}} \ \ \mathrm{d\rho}\right)^{-1} +$$ + +The integral above is set by the definite solution below + +$$ +\frac{\operatorname{B}\left(\frac{1}{2},{\alpha}_{n} + {\alpha}_{T} + 1\right)}{2} +$$ + +Assuming that $\alpha_{\text{n}} + \alpha_{\text{T}} > 0, \alpha_{\text{n}} + \alpha_{\text{T}} + 1 > 0$ + +$$ +f_{\text{BS}} = C_{\text{BS}} f_{\text{peak}}^{0.25}\beta_{\text{p}}\sqrt{\epsilon} +$$ + +--------------------- + ## Setting of maximum desirable bootstrap current fraction The variable `bootstrap_current_fraction_max` can be set to the value of maximum desirable bootstrap current fraction for a specific design. When optimising if the value of the calculated `bootstrap_current_fraction` for the model selected with `i_bootstrap_current` exceeds this value, then `bootstrap_current_fraction` is set to the value of `bootstrap_current_fraction_max`. @@ -556,4 +582,5 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN [^7]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 [^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. [^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001. -[^10]: G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, Jan. 1997, doi: https://doi.org/10.1063/1.53414. \ No newline at end of file +[^10]: G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, Jan. 1997, doi: https://doi.org/10.1063/1.53414. +[^11]: C.-P. Wong, J. C. Wesley, R. D. Stambaugh, and E. T. Cheng, “Toroidal reactor designs as a function of aspect ratio and elongation,” vol. 42, no. 5, pp. 547–556, May 2002, doi: https://doi.org/10.1088/0029-5515/42/5/307. \ No newline at end of file diff --git a/process/physics.py b/process/physics.py index 35287b44..8f93275d 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1,6 +1,7 @@ import math from typing import Tuple import numpy as np +import scipy import numba as nb from scipy.optimize import root_scalar from process.utilities.f2py_string_patch import f2py_compatible_to_string @@ -1749,6 +1750,17 @@ def physics(self): ) ) + current_drive_variables.bscf_wong = ( + current_drive_variables.cboot + * self.bootstrap_fraction_wong( + betap=physics_variables.betap, + density_index=physics_variables.alphan, + temperature_index=physics_variables.alphat, + inverse_aspect=physics_variables.eps, + elongation=physics_variables.kappa, + ) + ) + if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( current_drive_variables.bootstrap_current_fraction_max @@ -1791,6 +1803,10 @@ def physics(self): current_drive_variables.bootstrap_current_fraction = ( current_drive_variables.bscf_hoang ) + elif physics_variables.i_bootstrap_current == 9: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_wong + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -4922,6 +4938,14 @@ def outplas(self): current_drive_variables.bscf_hoang, "OP ", ) + + po.ovarrf( + self.outfile, + "Bootstrap fraction (Wong)", + "(bscf_wong)", + current_drive_variables.bscf_wong, + "OP ", + ) po.ovarrf( self.outfile, @@ -4994,6 +5018,12 @@ def outplas(self): self.outfile, " (Hoang et al bootstrap current fraction model used)", ) + elif physics_variables.i_bootstrap_current == 9: + po.ocmmnt( + self.outfile, + " (Wong et al bootstrap current fraction model used)", + ) + if physics_variables.i_diamagnetic_current == 0: po.ocmmnt( @@ -5740,23 +5770,22 @@ def bootstrap_fraction_hoang( inverse_aspect: float, ) -> float: """ - Calculate the bootstrap fraction using the Hoang et al formula. + Calculate the bootstrap fraction using the Hoang et al formula. - Parameters: - betap (float): Plasma poloidal beta. - pressure_index (float): Pressure profile index. - current_index (float): Current profile index. - inverse_aspect (float): Inverse aspect ratio. + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + current_index (float): Current profile index. + inverse_aspect (float): Inverse aspect ratio. - Returns: - float: The calculated bootstrap fraction. + Returns: + float: The calculated bootstrap fraction. - Notes: + Notes: - References: - - G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, - Jan. 1997, doi: https://doi.org/10.1063/1.53414. - ‌ + References: + - G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, + Jan. 1997, doi: https://doi.org/10.1063/1.53414. """ # Using the standard variable naming from the Hoang et.al. paper @@ -5766,6 +5795,40 @@ def bootstrap_fraction_hoang( return 0.4 * np.sqrt(inverse_aspect) * betap**0.9 * c_bs + @staticmethod + def bootstrap_fraction_wong( + betap: float, + density_index: float, + temperature_index: float, + inverse_aspect: float, + elongation: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Wong et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + density_index (float): Density profile index. + temperature_index (float): Temperature profile index. + inverse_aspect (float): Inverse aspect ratio. + elongation (float): Plasma elongation. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + + References: + - C.-P. Wong, J. C. Wesley, R. D. Stambaugh, and E. T. Cheng, “Toroidal reactor designs as a function of aspect ratio and elongation,” + vol. 42, no. 5, pp. 547–556, May 2002, doi: https://doi.org/10.1088/0029-5515/42/5/307. + """ + # Using the standard variable naming from the Wong et.al. paper + f_peak = 2.0/scipy.special.beta(0.5, density_index + temperature_index + 1) + + c_bs = 0.773 + 0.019 * elongation + + return c_bs * f_peak**0.25 * betap * np.sqrt(inverse_aspect) + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index b96ad80e..b0cd7f1d 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -54,6 +54,9 @@ module current_drive_variables real(dp) :: bscf_hoang !! Bootstrap current fraction, Hoang et al model + real(dp) :: bscf_wong + !! Bootstrap current fraction, Wong et al model + real(dp) :: cboot !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index dbdb1e9f..a629224e 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,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 ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 8, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 9, & '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 efaafd30..c609fb41 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -253,6 +253,7 @@ module physics_variables !! - =6 for ARIES scaling !! - =7 for Andrade et al scaling !! - =8 for Hoang et al scaling + !! - =9 for Wong et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) From 3f04b9110c5bd1a426937d8bec058ac0d784cb6e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 23 Oct 2024 17:33:16 +0100 Subject: [PATCH 07/10] Update bootstrap comparison model with new violin plot --- process/io/plot_proc.py | 34 ++++++++++++++++++++++++---------- 1 file changed, 24 insertions(+), 10 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 373b0294..606200a0 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2915,6 +2915,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_aries = mfile_data.data["bscf_aries"].get_scan(scan) boot_andrade = mfile_data.data["bscf_andrade"].get_scan(scan) boot_hoang = mfile_data.data["bscf_hoang"].get_scan(scan) + boot_wong = mfile_data.data["bscf_wong"].get_scan(scan) # Data for the box plot data = [ @@ -2926,6 +2927,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_aries, boot_andrade, boot_hoang, + boot_wong, ] labels = [ "IPDG", @@ -2936,30 +2938,42 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): "ARIES", "Andrade", "Hoang", + "Wong", ] x = np.ones(len(data)) - # Labels for the box plot - plt.scatter(x, data, color="black") + # Create the violin plot + axis.violinplot(data, showextrema=False) + # Create the box plot - axis.boxplot(data) - # Calculate average and standard deviation + axis.boxplot(data, showfliers=True) + + # Scatter plot for each data point + colors = plt.cm.plasma(np.linspace(0, 1, len(data))) + for i, value in enumerate(data): + axis.scatter(x[i], value, color=colors[i], label=labels[i], alpha=1.0) + axis.legend(loc='upper left', bbox_to_anchor=(1, 1)) + + # Calculate average, standard deviation, and median avg_bootstrap = np.mean(data) std_bootstrap = np.std(data) + median_bootstrap = np.median(data) - # Plot average and standard deviation as text + # Plot average, standard deviation, and median as text axis.text( - 1.1, 0.9, f"Average: {avg_bootstrap:.2f}", transform=axis.transAxes, fontsize=10 + 0.6, 0.9, f"Average: {avg_bootstrap:.4f}", transform=axis.transAxes, fontsize=9 ) axis.text( - 1.1, 0.8, f"Std Dev: {std_bootstrap:.2f}", transform=axis.transAxes, fontsize=10 + 0.6, 0.85, f"Standard Dev: {std_bootstrap:.4f}", transform=axis.transAxes, fontsize=9 ) - axis.set_xticks([]) - for i, value in enumerate(data): - axis.text(x[i] + 0.1, value, labels[i], fontsize=8, verticalalignment="center") + axis.text(0.6, 0.8, f"Median: {median_bootstrap:.4f}", transform=axis.transAxes, fontsize=9) + axis.set_title("Bootstrap Current Fraction Comparison") axis.set_ylabel("Bootstrap Current Fraction") + axis.set_xlim([0.5, 1.5]) + axis.set_xticks([]) + axis.set_xticklabels([]) def main_plot( From caebdda8d87a9cc85783d67b1b01f175c9b40703 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 23 Oct 2024 17:33:16 +0100 Subject: [PATCH 08/10] Update bootstrap comparison model with new violin plot --- process/io/plot_proc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 606200a0..02c414b2 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2947,7 +2947,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): axis.violinplot(data, showextrema=False) # Create the box plot - axis.boxplot(data, showfliers=True) + axis.boxplot(data, showfliers=True, showmeans=True, meanline=True, widths=0.3) # Scatter plot for each data point colors = plt.cm.plasma(np.linspace(0, 1, len(data))) From b434c08b5335090bcc8a5568a4fac1657536f0b6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 24 Oct 2024 09:35:10 +0100 Subject: [PATCH 09/10] Add new Gi bootstrap scaling formula and add to docs, no tests --- .../plasma_current/bootstrap_current.md | 17 ++++- process/io/plot_proc.py | 19 ++++- process/physics.py | 75 +++++++++++++++++-- source/fortran/current_drive_variables.f90 | 5 +- source/fortran/input.f90 | 2 +- source/fortran/physics_variables.f90 | 1 + 6 files changed, 107 insertions(+), 12 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index a9f154cd..7c36f3e7 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -551,6 +551,20 @@ $$ --------------------- +### Gi Scaling | `bootstrap_fraction_gi()` + +Is selected by setting `i_bootstrap_current = 10`[^12] + +$$ +C_{\text{BS}} = 0.474 \epsilon^{-0.1} \alpha_{\text{p}}^{0.974} \alpha_{\text{T}}^{-0.416} Z_{\text{eff}}^{0.178} \left(\frac{q_{95}}{q_0}\right)^{-0.133} +$$ + +$$ +f_{\text{BS}} = C_{\text{BS}} \beta_{\text{p}}\sqrt{\epsilon} +$$ + +--------------------- + ## Setting of maximum desirable bootstrap current fraction The variable `bootstrap_current_fraction_max` can be set to the value of maximum desirable bootstrap current fraction for a specific design. When optimising if the value of the calculated `bootstrap_current_fraction` for the model selected with `i_bootstrap_current` exceeds this value, then `bootstrap_current_fraction` is set to the value of `bootstrap_current_fraction_max`. @@ -583,4 +597,5 @@ Fusion Engineering and Design, Volume 89, Issue 11, 2014, Pages 2709-2715, ISSN [^8]: Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants, ”Fusion Engineering and Design, vol. 85, no. 2, pp. 243–265, Apr. 2010, doi: https://doi.org/10.1016/j.fusengdes.2010.02.015. [^9]: M. C. R. Andrade and G. O. Ludwig, “Scaling of bootstrap current on equilibrium and plasma profile parameters in tokamak plasmas,” Plasma Physics and Controlled Fusion, vol. 50, no. 6, pp. 065001–065001, Apr. 2008, doi: https://doi.org/10.1088/0741-3335/50/6/065001. [^10]: G. T. Hoang and R. V. Budny, “The bootstrap fraction in TFTR,” AIP conference proceedings, Jan. 1997, doi: https://doi.org/10.1063/1.53414. -[^11]: C.-P. Wong, J. C. Wesley, R. D. Stambaugh, and E. T. Cheng, “Toroidal reactor designs as a function of aspect ratio and elongation,” vol. 42, no. 5, pp. 547–556, May 2002, doi: https://doi.org/10.1088/0029-5515/42/5/307. \ No newline at end of file +[^11]: C.-P. Wong, J. C. Wesley, R. D. Stambaugh, and E. T. Cheng, “Toroidal reactor designs as a function of aspect ratio and elongation,” vol. 42, no. 5, pp. 547–556, May 2002, doi: https://doi.org/10.1088/0029-5515/42/5/307. +[^12]: K. Gi, M. Nakamura, Kenji Tobita, and Y. Ono, “Bootstrap current fraction scaling for a tokamak reactor design study,” Fusion Engineering and Design, vol. 89, no. 11, pp. 2709–2715, Aug. 2014, doi: https://doi.org/10.1016/j.fusengdes.2014.07.009. \ No newline at end of file diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 02c414b2..5e544013 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2916,6 +2916,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_andrade = mfile_data.data["bscf_andrade"].get_scan(scan) boot_hoang = mfile_data.data["bscf_hoang"].get_scan(scan) boot_wong = mfile_data.data["bscf_wong"].get_scan(scan) + boot_gi = mfile_data.data["bscf_gi"].get_scan(scan) # Data for the box plot data = [ @@ -2928,6 +2929,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): boot_andrade, boot_hoang, boot_wong, + boot_gi, ] labels = [ "IPDG", @@ -2939,6 +2941,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): "Andrade", "Hoang", "Wong", + "Gi", ] x = np.ones(len(data)) @@ -2953,7 +2956,7 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): colors = plt.cm.plasma(np.linspace(0, 1, len(data))) for i, value in enumerate(data): axis.scatter(x[i], value, color=colors[i], label=labels[i], alpha=1.0) - axis.legend(loc='upper left', bbox_to_anchor=(1, 1)) + axis.legend(loc="upper left", bbox_to_anchor=(1, 1)) # Calculate average, standard deviation, and median avg_bootstrap = np.mean(data) @@ -2965,9 +2968,19 @@ def plot_bootstrap_comparison(axis, mfile_data, scan): 0.6, 0.9, f"Average: {avg_bootstrap:.4f}", transform=axis.transAxes, fontsize=9 ) axis.text( - 0.6, 0.85, f"Standard Dev: {std_bootstrap:.4f}", transform=axis.transAxes, fontsize=9 + 0.6, + 0.85, + f"Standard Dev: {std_bootstrap:.4f}", + transform=axis.transAxes, + fontsize=9, + ) + axis.text( + 0.6, + 0.8, + f"Median: {median_bootstrap:.4f}", + transform=axis.transAxes, + fontsize=9, ) - axis.text(0.6, 0.8, f"Median: {median_bootstrap:.4f}", transform=axis.transAxes, fontsize=9) axis.set_title("Bootstrap Current Fraction Comparison") axis.set_ylabel("Bootstrap Current Fraction") diff --git a/process/physics.py b/process/physics.py index 8f93275d..41088f88 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1739,7 +1739,6 @@ def physics(self): inverse_aspect=physics_variables.eps, ) ) - current_drive_variables.bscf_hoang = ( current_drive_variables.cboot * self.bootstrap_fraction_hoang( @@ -1749,7 +1748,6 @@ def physics(self): inverse_aspect=physics_variables.eps, ) ) - current_drive_variables.bscf_wong = ( current_drive_variables.cboot * self.bootstrap_fraction_wong( @@ -1760,6 +1758,18 @@ def physics(self): elongation=physics_variables.kappa, ) ) + current_drive_variables.bscf_gi = ( + current_drive_variables.cboot + * self.bootstrap_fraction_gi( + betap=physics_variables.betap, + pressure_index=physics_variables.alphap, + temperature_index=physics_variables.alphat, + inverse_aspect=physics_variables.eps, + effective_charge=physics_variables.zeff, + q95=physics_variables.q95, + q0=physics_variables.q0, + ) + ) if current_drive_variables.bootstrap_current_fraction_max < 0.0e0: current_drive_variables.bootstrap_current_fraction = abs( @@ -1806,7 +1816,11 @@ def physics(self): elif physics_variables.i_bootstrap_current == 9: current_drive_variables.bootstrap_current_fraction = ( current_drive_variables.bscf_wong - ) + ) + elif physics_variables.i_bootstrap_current == 10: + current_drive_variables.bootstrap_current_fraction = ( + current_drive_variables.bscf_gi + ) else: error_handling.idiags[0] = physics_variables.i_bootstrap_current error_handling.report_error(75) @@ -4938,7 +4952,6 @@ def outplas(self): current_drive_variables.bscf_hoang, "OP ", ) - po.ovarrf( self.outfile, "Bootstrap fraction (Wong)", @@ -4946,6 +4959,13 @@ def outplas(self): current_drive_variables.bscf_wong, "OP ", ) + po.ovarrf( + self.outfile, + "Bootstrap fraction (Gi)", + "(bscf_gi)", + current_drive_variables.bscf_gi, + "OP ", + ) po.ovarrf( self.outfile, @@ -5022,8 +5042,12 @@ def outplas(self): po.ocmmnt( self.outfile, " (Wong et al bootstrap current fraction model used)", - ) - + ) + elif physics_variables.i_bootstrap_current == 10: + po.ocmmnt( + self.outfile, + " (Gi et al bootstrap current fraction model used)", + ) if physics_variables.i_diamagnetic_current == 0: po.ocmmnt( @@ -5829,6 +5853,45 @@ def bootstrap_fraction_wong( return c_bs * f_peak**0.25 * betap * np.sqrt(inverse_aspect) + @staticmethod + def bootstrap_fraction_gi( + betap: float, + pressure_index: float, + temperature_index: float, + inverse_aspect: float, + effective_charge: float, + q95: float, + q0: float, + ) -> float: + """ + Calculate the bootstrap fraction using the Gi et al formula. + + Parameters: + betap (float): Plasma poloidal beta. + pressure_index (float): Pressure profile index. + temperature_index (float): Temperature profile index. + inverse_aspect (float): Inverse aspect ratio. + effective_charge (float): Plasma effective charge. + q95 (float): Safety factor at 95% of the plasma radius. + q0 (float): Safety factor at the magnetic axis. + + Returns: + float: The calculated bootstrap fraction. + + Notes: + + References: + - K. Gi, M. Nakamura, Kenji Tobita, and Y. Ono, “Bootstrap current fraction scaling for a tokamak reactor design study,” + Fusion Engineering and Design, vol. 89, no. 11, pp. 2709–2715, Aug. 2014, + doi: https://doi.org/10.1016/j.fusengdes.2014.07.009. + """ + + # Using the standard variable naming from the Gi et.al. paper + + c_bs = 0.474 * inverse_aspect**-0.1 * pressure_index**0.974 * temperature_index**-0.416 * effective_charge**0.178 * (q95/q0)**-0.133 + + return c_bs * np.sqrt(inverse_aspect) * betap + def fhfac(self, is_): """Function to find H-factor for power balance author: P J Knight, CCFE, Culham Science Centre diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index b0cd7f1d..683bc4b7 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -57,8 +57,11 @@ module current_drive_variables real(dp) :: bscf_wong !! Bootstrap current fraction, Wong et al model + real(dp) :: bscf_gi + !! Bootstrap current fraction, Gi et al model + real(dp) :: cboot - !! bootstrap current fraction multiplier (`i_bootstrap_current=1`) + !! bootstrap current fraction multiplier real(dp) :: cnbeam !! neutral beam current (A) diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index a629224e..40d2dc9c 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -621,7 +621,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 ('i_bootstrap_current') - call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 9, & + call parse_int_variable('i_bootstrap_current', i_bootstrap_current, 1, 10, & '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 c609fb41..41d4a833 100644 --- a/source/fortran/physics_variables.f90 +++ b/source/fortran/physics_variables.f90 @@ -254,6 +254,7 @@ module physics_variables !! - =7 for Andrade et al scaling !! - =8 for Hoang et al scaling !! - =9 for Wong et al scaling + !! - =10 for Gi et al scaling integer :: iculbl !! switch for beta limit scaling (`constraint equation 24`) From 9e57e4c90b4749139402de689776d9a70f6f32bb Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 24 Oct 2024 10:17:07 +0100 Subject: [PATCH 10/10] Add tests and notes for the ARIES bootstrap model --- .../plasma_current/bootstrap_current.md | 2 + process/physics.py | 5 +- tests/unit/test_physics.py | 48 +++++++++++++++++++ 3 files changed, 53 insertions(+), 2 deletions(-) diff --git a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md index 7c36f3e7..ce5438c9 100644 --- a/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/bootstrap_current.md @@ -475,6 +475,8 @@ The model includes the toroidal diamagnetic current in the calculation due to th Is selected by setting `i_bootstrap_current = 6`[^8] +The source reference[^8] does not provide any info about the derivation of the formula. It is only stated like that shown below. + $$ a_1 = 1.10-1.165l_{\text{i}}+0.47l_{\text{i}}^2 $$ diff --git a/process/physics.py b/process/physics.py index 41088f88..63029c6e 100644 --- a/process/physics.py +++ b/process/physics.py @@ -5718,7 +5718,7 @@ def bootstrap_fraction_aries( betap: float, rli: float, core_density: float, - average_desnity: float, + average_density: float, inverse_aspect: float, ) -> float: """ @@ -5735,6 +5735,7 @@ def bootstrap_fraction_aries( float: The calculated bootstrap fraction. Notes: + - The source reference does not provide any info about the derivation of the formula. It is only stated References: - Zoran Dragojlovic et al., “An advanced computational algorithm for systems analysis of tokamak power plants,” @@ -5746,7 +5747,7 @@ def bootstrap_fraction_aries( a_1 = 1.10 - 1.165 * rli + 0.47 * rli**2 b_1 = 0.806 - 0.885 * rli + 0.297 * rli**2 - c_bs = a_1 + b_1 * (core_density / average_desnity) + c_bs = a_1 + b_1 * (core_density / average_density) return c_bs * np.sqrt(inverse_aspect) * betap diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 402a6732..2832702e 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -564,6 +564,54 @@ def test_bootstrap_fraction_sakai(bootstrapfractionsakaiparam, monkeypatch, phys assert bfs == pytest.approx(bootstrapfractionsakaiparam.expected_bfs) +class BootstrapFractionAriesParam(NamedTuple): + betap: Any = None + + rli: Any = None + + core_density: Any = None + + average_density: Any = None + + inverse_aspect: Any = None + + expected_bfs: Any = None + + +@pytest.mark.parametrize( + "bootstrapfractionariesparam", + ( + BootstrapFractionAriesParam( + betap=1.2708883332338736, + rli=1.4279108047138775, + core_density=1.0695994460047332E+20, + average_density=8.1317358967210131E+19, + inverse_aspect=1 / 3, + expected_bfs=4.3237405809568441E-01, + ), + ), +) +def test_bootstrap_fraction_aries(bootstrapfractionariesparam, physics): + """ + Automatically generated Regression Unit Test for bootstrap_fraction_aries. + + This test was generated using data from tests/regression/input_files/large_tokamak.IN.DAT. + + :param bootstrapfractionsauterparam: the data used to mock and assert in this test. + :type bootstrapfractionsauterparam: bootstrapfractionsauterparam + """ + + bfs = physics.bootstrap_fraction_aries( + betap=bootstrapfractionariesparam.betap, + rli=bootstrapfractionariesparam.rli, + core_density=bootstrapfractionariesparam.core_density, + average_density=bootstrapfractionariesparam.average_density, + inverse_aspect=bootstrapfractionariesparam.inverse_aspect, + ) + + assert bfs == pytest.approx(bootstrapfractionariesparam.expected_bfs) + + class PlasmaCurrentParam(NamedTuple): normalised_total_beta: Any = None