From d3b1b6b7137e1f210f9dbbbe3e8f2b6d9722f5a3 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 6 Nov 2025 14:46:53 +0000 Subject: [PATCH 1/6] Add bootstrap current density profile variable using Sauter scaling --- process/data_structure/physics_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 49f5d21484..58b2be716c 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 From e4ea32bf9da543c160e57faece057c89cc9c6274 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 6 Nov 2025 16:27:38 +0000 Subject: [PATCH 2/6] Enhance bootstrap current density profile handling in physics module and plotting function --- process/data_structure/physics_variables.py | 2 +- process/io/plot_proc.py | 28 ++++++++++++++++++--- process/physics.py | 12 ++++++++- 3 files changed, 37 insertions(+), 5 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 58b2be716c..ef365893bc 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -837,7 +837,7 @@ """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)""" +"""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""" 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..de28feeedf 100644 --- a/process/physics.py +++ b/process/physics.py @@ -6535,6 +6535,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, + "Sauter et al bootstrap current density profile at point", + f"(j_plasma_bootstrap_sauter_profile{point})", + physics_variables.j_plasma_bootstrap_sauter_profile[point], + "OP ", + ) po.ovarrf( self.outfile, @@ -7253,7 +7263,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: np.log(ne[radial_elements]) - np.log(ne[radial_elements - 1]) ) / drho - jboot = ( + jboot = physics_variables.j_plasma_bootstrap_sauter_profile = ( 0.5 * ( _calculate_l31_coefficient( From 8184f9736d2738aef4970beda3fd709f2d2d412e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 6 Nov 2025 16:59:05 +0000 Subject: [PATCH 3/6] Refactor bootstrap current calculations to prevent division by zero and improve partial derivative computation using numpy gradient --- process/physics.py | 24 ++++-------------------- 1 file changed, 4 insertions(+), 20 deletions(-) diff --git a/process/physics.py b/process/physics.py index de28feeedf..8d57a18108 100644 --- a/process/physics.py +++ b/process/physics.py @@ -7232,16 +7232,6 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: # 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) @@ -7252,16 +7242,10 @@ 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 = physics_variables.j_plasma_bootstrap_sauter_profile = ( 0.5 From 1a541e3b25a6c4db53a35ae5e01e7dc35419f70a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 6 Nov 2025 17:05:38 +0000 Subject: [PATCH 4/6] Refactor bootstrap current calculations to use total ion mass and profile size directly, improving clarity and consistency --- process/physics.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/process/physics.py b/process/physics.py index 8d57a18108..d516a19819 100644 --- a/process/physics.py +++ b/process/physics.py @@ -7190,8 +7190,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 @@ -7227,14 +7225,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) # 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] @@ -7252,7 +7250,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: * ( _calculate_l31_coefficient( radial_elements, - nr, + plasma_profile.profile_size, physics_variables.rmajor, physics_variables.b_plasma_toroidal_on_axis, physics_variables.triang, @@ -7268,7 +7266,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, @@ -7284,7 +7282,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, From 5ff48ffb3fde78a79ca8cdd886c11e53416ed78f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Nov 2025 14:09:51 +0000 Subject: [PATCH 5/6] Rename fuel mass variable to total ion mass for clarity in bootstrap fraction calculations --- process/physics.py | 1 - tests/unit/test_physics.py | 10 ++++++---- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/process/physics.py b/process/physics.py index d516a19819..fe08e76b8f 100644 --- a/process/physics.py +++ b/process/physics.py @@ -7190,7 +7190,6 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: The code was supplied by Emiliano Fable, IPP Garching (private communication). """ - # Radial points from 0 to 1 seperated by 1/profile_size roa = plasma_profile.neprofile.profile_x diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 270032a9da..702883ac9d 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( From fe94e801c6ccc21820b24d51eb0aadbb7b74e8f0 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 2 Dec 2025 08:50:24 +0000 Subject: [PATCH 6/6] Requested changes --- process/physics.py | 13 ++++++++----- tests/unit/test_physics.py | 2 +- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/process/physics.py b/process/physics.py index fe08e76b8f..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 = ( @@ -6540,7 +6543,7 @@ def outplas(self): ): po.ovarrf( self.mfile, - "Sauter et al bootstrap current density profile at point", + 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 ", @@ -7244,7 +7247,7 @@ def bootstrap_fraction_sauter(plasma_profile: float) -> float: dlogti_drho = np.gradient(np.log(tempi), rho)[radial_elements - 1] dlogne_drho = np.gradient(np.log(ne), rho)[radial_elements - 1] - jboot = physics_variables.j_plasma_bootstrap_sauter_profile = ( + jboot = ( 0.5 * ( _calculate_l31_coefficient( @@ -7307,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 702883ac9d..4b8fa90aed 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -563,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)