🎨 Plot sauter bootstrap profile

process/data_structure/physics_variables.py

@@ -836,6 +836,9 @@
 j_plasma_on_axis: float = None
 """Central plasma current density (A/m2)"""
 
+j_plasma_bootstrap_sauter_profile: list[float] = None
+"""Profile of bootstrap current density in plasma using Sauter et al scaling (A/m2)"""
+
 n_plasma_profile_elements: int = None
 """Number of elements in plasma profile"""
 
@@ -1520,6 +1523,7 @@ def init_physics_variables():
     global pres_plasma_ion_total_profile
     global pres_plasma_fuel_profile
     global j_plasma_on_axis
+    global j_plasma_bootstrap_sauter_profile
     global n_plasma_profile_elements
     global f_dd_branching_trit
     global pden_plasma_alpha_mw
@@ -1779,6 +1783,7 @@ def init_physics_variables():
     pres_plasma_ion_total_profile = []
     pres_plasma_fuel_profile = []
     j_plasma_on_axis = 0.0
+    j_plasma_bootstrap_sauter_profile = []
     n_plasma_profile_elements = 500
     f_dd_branching_trit = 0.0
     pden_plasma_alpha_mw = 0.0

process/io/plot_proc.py

@@ -3963,12 +3963,17 @@ def plot_n_profiles(prof, demo_ranges, mfile_data, scan):
     # ---
 
 
-def plot_jprofile(prof):
+def plot_jprofile(prof, mfile_data, scan):
     """Function to plot density profile
     Arguments:
       prof --> axis object to add plot to
     """
 
+    j_plasma_bootstrap_sauter_profile = [
+        mfile_data.data[f"j_plasma_bootstrap_sauter_profile{i}"].get_scan(scan) / 1000.0
+        for i in range(498)
+    ]
+
     prof.set_xlabel(r"$\rho \quad [r/a]$")
     prof.set_ylabel(r"Current density $[kA/m^2]$")
     prof.set_title("$J$ profile")
@@ -3978,7 +3983,16 @@ def plot_jprofile(prof):
     rho = np.linspace(0, 1)
     y2 = (j_plasma_0 * (1 - rho**2) ** alphaj) / 1e3
 
-    prof.plot(rho, y2, label="$n_{i}$", color="red")
+    prof.plot(rho, y2, color="red")
+
+    prof.plot(
+        np.linspace(0, 1, 498),
+        j_plasma_bootstrap_sauter_profile,
+        label="Sauter Bootstrap",
+        color="green",
+        linestyle="--",
+    )
+    prof.legend()
 
     textstr_j = "\n".join((
         r"$j_0$: " + f"{y2[0]:.3f} kA m$^{{-2}}$\n",
@@ -4004,6 +4018,14 @@ def plot_jprofile(prof):
         ha="left",
         transform=plt.gcf().transFigure,
     )
+    prof.text(
+        0.05,
+        0.02,
+        "*Bootstrap profile is for representation only",
+        fontsize=10,
+        ha="left",
+        transform=plt.gcf().transFigure,
+    )
     prof.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2)
 
 
@@ -12546,7 +12568,7 @@ def main_plot(
     # Plot current density profile
     ax12 = fig4.add_subplot(4, 3, 10)
     ax12.set_position([0.075, 0.105, 0.25, 0.15])
-    plot_jprofile(ax12)
+    plot_jprofile(ax12, m_file_data, scan)
 
     # Plot q profile
     ax13 = fig4.add_subplot(4, 3, 12)

process/physics.py

@@ -2052,9 +2052,12 @@ def physics(self):
             )
         )
 
