diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 49f5d21484..ef365893bc 100644 --- a/process/data_structure/physics_variables.py +++ b/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 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 22612945c9..3e29e304aa 100644 --- a/process/io/plot_proc.py +++ b/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) diff --git a/process/physics.py b/process/physics.py index d983159135..c93c997650 100644 --- a/process/physics.py +++ b/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,24 +7227,14 @@ 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] @@ -7242,23 +7242,17 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: # 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( diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 270032a9da..4b8fa90aed 100644 --- a/tests/unit/test_physics.py +++ b/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)