-        current_drive_variables.f_c_plasma_bootstrap_sauter = (
+        (
+            current_drive_variables.f_c_plasma_bootstrap_sauter,
+            physics_variables.j_plasma_bootstrap_sauter_profile,
+        ) = self.bootstrap_fraction_sauter(self.plasma_profile)
+        current_drive_variables.f_c_plasma_bootstrap_sauter *= (
             current_drive_variables.cboot
-            * self.bootstrap_fraction_sauter(self.plasma_profile)
         )
 
         current_drive_variables.f_c_plasma_bootstrap_sakai = (
@@ -6535,6 +6538,16 @@ def outplas(self):
                 current_drive_variables.f_c_plasma_bootstrap_sauter,
                 "OP ",
             )
+            for point in range(
+                len(physics_variables.j_plasma_bootstrap_sauter_profile)
+            ):
+                po.ovarrf(
+                    self.mfile,
+                    f"Sauter et al bootstrap current density profile at point {point}",
+                    f"(j_plasma_bootstrap_sauter_profile{point})",
+                    physics_variables.j_plasma_bootstrap_sauter_profile[point],
+                    "OP ",
+                )
 
             po.ovarrf(
                 self.outfile,
@@ -7180,9 +7193,6 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
         The code was supplied by Emiliano Fable, IPP Garching (private communication).
         """
 
-        # Number of radial data points along the profile
-        nr = plasma_profile.profile_size
-
         # Radial points from 0 to 1 seperated by 1/profile_size
         roa = plasma_profile.neprofile.profile_x
 
@@ -7217,48 +7227,32 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
             + (physics_variables.q95 - physics_variables.q0) * roa**2
         )
         # Create new array of average mass of fuel portion of ions
-        amain = np.full_like(inverse_q, physics_variables.m_fuel_amu)
+        amain = np.full_like(inverse_q, physics_variables.m_ions_total_amu)
 
         # Create new array of average main ion charge
         zmain = np.full_like(inverse_q, 1.0 + physics_variables.f_plasma_fuel_helium3)
 
-        # Prevent division by zero
-        if ne[nr - 1] == 0.0:
-            ne[nr - 1] = 1e-4 * ne[nr - 2]
-            ni[nr - 1] = 1e-4 * ni[nr - 2]
-
-        # Prevent division by zero
-        if tempe[nr - 1] == 0.0:
-            tempe[nr - 1] = 1e-4 * tempe[nr - 2]
-            tempi[nr - 1] = 1e-4 * tempi[nr - 2]
-
         # Calculate total bootstrap current (MA) by summing along profiles
         # Looping from 2 because _calculate_l31_coefficient() etc should return 0 @ j == 1
-        radial_elements = np.arange(2, nr)
+        radial_elements = np.arange(2, plasma_profile.profile_size)
 
         # Change in localised minor radius to be used as delta term in derivative
         drho = rho[radial_elements] - rho[radial_elements - 1]
 
         # Area of annulus, assuming circular plasma cross-section
         da = 2 * np.pi * rho[radial_elements - 1] * drho  # area of annulus
 
-        # Create the partial derivatives
-        dlogte_drho = (
-            np.log(tempe[radial_elements]) - np.log(tempe[radial_elements - 1])
-        ) / drho
-        dlogti_drho = (
-            np.log(tempi[radial_elements]) - np.log(tempi[radial_elements - 1])
-        ) / drho
-        dlogne_drho = (
-            np.log(ne[radial_elements]) - np.log(ne[radial_elements - 1])
-        ) / drho
+        # Create the partial derivatives using numpy gradient (central differences)
+        dlogte_drho = np.gradient(np.log(tempe), rho)[radial_elements - 1]
+        dlogti_drho = np.gradient(np.log(tempi), rho)[radial_elements - 1]
+        dlogne_drho = np.gradient(np.log(ne), rho)[radial_elements - 1]
 
         jboot = (
             0.5
             * (
                 _calculate_l31_coefficient(
                     radial_elements,
-                    nr,
+                    plasma_profile.profile_size,
                     physics_variables.rmajor,
                     physics_variables.b_plasma_toroidal_on_axis,
                     physics_variables.triang,
@@ -7274,7 +7268,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
                 * dlogne_drho
                 + _calculate_l31_32_coefficient(
                     radial_elements,
-                    nr,
+                    plasma_profile.profile_size,
                     physics_variables.rmajor,
                     physics_variables.b_plasma_toroidal_on_axis,
                     physics_variables.triang,
@@ -7290,7 +7284,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
                 * dlogte_drho
                 + _calculate_l34_alpha_31_coefficient(
                     radial_elements,
-                    nr,
+                    plasma_profile.profile_size,
                     physics_variables.rmajor,
                     physics_variables.b_plasma_toroidal_on_axis,
                     physics_variables.triang,
@@ -7316,7 +7310,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float:
             )
         )  # A/m2
 
-        return np.sum(da * jboot, axis=0) / physics_variables.plasma_current
+        return (np.sum(da * jboot, axis=0) / physics_variables.plasma_current), jboot
 
     @staticmethod
     def bootstrap_fraction_sakai(

tests/unit/test_physics.py

@@ -341,7 +341,7 @@ class BootstrapFractionSauterParam(NamedTuple):
 
     q0: Any = None
 
-    m_fuel_amu: Any = None
+    m_ions_total_amu: Any = None
 
     zeff: Any = None
 
@@ -394,7 +394,7 @@ class BootstrapFractionSauterParam(NamedTuple):
             temp_plasma_ion_vol_avg_kev=12.570861186498382,
             triang=0.5,
             q0=1,
-            m_fuel_amu=2.5,
+            m_ions_total_amu=2.5,
             zeff=2.5211399464385624,
             radius_plasma_pedestal_density_norm=0.9400000000000001,
             b_plasma_toroidal_on_axis=5.326133750416047,
@@ -414,7 +414,7 @@ class BootstrapFractionSauterParam(NamedTuple):
             alphan=1,
             radius_plasma_pedestal_temp_norm=0.9400000000000001,
             alphat=1.45,
-            expected_bfs=0.4110838247346975,
+            expected_bfs=0.4052168782500341,
         ),
     ),
 )
@@ -460,7 +460,9 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph
     monkeypatch.setattr(physics_variables, "q0", bootstrapfractionsauterparam.q0)
 
     monkeypatch.setattr(
-        physics_variables, "m_fuel_amu", bootstrapfractionsauterparam.m_fuel_amu
+        physics_variables,
+        "m_ions_total_amu",
+        bootstrapfractionsauterparam.m_ions_total_amu,
     )
 
     monkeypatch.setattr(
@@ -561,7 +563,7 @@ def test_bootstrap_fraction_sauter(bootstrapfractionsauterparam, monkeypatch, ph
         physics_variables, "alphat", bootstrapfractionsauterparam.alphat
     )
     physics.plasma_profile.run()
-    bfs = physics.bootstrap_fraction_sauter(physics.plasma_profile)
+    bfs, _ = physics.bootstrap_fraction_sauter(physics.plasma_profile)
 
     assert bfs == pytest.approx(bootstrapfractionsauterparam.expected_bfs)