From ba4e3b6633cecc51e1b04b72e32d4b9d8d142a76 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 13:29:10 +0100 Subject: [PATCH 01/29] :sparkle: Add b_plasma_toroidal_profile variable and initialize in init_physics_variables function --- 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 c30155c8a4..68c136904e 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -288,6 +288,9 @@ b_plasma_toroidal_on_axis: float = None """Plasma toroidal field on axis (T) (`iteration variable 2`)""" +b_plasma_toroidal_profile: list[float] = None +"""toroidal field profile in plasma (T)""" + b_plasma_total: float = None """Sum of plasma total toroidal + poloidal field (T)""" @@ -1391,6 +1394,7 @@ def init_physics_variables(): global betbm0 global b_plasma_poloidal_average global b_plasma_toroidal_on_axis + global b_plasma_toroidal_profile global b_plasma_total global burnup global burnup_in @@ -1643,6 +1647,7 @@ def init_physics_variables(): betbm0 = 1.5 b_plasma_poloidal_average = 0.0 b_plasma_toroidal_on_axis = 5.68 + b_plasma_toroidal_profile = [] b_plasma_total = 0.0 burnup = 0.0 burnup_in = 0.0 From a76ec210e38c434fbd0cb607eda4a76585e779ce Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 13:41:22 +0100 Subject: [PATCH 02/29] :sparkle: Add calculation and output for toroidal field profile across the plasma --- process/physics.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/process/physics.py b/process/physics.py index ba0df63b78..3ea5186416 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1735,6 +1735,19 @@ def physics(self): + physics_variables.b_plasma_poloidal_average**2 ) + # Calculate the toroidal field across the plasma + # Calculate the toroidal field profile across the plasma (1/R dependence) + rho = np.linspace( + physics_variables.rmajor - physics_variables.rminor, + physics_variables.rmajor + physics_variables.rminor, + physics_variables.n_plasma_profile_elements, + ) + # Avoid division by zero at the magnetic axis + rho = np.where(rho == 0, 1e-10, rho) + physics_variables.b_plasma_toroidal_profile = ( + physics_variables.rmajor * physics_variables.bt / rho + ) + # ============================================ # ----------------------------------------------------- @@ -4214,6 +4227,13 @@ def outplas(self): "(b_plasma_toroidal_on_axis)", physics_variables.b_plasma_toroidal_on_axis, ) + for i in range(len(physics_variables.b_plasma_toroidal_profile)): + po.ovarre( + self.mfile, + f"Toroidal field in plasma at point {i}", + f"b_plasma_toroidal_profile{i}", + physics_variables.b_plasma_toroidal_profile[i], + ) po.ovarrf( self.outfile, "Average poloidal field (T)", From 8629a6bbe301a5f9fd4a2776e040b327a71730c8 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 13:56:31 +0100 Subject: [PATCH 03/29] :sparkle: Add function to plot toroidal magnetic field profiles in plasma --- process/io/plot_proc.py | 59 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index e149229d3d..7128cd8f33 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10695,6 +10695,56 @@ def plot_plasma_poloidal_pressure_contours( axis.set_title("Plasma Poloidal Pressure Contours") +def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): + # Plot magnetic field profiles inside the plasma boundary + + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + + # Get toroidal magnetic field profile (in Tesla) + b_plasma_toroidal_profile = [ + mfile_data.data[f"b_plasma_toroidal_profile{i}"].get_scan(scan) + for i in range(n_plasma_profile_elements) + ] + + # Get major and minor radius for x-axis in metres + rmajor = mfile_data.data["rmajor"].get_scan(scan) + rminor = mfile_data.data["rminor"].get_scan(scan) + + # X-axis: radial position in metres across plasma + x_range_plasma = np.linspace( + rmajor - rminor, rmajor + rminor, len(b_plasma_toroidal_profile) + ) + + # Plot magnetic field first (background) + axis.plot( + x_range_plasma, + b_plasma_toroidal_profile, + color="blue", + label="Toroidal B-field [T]", + ) + + # Plot plasma on top of magnetic field, displaced vertically by bt + plot_plasma(axis, mfile_data, scan, colour_scheme=1) + + # Plot plasma centre dot + axis.plot(rmajor, 0, marker="o", color="red", markersize=8, label="Plasma Centre") + + # Plot vertical lines at plasma edge + axis.axvline( + rmajor - rminor, color="green", linestyle="--", linewidth=1, label="Plasma Edge" + ) + axis.axvline(rmajor + rminor, color="green", linestyle="--", linewidth=1) + axis.set_xlabel("Radial Position [m]") + axis.set_ylabel("Toroidal Magnetic Field [T]") + axis.set_title("Toroidal Magnetic Field Profile in Plasma") + axis.grid(True, linestyle="--", alpha=0.5) + axis.minorticks_on() + axis.legend() + axis.set_xlim(rmajor - 1.25 * rminor, rmajor + 1.25 * rminor) + + def main_plot( fig0, fig1, @@ -10716,6 +10766,7 @@ def main_plot( fig17, fig18, fig19, + fig20, m_file_data, scan, imp="../data/lz_non_corona_14_elements/", @@ -10901,6 +10952,10 @@ def main_plot( fig19.add_subplot(111, aspect="equal"), m_file_data, scan, fig19 ) + plot_magnetic_fields_in_plasma( + fig20.add_subplot(111, aspect="equal"), m_file_data, scan + ) + def main(args=None): # TODO The use of globals here isn't ideal, but is required to get main() @@ -11210,6 +11265,7 @@ def main(args=None): page17 = plt.figure(figsize=(12, 9), dpi=80) page18 = plt.figure(figsize=(12, 9), dpi=80) page19 = plt.figure(figsize=(12, 9), dpi=80) + page20 = plt.figure(figsize=(12, 9), dpi=80) # run main_plot main_plot( @@ -11233,6 +11289,7 @@ def main(args=None): page17, page18, page19, + page20, m_file, scan=scan, demo_ranges=demo_ranges, @@ -11261,6 +11318,7 @@ def main(args=None): pdf.savefig(page17) pdf.savefig(page18) pdf.savefig(page19) + pdf.savefig(page20) # show fig if option used if args.show: @@ -11286,6 +11344,7 @@ def main(args=None): plt.close(page17) plt.close(page18) plt.close(page19) + plt.close(page20) if __name__ == "__main__": From 69db254a0e86efcf26ca2c71f11003bb19ffa6ee Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 14:02:25 +0100 Subject: [PATCH 04/29] :sparkle: Add beta_toroidal_profile variable and initialize in init_physics_variables function --- 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 68c136904e..a945a3d669 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -244,6 +244,9 @@ beta_toroidal: float = None """toroidal beta""" +beta_toroidal_profile: list[float] = None +"""toroidal beta profile""" + beta_thermal: float = None """thermal beta""" @@ -1383,6 +1386,7 @@ def init_physics_variables(): global beta_poloidal global beta_poloidal_eps global beta_toroidal + global beta_toroidal_profile global beta_thermal global beta_thermal_poloidal global beta_thermal_toroidal @@ -1636,6 +1640,7 @@ def init_physics_variables(): beta_poloidal = 0.0 beta_poloidal_eps = 0.0 beta_toroidal = 0.0 + beta_toroidal_profile = [] beta_thermal = 0.0 beta_thermal_poloidal = 0.0 beta_thermal_toroidal = 0.0 From 47bc10dc38b0bec18efb8ba1f1a47419e28ba88d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 14:06:59 +0100 Subject: [PATCH 05/29] :sparkle: Add calculate_plasma_beta method to compute plasma beta from pressure and magnetic field --- process/physics.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/process/physics.py b/process/physics.py index 3ea5186416..7f5b6e1ed8 100644 --- a/process/physics.py +++ b/process/physics.py @@ -3842,6 +3842,22 @@ def calculate_plasma_current( return b_plasma_poloidal_average, qstar, plasma_current + def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float: + """ + Calculate the plasma beta for a given pressure and field. + + :param pres_plasma: Plasma pressure (in Pascals). + :type pres_plasma: float + :param b_field: Magnetic field strength (in Tesla). + :type b_field: float + :returns: The plasma beta (dimensionless). + :rtype: float + + Plasma beta is the ratio of plasma pressure to magnetic pressure. + """ + + return 2 * constants.RMU0 * pres_plasma / (b_field**2) + def outtim(self): po.oheadr(self.outfile, "Times") From a01df3e858fa6933ad9f382184846580bdc8088e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 30 Sep 2025 15:06:10 +0100 Subject: [PATCH 06/29] :sparkle: Add plot_beta_profiles function to visualize beta profiles and update beta_toroidal_profile calculation --- process/io/plot_proc.py | 35 ++++++++++++++++++++++++++++++++++- process/physics.py | 25 ++++++++++++++++++++++++- 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 7128cd8f33..a5ca3cdeb2 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10745,6 +10745,36 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): axis.set_xlim(rmajor - 1.25 * rminor, rmajor + 1.25 * rminor) +def plot_beta_profiles(axis, mfile_data, scan): + # Plot the beta profiles on the given axis + + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + + beta_plasma_toroidal_profile = [ + mfile_data.data[f"beta_toroidal_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) + ] + + axis.plot( + np.linspace(-1, 1, 2 * n_plasma_profile_elements), + beta_plasma_toroidal_profile, + color="blue", + label="Beta Toroidal", + ) + + axis.set_xlabel("$\\rho$ [r/a]") + axis.set_ylabel("$\\beta$ [%]") + axis.minorticks_on() + axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5) + axis.set_title("Beta Profile") + axis.legend() + axis.axvline(x=0, color="black", linestyle="--", linewidth=1) + axis.grid(True, linestyle="--", alpha=0.5) + + + def main_plot( fig0, fig1, @@ -10953,9 +10983,12 @@ def main_plot( ) plot_magnetic_fields_in_plasma( - fig20.add_subplot(111, aspect="equal"), m_file_data, scan + fig20.add_subplot(122, aspect="equal"), m_file_data, scan ) + plot_beta_profiles(fig20.add_subplot(221), m_file_data, scan) + + def main(args=None): # TODO The use of globals here isn't ideal, but is required to get main() diff --git a/process/physics.py b/process/physics.py index 7f5b6e1ed8..471c93e53e 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1737,10 +1737,11 @@ def physics(self): # Calculate the toroidal field across the plasma # Calculate the toroidal field profile across the plasma (1/R dependence) + # Double element size to include both sides of the plasma rho = np.linspace( physics_variables.rmajor - physics_variables.rminor, physics_variables.rmajor + physics_variables.rminor, - physics_variables.n_plasma_profile_elements, + 2 * physics_variables.n_plasma_profile_elements, ) # Avoid division by zero at the magnetic axis rho = np.where(rho == 0, 1e-10, rho) @@ -1760,6 +1761,20 @@ def physics(self): / physics_variables.b_plasma_toroidal_on_axis**2 ) + # Mirror the pressure profiles to match the doubled toroidal field profile + pres_profile_total = np.concatenate([ + physics_variables.pres_plasma_total_profile[::-1], + physics_variables.pres_plasma_total_profile, + ]) + + physics_variables.beta_toroidal_profile = np.array([ + self.calculate_plasma_beta( + pres_plasma=pres_profile_total[i], + b_field=physics_variables.b_plasma_toroidal_profile[i], + ) + for i in range(len(physics_variables.b_plasma_toroidal_profile)) + ]) + # Calculate physics_variables.beta poloidal [-] physics_variables.beta_poloidal = calculate_poloidal_beta( physics_variables.b_plasma_total, @@ -4408,6 +4423,14 @@ def outplas(self): physics_variables.beta_toroidal, "OP ", ) + for i in range(len(physics_variables.beta_toroidal_profile)): + po.ovarre( + self.mfile, + f"Beta toroidal profile at point {i}", + f"beta_toroidal_profile{i}", + physics_variables.beta_toroidal_profile[i], + ) + po.ovarre( self.outfile, "Fast alpha beta", From 975c57dfe79e5d42d38744fa4c07cb419fe6a775 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 20:42:44 +0100 Subject: [PATCH 07/29] :sparkle: Add poloidal field profile variable and initialize in init_physics_variables function --- process/data_structure/physics_variables.py | 5 +++++ process/io/plot_proc.py | 4 +--- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index a945a3d669..6693be77d2 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -294,6 +294,9 @@ b_plasma_toroidal_profile: list[float] = None """toroidal field profile in plasma (T)""" +b_plasma_poloidal_profile: list[float] = None +"""poloidal field profile in plasma (T)""" + b_plasma_total: float = None """Sum of plasma total toroidal + poloidal field (T)""" @@ -1399,6 +1402,7 @@ def init_physics_variables(): global b_plasma_poloidal_average global b_plasma_toroidal_on_axis global b_plasma_toroidal_profile + global b_plasma_poloidal_profile global b_plasma_total global burnup global burnup_in @@ -1653,6 +1657,7 @@ def init_physics_variables(): b_plasma_poloidal_average = 0.0 b_plasma_toroidal_on_axis = 5.68 b_plasma_toroidal_profile = [] + b_plasma_poloidal_profile = [] b_plasma_total = 0.0 burnup = 0.0 burnup_in = 0.0 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index a5ca3cdeb2..88d1f58d35 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10763,7 +10763,7 @@ def plot_beta_profiles(axis, mfile_data, scan): color="blue", label="Beta Toroidal", ) - + axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("$\\beta$ [%]") axis.minorticks_on() @@ -10772,7 +10772,6 @@ def plot_beta_profiles(axis, mfile_data, scan): axis.legend() axis.axvline(x=0, color="black", linestyle="--", linewidth=1) axis.grid(True, linestyle="--", alpha=0.5) - def main_plot( @@ -10987,7 +10986,6 @@ def main_plot( ) plot_beta_profiles(fig20.add_subplot(221), m_file_data, scan) - def main(args=None): From 056770620a7b56a87afd8fc87aacb1610fa8bec8 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 21:44:03 +0100 Subject: [PATCH 08/29] :sparkle: Add calculation and plotting of poloidal magnetic field profile in plasma --- process/io/plot_proc.py | 12 +++++++++ process/physics.py | 56 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 68 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 88d1f58d35..62a58e6666 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10708,6 +10708,11 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): for i in range(n_plasma_profile_elements) ] + b_plasma_poloidal_profile = [ + mfile_data.data[f"b_plasma_poloidal_profile{i}"].get_scan(scan) + for i in range(n_plasma_profile_elements) + ] + # Get major and minor radius for x-axis in metres rmajor = mfile_data.data["rmajor"].get_scan(scan) rminor = mfile_data.data["rminor"].get_scan(scan) @@ -10725,6 +10730,13 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): label="Toroidal B-field [T]", ) + axis.plot( + x_range_plasma, + b_plasma_poloidal_profile, + color="red", + label="Poloidal B-field [T]", + ) + # Plot plasma on top of magnetic field, displaced vertically by bt plot_plasma(axis, mfile_data, scan, colour_scheme=1) diff --git a/process/physics.py b/process/physics.py index 471c93e53e..6894f430e0 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1749,6 +1749,15 @@ def physics(self): physics_variables.rmajor * physics_variables.bt / rho ) + physics_variables.b_plasma_poloidal_profile = ( + self.calculate_poloidal_field_profile( + j_plasma_on_axis=physics_variables.j_plasma_0, + alphaj=physics_variables.alphaj, + n_profile_elements=physics_variables.n_plasma_profile_elements, + rminor=physics_variables.rminor, + ) + ) + # ============================================ # ----------------------------------------------------- @@ -3873,6 +3882,46 @@ def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float: return 2 * constants.RMU0 * pres_plasma / (b_field**2) + def calculate_poloidal_field_profile( + self, + j_plasma_on_axis: float, + alphaj: float, + n_profile_elements: int, + rminor: float, + ) -> np.ndarray: + """ + This method uses Ampère's law and the circular plasma approximation with a parabolic current profile + to compute the poloidal magnetic field profile. + + :param float j_plasma_on_axis: On-axis plasma current density (A/m^2). + :param float alphaj: Current profile index (dimensionless). + :param int n_profile_elements: Number of radial profile elements (half the total number of points across the diameter). + :param float rminor: Plasma minor radius (m). + + :returns: Poloidal magnetic field profile (T) across the full plasma diameter (length 2 * n_profile_elements). + :rtype: numpy.ndarray + + """ + n_profile_elements = n_profile_elements // 2 + rho = np.linspace( + 0, rminor, n_profile_elements + ) # Normalized minor radius (0 to 1) + r = rho * rminor + + j_plasma_profile = j_plasma_on_axis * (1 - (rho / rminor) ** 2) ** alphaj + + integrand = 2 * np.pi * (r * j_plasma_profile) # dI/dr' + dr = r[1] - r[0] + c_enclosed = np.cumsum(integrand * dr) # cumulative current + + # Poloidal field (avoid division by zero at axis) + b_poloidal = np.zeros_like(r) + mask = rho > 0 + b_poloidal[mask] = constants.RMU0 / (2 * np.pi * r[mask]) * c_enclosed[mask] + + # Mirror the b_poloidal profile and concatenate to match the doubled toroidal field profile + return np.concatenate([b_poloidal[::-1], b_poloidal]) + def outtim(self): po.oheadr(self.outfile, "Times") @@ -4265,6 +4314,13 @@ def outplas(self): f"b_plasma_toroidal_profile{i}", physics_variables.b_plasma_toroidal_profile[i], ) + for i in range(len(physics_variables.b_plasma_poloidal_profile)): + po.ovarre( + self.mfile, + f"Poloidal field in plasma at point {i}", + f"b_plasma_poloidal_profile{i}", + physics_variables.b_plasma_poloidal_profile[i], + ) po.ovarrf( self.outfile, "Average poloidal field (T)", From 39748a6494a11ba0296633a23dfe57dc94f0b592 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 21:47:44 +0100 Subject: [PATCH 09/29] :sparkle: Add b_plasma_total_profile variable and initialize in init_physics_variables function --- 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 6693be77d2..2d237a72af 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -297,6 +297,9 @@ b_plasma_poloidal_profile: list[float] = None """poloidal field profile in plasma (T)""" +b_plasma_total_profile: list[float] = None +"""total field profile in plasma (T)""" + b_plasma_total: float = None """Sum of plasma total toroidal + poloidal field (T)""" @@ -1403,6 +1406,7 @@ def init_physics_variables(): global b_plasma_toroidal_on_axis global b_plasma_toroidal_profile global b_plasma_poloidal_profile + global b_plasma_total_profile global b_plasma_total global burnup global burnup_in @@ -1658,6 +1662,7 @@ def init_physics_variables(): b_plasma_toroidal_on_axis = 5.68 b_plasma_toroidal_profile = [] b_plasma_poloidal_profile = [] + b_plasma_total_profile = [] b_plasma_total = 0.0 burnup = 0.0 burnup_in = 0.0 From 63409051cf4a72f3b04bb1cbdae04390728687aa Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 22:49:59 +0100 Subject: [PATCH 10/29] :sparkle: Update magnetic field profile calculations and adjust profile size in PlasmaProfile class --- process/io/plot_proc.py | 25 ++++++++++++++++--------- process/physics.py | 15 ++++++++++++++- process/plasma_profiles.py | 4 ++-- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 62a58e6666..75e169ccab 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10705,38 +10705,45 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): # Get toroidal magnetic field profile (in Tesla) b_plasma_toroidal_profile = [ mfile_data.data[f"b_plasma_toroidal_profile{i}"].get_scan(scan) - for i in range(n_plasma_profile_elements) + for i in range(2 * n_plasma_profile_elements) ] b_plasma_poloidal_profile = [ mfile_data.data[f"b_plasma_poloidal_profile{i}"].get_scan(scan) - for i in range(n_plasma_profile_elements) + for i in range(2 * n_plasma_profile_elements) + ] + + b_plasma_total_profile = [ + mfile_data.data[f"b_plasma_total_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) ] # Get major and minor radius for x-axis in metres rmajor = mfile_data.data["rmajor"].get_scan(scan) rminor = mfile_data.data["rminor"].get_scan(scan) - # X-axis: radial position in metres across plasma - x_range_plasma = np.linspace( - rmajor - rminor, rmajor + rminor, len(b_plasma_toroidal_profile) - ) - # Plot magnetic field first (background) axis.plot( - x_range_plasma, + np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_toroidal_profile)), b_plasma_toroidal_profile, color="blue", label="Toroidal B-field [T]", ) axis.plot( - x_range_plasma, + np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_poloidal_profile)), b_plasma_poloidal_profile, color="red", label="Poloidal B-field [T]", ) + axis.plot( + np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_total_profile)), + b_plasma_total_profile, + color="green", + label="Total B-field [T]", + ) + # Plot plasma on top of magnetic field, displaced vertically by bt plot_plasma(axis, mfile_data, scan, colour_scheme=1) diff --git a/process/physics.py b/process/physics.py index 6894f430e0..4152931460 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1753,11 +1753,17 @@ def physics(self): self.calculate_poloidal_field_profile( j_plasma_on_axis=physics_variables.j_plasma_0, alphaj=physics_variables.alphaj, - n_profile_elements=physics_variables.n_plasma_profile_elements, + n_profile_elements=physics_variables.n_plasma_profile_elements * 2, rminor=physics_variables.rminor, ) ) + # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components + physics_variables.b_plasma_total_profile = np.sqrt( + np.square(physics_variables.b_plasma_toroidal_profile) + + np.square(physics_variables.b_plasma_poloidal_profile) + ) + # ============================================ # ----------------------------------------------------- @@ -4321,6 +4327,13 @@ def outplas(self): f"b_plasma_poloidal_profile{i}", physics_variables.b_plasma_poloidal_profile[i], ) + for i in range(len(physics_variables.b_plasma_total_profile)): + po.ovarre( + self.mfile, + f"Total field in plasma at point {i}", + f"b_plasma_total_profile{i}", + physics_variables.b_plasma_total_profile[i], + ) po.ovarrf( self.outfile, "Average poloidal field (T)", diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index 2d154fa4eb..785a8fa480 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -41,8 +41,8 @@ def __init__(self) -> None: neprofile (NeProfile): An instance of the NeProfile class. teprofile (TeProfile): An instance of the TeProfile class. """ - # Default profile_size = 501, but it's possible to experiment with this value. - self.profile_size = 501 + # Default profile_size = 500, but it's possible to experiment with this value. + self.profile_size = 500 physics_variables.n_plasma_profile_elements = self.profile_size self.outfile = constants.NOUT self.neprofile = profiles.NeProfile(self.profile_size) From ac7b403b9b874f9758b38d25aaa259761969b862 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 22:53:48 +0100 Subject: [PATCH 11/29] Add beta_poloidal_profile variable and initialize in init_physics_variables function --- 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 2d237a72af..19534a06fb 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -247,6 +247,9 @@ beta_toroidal_profile: list[float] = None """toroidal beta profile""" +beta_poloidal_profile: list[float] = None +"""poloidal beta profile""" + beta_thermal: float = None """thermal beta""" @@ -1393,6 +1396,7 @@ def init_physics_variables(): global beta_poloidal_eps global beta_toroidal global beta_toroidal_profile + global beta_poloidal_profile global beta_thermal global beta_thermal_poloidal global beta_thermal_toroidal @@ -1649,6 +1653,7 @@ def init_physics_variables(): beta_poloidal_eps = 0.0 beta_toroidal = 0.0 beta_toroidal_profile = [] + beta_poloidal_profile = [] beta_thermal = 0.0 beta_thermal_poloidal = 0.0 beta_thermal_toroidal = 0.0 From 489cdc7dc495b7dd9624dc5eccd02eef96b9a28f Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 22:54:24 +0100 Subject: [PATCH 12/29] :sparkle: Add beta_total_profile variable and initialize in init_physics_variables function --- process/data_structure/physics_variables.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 19534a06fb..81638437c7 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -250,6 +250,8 @@ beta_poloidal_profile: list[float] = None """poloidal beta profile""" +beta_total_profile: list[float] = None +"""total beta profile""" beta_thermal: float = None """thermal beta""" @@ -1396,6 +1398,7 @@ def init_physics_variables(): global beta_poloidal_eps global beta_toroidal global beta_toroidal_profile + global beta_total_profile global beta_poloidal_profile global beta_thermal global beta_thermal_poloidal @@ -1654,6 +1657,7 @@ def init_physics_variables(): beta_toroidal = 0.0 beta_toroidal_profile = [] beta_poloidal_profile = [] + beta_total_profile = [] beta_thermal = 0.0 beta_thermal_poloidal = 0.0 beta_thermal_toroidal = 0.0 From db52163242765f74f8508428c680c4e18638e8dd Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Tue, 30 Sep 2025 23:06:55 +0100 Subject: [PATCH 13/29] :sparkle: Add beta poloidal and total profile calculations and plotting in plot_beta_profiles function --- process/io/plot_proc.py | 25 +++++++++++++++++++++++++ process/physics.py | 34 ++++++++++++++++++++++++++++++++-- 2 files changed, 57 insertions(+), 2 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 75e169ccab..5641d249cf 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10776,6 +10776,16 @@ def plot_beta_profiles(axis, mfile_data, scan): for i in range(2 * n_plasma_profile_elements) ] + beta_poloidal_profile = [ + mfile_data.data[f"beta_poloidal_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) + ] + + beta_total_profile = [ + mfile_data.data[f"beta_total_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) + ] + axis.plot( np.linspace(-1, 1, 2 * n_plasma_profile_elements), beta_plasma_toroidal_profile, @@ -10783,6 +10793,20 @@ def plot_beta_profiles(axis, mfile_data, scan): label="Beta Toroidal", ) + axis.plot( + np.linspace(-1, 1, 2 * n_plasma_profile_elements), + beta_poloidal_profile, + color="red", + label="Beta Poloidal", + ) + + axis.plot( + np.linspace(-1, 1, 2 * n_plasma_profile_elements), + beta_total_profile, + color="green", + label="Beta Total", + ) + axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("$\\beta$ [%]") axis.minorticks_on() @@ -10791,6 +10815,7 @@ def plot_beta_profiles(axis, mfile_data, scan): axis.legend() axis.axvline(x=0, color="black", linestyle="--", linewidth=1) axis.grid(True, linestyle="--", alpha=0.5) + axis.set_ylim(bottom=0, top=1.5) def main_plot( diff --git a/process/physics.py b/process/physics.py index 4152931460..42b6a88ae1 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1790,6 +1790,22 @@ def physics(self): for i in range(len(physics_variables.b_plasma_toroidal_profile)) ]) + physics_variables.beta_poloidal_profile = np.array([ + self.calculate_plasma_beta( + pres_plasma=pres_profile_total[i], + b_field=physics_variables.b_plasma_poloidal_profile[i], + ) + for i in range(len(physics_variables.b_plasma_poloidal_profile)) + ]) + + physics_variables.beta_total_profile = np.array([ + self.calculate_plasma_beta( + pres_plasma=pres_profile_total[i], + b_field=physics_variables.b_plasma_total_profile[i], + ) + for i in range(len(physics_variables.b_plasma_total_profile)) + ]) + # Calculate physics_variables.beta poloidal [-] physics_variables.beta_poloidal = calculate_poloidal_beta( physics_variables.b_plasma_total, @@ -3885,7 +3901,8 @@ def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float: Plasma beta is the ratio of plasma pressure to magnetic pressure. """ - + if b_field == 0: + return 0.0 return 2 * constants.RMU0 * pres_plasma / (b_field**2) def calculate_poloidal_field_profile( @@ -4499,7 +4516,20 @@ def outplas(self): f"beta_toroidal_profile{i}", physics_variables.beta_toroidal_profile[i], ) - + for i in range(len(physics_variables.beta_poloidal_profile)): + po.ovarre( + self.mfile, + f"Beta poloidal profile at point {i}", + f"beta_poloidal_profile{i}", + physics_variables.beta_poloidal_profile[i], + ) + for i in range(len(physics_variables.beta_total_profile)): + po.ovarre( + self.mfile, + f"Beta total profile at point {i}", + f"beta_total_profile{i}", + physics_variables.beta_total_profile[i], + ) po.ovarre( self.outfile, "Fast alpha beta", From f33a1178545ec786085be65a70110cf50aede16e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 2 Oct 2025 11:14:31 +0100 Subject: [PATCH 14/29] :sparkle: Refactor beta profile plotting: separate poloidal profile into its own function --- process/io/plot_proc.py | 43 ++++++++++++++++++++++++++++++----------- process/physics.py | 4 ++-- 2 files changed, 34 insertions(+), 13 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 5641d249cf..dc10ea1a6c 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10776,10 +10776,6 @@ def plot_beta_profiles(axis, mfile_data, scan): for i in range(2 * n_plasma_profile_elements) ] - beta_poloidal_profile = [ - mfile_data.data[f"beta_poloidal_profile{i}"].get_scan(scan) - for i in range(2 * n_plasma_profile_elements) - ] beta_total_profile = [ mfile_data.data[f"beta_total_profile{i}"].get_scan(scan) @@ -10793,12 +10789,6 @@ def plot_beta_profiles(axis, mfile_data, scan): label="Beta Toroidal", ) - axis.plot( - np.linspace(-1, 1, 2 * n_plasma_profile_elements), - beta_poloidal_profile, - color="red", - label="Beta Poloidal", - ) axis.plot( np.linspace(-1, 1, 2 * n_plasma_profile_elements), @@ -10815,7 +10805,37 @@ def plot_beta_profiles(axis, mfile_data, scan): axis.legend() axis.axvline(x=0, color="black", linestyle="--", linewidth=1) axis.grid(True, linestyle="--", alpha=0.5) - axis.set_ylim(bottom=0, top=1.5) + +def plot_beta_profiles_poloidal(axis, mfile_data, scan): + # Plot the beta profiles on the given axis + + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + + beta_poloidal_profile = [ + mfile_data.data[f"beta_poloidal_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) + ] + + + axis.plot( + np.linspace(-1, 1, 2 * n_plasma_profile_elements), + beta_poloidal_profile, + color="red", + label="Beta Poloidal", + ) + + + axis.set_xlabel("$\\rho$ [r/a]") + axis.set_ylabel("$\\beta$ [%]") + axis.minorticks_on() + axis.grid(which="minor", linestyle=":", linewidth=0.5, alpha=0.5) + axis.set_title("Beta Profile") + axis.legend() + axis.set_yscale("log") + axis.axvline(x=0, color="black", linestyle="--", linewidth=1) + axis.grid(True, linestyle="--", alpha=0.5) def main_plot( @@ -11030,6 +11050,7 @@ def main_plot( ) plot_beta_profiles(fig20.add_subplot(221), m_file_data, scan) + plot_beta_profiles_poloidal(fig20.add_subplot(223), m_file_data, scan) def main(args=None): diff --git a/process/physics.py b/process/physics.py index 42b6a88ae1..7e9a00c84d 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1746,12 +1746,12 @@ def physics(self): # Avoid division by zero at the magnetic axis rho = np.where(rho == 0, 1e-10, rho) physics_variables.b_plasma_toroidal_profile = ( - physics_variables.rmajor * physics_variables.bt / rho + physics_variables.rmajor * physics_variables.b_plasma_toroidal_on_axis / rho ) physics_variables.b_plasma_poloidal_profile = ( self.calculate_poloidal_field_profile( - j_plasma_on_axis=physics_variables.j_plasma_0, + j_plasma_on_axis=physics_variables.j_plasma_on_axis, alphaj=physics_variables.alphaj, n_profile_elements=physics_variables.n_plasma_profile_elements * 2, rminor=physics_variables.rminor, From f688fab35fba05be4c24f703aa42e319ca67f86f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 15:26:44 +0100 Subject: [PATCH 15/29] Calculate new mean poloidal field --- process/io/plot_proc.py | 9 +++------ process/physics.py | 6 +++++- process/plasma_profiles.py | 2 +- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index dc10ea1a6c..f36456fa0e 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10776,7 +10776,6 @@ def plot_beta_profiles(axis, mfile_data, scan): for i in range(2 * n_plasma_profile_elements) ] - beta_total_profile = [ mfile_data.data[f"beta_total_profile{i}"].get_scan(scan) for i in range(2 * n_plasma_profile_elements) @@ -10789,7 +10788,6 @@ def plot_beta_profiles(axis, mfile_data, scan): label="Beta Toroidal", ) - axis.plot( np.linspace(-1, 1, 2 * n_plasma_profile_elements), beta_total_profile, @@ -10805,7 +10803,8 @@ def plot_beta_profiles(axis, mfile_data, scan): axis.legend() axis.axvline(x=0, color="black", linestyle="--", linewidth=1) axis.grid(True, linestyle="--", alpha=0.5) - + + def plot_beta_profiles_poloidal(axis, mfile_data, scan): # Plot the beta profiles on the given axis @@ -10818,7 +10817,6 @@ def plot_beta_profiles_poloidal(axis, mfile_data, scan): for i in range(2 * n_plasma_profile_elements) ] - axis.plot( np.linspace(-1, 1, 2 * n_plasma_profile_elements), beta_poloidal_profile, @@ -10826,7 +10824,6 @@ def plot_beta_profiles_poloidal(axis, mfile_data, scan): label="Beta Poloidal", ) - axis.set_xlabel("$\\rho$ [r/a]") axis.set_ylabel("$\\beta$ [%]") axis.minorticks_on() @@ -10835,7 +10832,7 @@ def plot_beta_profiles_poloidal(axis, mfile_data, scan): axis.legend() axis.set_yscale("log") axis.axvline(x=0, color="black", linestyle="--", linewidth=1) - axis.grid(True, linestyle="--", alpha=0.5) + axis.grid(True, linestyle="--", alpha=0.5) def main_plot( diff --git a/process/physics.py b/process/physics.py index 7e9a00c84d..6452a17db9 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1617,7 +1617,7 @@ def physics(self): # Calculate plasma current ( - physics_variables.b_plasma_poloidal_average, + _, physics_variables.qstar, physics_variables.plasma_current, ) = self.calculate_plasma_current( @@ -1757,6 +1757,10 @@ def physics(self): rminor=physics_variables.rminor, ) ) + # Calculate the average poloidal field from the profile + physics_variables.b_plasma_poloidal_average = np.mean( + physics_variables.b_plasma_poloidal_profile + ) # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components physics_variables.b_plasma_total_profile = np.sqrt( diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index 785a8fa480..09c33f218b 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -42,7 +42,7 @@ def __init__(self) -> None: teprofile (TeProfile): An instance of the TeProfile class. """ # Default profile_size = 500, but it's possible to experiment with this value. - self.profile_size = 500 + self.profile_size = 501 physics_variables.n_plasma_profile_elements = self.profile_size self.outfile = constants.NOUT self.neprofile = profiles.NeProfile(self.profile_size) From 138d9e093bfcdcc8b6723175411773abed01c7d9 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 15:32:28 +0100 Subject: [PATCH 16/29] :sparkle: Add b_plasma_poloidal_edge variable and initialize in init_physics_variables function --- process/data_structure/physics_variables.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 81638437c7..937e99a856 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -292,6 +292,8 @@ b_plasma_poloidal_average: float = None """Plasma average poloidal field (T)""" +b_plasma_poloidal_edge: float = None +"""Plasma poloidal field at plasma edge (T)""" b_plasma_toroidal_on_axis: float = None """Plasma toroidal field on axis (T) (`iteration variable 2`)""" @@ -1410,6 +1412,7 @@ def init_physics_variables(): global beta_norm_toroidal global betbm0 global b_plasma_poloidal_average + global b_plasma_poloidal_edge global b_plasma_toroidal_on_axis global b_plasma_toroidal_profile global b_plasma_poloidal_profile @@ -1668,6 +1671,7 @@ def init_physics_variables(): beta_norm_toroidal = 0.0 betbm0 = 1.5 b_plasma_poloidal_average = 0.0 + b_plasma_poloidal_edge = 0.0 b_plasma_toroidal_on_axis = 5.68 b_plasma_toroidal_profile = [] b_plasma_poloidal_profile = [] From 7a12caef7074076a173d13ebe5370503896936a8 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 16:42:50 +0100 Subject: [PATCH 17/29] :memo: Update plasma beta documentation: enhance definitions and clarify volume averaged beta equations --- .../physics-models/plasma_beta/plasma_beta.md | 26 ++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/documentation/physics-models/plasma_beta/plasma_beta.md b/documentation/physics-models/plasma_beta/plasma_beta.md index d49bff6c36..1ecc214ae0 100644 --- a/documentation/physics-models/plasma_beta/plasma_beta.md +++ b/documentation/physics-models/plasma_beta/plasma_beta.md @@ -5,9 +5,19 @@ The efficiency of confinement of plasma pressure by the magnetic field is represented by the ratio: $$ -\beta = \frac{2\mu_0p}{B^2} +\beta(\rho) = \frac{2\mu_0p(\rho)}{\left(B(\rho)\right)^2} $$ +Where $\beta$ is a function of normalised minor radius across the plasma ($\rho$), due to the change in pressure and magnetic field strength. + +The standard $\beta$ term used for comparison and to represent the plasma as a whole in many calculations is the volume averaged value given by: + +$$ +\langle \beta \rangle = \frac{2\mu_0 \langle p \rangle}{\langle B \rangle^2} +$$ + +Where $\langle p \rangle$ is the volume averaged plasma pressure and $\langle B \rangle$ is the average field. + There are several different measures of this type, arising from different choices of definition and from the need to quantify different equilibrium properties. In its expanded form of magnetic field components: @@ -92,6 +102,20 @@ $$ ------------------------ +## Definitions + +### Volume averaged thermal toroidal beta + +$$ +\overbrace{\langle \beta_t \rangle}^{\texttt{beta_toroidal_thermal_vol_avg}} = \frac{2\mu_0 \langle p_{\text{thermal}} \rangle}{\underbrace{B_{\text{T,on-axis}}^2}_{\texttt{b_plasma_toroidal_on_axis}}} +$$ + +### Volume averaged thermal poloidal beta + +$$ +\langle \beta_p \rangle = \frac{2\mu_0 \langle p_{\text{thermal}} \rangle}{\langle B_{\text{P,average}} \rangle^2} +$$ + ## Troyon Beta Limit The Troyon plasma beta limit is given by[^0][^1]: From a9944d2745bf5957a64564cc66bc2250283c1e3e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 16:45:50 +0100 Subject: [PATCH 18/29] =?UTF-8?q?=F0=9F=94=84=20Rename=20volume=20averaged?= =?UTF-8?q?=20plasma=20pressure=20variable=20to=20thermal=20plasma=20press?= =?UTF-8?q?ure=20across=20multiple=20files?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/data/large_tokamak_nof_2_MFILE.DAT | 2 +- examples/data/large_tokamak_nof_MFILE.DAT | 2 +- process/data_structure/physics_variables.py | 4 ++-- process/physics.py | 8 ++++---- process/plasma_profiles.py | 2 +- tests/integration/data/large_tokamak_MFILE.DAT | 2 +- 6 files changed, 10 insertions(+), 10 deletions(-) diff --git a/examples/data/large_tokamak_nof_2_MFILE.DAT b/examples/data/large_tokamak_nof_2_MFILE.DAT index e1d77b2528..db590328c7 100644 --- a/examples/data/large_tokamak_nof_2_MFILE.DAT +++ b/examples/data/large_tokamak_nof_2_MFILE.DAT @@ -232,7 +232,7 @@ Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.02466460344915771e+20 OP Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 8.63127495812675338e+19 OP Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 7.82668759018785902e+05 OP - Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.26860509860517632e+05 OP + Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.26860509860517632e+05 OP Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19833068191158598e+00 OP Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.06084473006886912e+19 OP Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.58453241268011663e+19 OP diff --git a/examples/data/large_tokamak_nof_MFILE.DAT b/examples/data/large_tokamak_nof_MFILE.DAT index cdd794dd89..6a1b0b18b1 100644 --- a/examples/data/large_tokamak_nof_MFILE.DAT +++ b/examples/data/large_tokamak_nof_MFILE.DAT @@ -304,7 +304,7 @@ Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.11814263037413868e+20 OP Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 9.41418090734859715e+19 OP Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 9.10652371594472905e+05 OP - Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.63957209157818230e+05 OP + Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.63957209157818230e+05 OP Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19999999999999951e+00 OP Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.34892735750798868e+19 OP Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.49714589130413015e+19 OP diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 937e99a856..fa823fb8da 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -843,8 +843,8 @@ n_plasma_profile_elements: int = None """Number of elements in plasma profile""" -pres_plasma_vol_avg: float = None -"""Volume averaged plasma pressure (Pa)""" +pres_plasma_thermal_vol_avg: float = None +"""Volume averaged thermal plasma pressure (Pa)""" f_dd_branching_trit: float = None diff --git a/process/physics.py b/process/physics.py index 6452a17db9..1282b65308 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2075,7 +2075,7 @@ def physics(self): * self.bootstrap_fraction_andrade( beta_poloidal=physics_variables.beta_poloidal, core_pressure=physics_variables.pres_plasma_on_axis, - average_pressure=physics_variables.pres_plasma_vol_avg, + average_pressure=physics_variables.pres_plasma_thermal_vol_avg, inverse_aspect=physics_variables.eps, ) ) @@ -2666,7 +2666,7 @@ def physics(self): self.calculate_beta_norm_max_thloreus( c_beta=physics_variables.c_beta, pres_plasma_on_axis=physics_variables.pres_plasma_on_axis, - pres_plasma_vol_avg=physics_variables.pres_plasma_vol_avg, + pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, ) ) @@ -4774,8 +4774,8 @@ def outplas(self): po.ovarre( self.outfile, "Volume averaged plasma pressure (Pa)", - "(pres_plasma_vol_avg)", - physics_variables.pres_plasma_vol_avg, + "(pres_plasma_thermal_vol_avg)", + physics_variables.pres_plasma_thermal_vol_avg, "OP ", ) # As array output is not currently supported, each element is output as a float instance diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index 09c33f218b..e18ce0da37 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -320,7 +320,7 @@ def calculate_profile_factors(self) -> None: # Shall assume that the pressure profile is parabolic. Can find volume average from # profile index and core value the same as for density and temperature - physics_variables.pres_plasma_vol_avg = ( + physics_variables.pres_plasma_thermal_vol_avg = ( physics_variables.pres_plasma_on_axis / (physics_variables.alphap + 1) ) diff --git a/tests/integration/data/large_tokamak_MFILE.DAT b/tests/integration/data/large_tokamak_MFILE.DAT index 1cf005fc1a..259813c18a 100644 --- a/tests/integration/data/large_tokamak_MFILE.DAT +++ b/tests/integration/data/large_tokamak_MFILE.DAT @@ -466,7 +466,7 @@ Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.05329687553494565e+20 OP Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_____________ 8.86821331867437138e+19 OP Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 8.42056568743616925e+05 OP - Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_vol_avg)_____________ 2.44074367751773010e+05 OP + Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.44074367751773010e+05 OP Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19999998892791648e+00 OP Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.10762537507155067e+19 OP Fuel_ion_number_density_(/m3)____________________________________________ (nd_fuel_ions)_________________ 6.43726858733289964e+19 OP From f10387f744d9c0f50c7f85e7fdcfddb8289cd9bb Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 16:49:01 +0100 Subject: [PATCH 19/29] =?UTF-8?q?=F0=9F=94=84=20Rename=20plasma=20pressure?= =?UTF-8?q?=20variable=20to=20thermal=20plasma=20pressure=20in=20multiple?= =?UTF-8?q?=20files?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- examples/data/large_tokamak_nof_2_MFILE.DAT | 2 +- examples/data/large_tokamak_nof_MFILE.DAT | 2 +- process/data_structure/physics_variables.py | 6 +++--- process/physics.py | 8 ++++---- process/plasma_profiles.py | 15 ++++++++------- tests/unit/test_plasma_profiles.py | 12 ++++++------ 6 files changed, 23 insertions(+), 22 deletions(-) diff --git a/examples/data/large_tokamak_nof_2_MFILE.DAT b/examples/data/large_tokamak_nof_2_MFILE.DAT index db590328c7..8920191321 100644 --- a/examples/data/large_tokamak_nof_2_MFILE.DAT +++ b/examples/data/large_tokamak_nof_2_MFILE.DAT @@ -231,7 +231,7 @@ Volume_averaged_electron_number_density_(/m3)____________________________ (nd_plasma_electrons_vol_avg)_________________________ 7.79622390002983731e+19 Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.02466460344915771e+20 OP Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 8.63127495812675338e+19 OP - Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 7.82668759018785902e+05 OP + Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_thermal_on_axis)___________________________ 7.82668759018785902e+05 OP Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.26860509860517632e+05 OP Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19833068191158598e+00 OP Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.06084473006886912e+19 OP diff --git a/examples/data/large_tokamak_nof_MFILE.DAT b/examples/data/large_tokamak_nof_MFILE.DAT index 6a1b0b18b1..0e3791a63a 100644 --- a/examples/data/large_tokamak_nof_MFILE.DAT +++ b/examples/data/large_tokamak_nof_MFILE.DAT @@ -303,7 +303,7 @@ Volume_averaged_electron_number_density_(/m3)____________________________ (nd_plasma_electrons_vol_avg)_________________________ 8.50078961768999158e+19 Electron_number_density_on_axis_(/m3)____________________________________ (nd_plasma_electron_on_axis)__________________________ 1.11814263037413868e+20 OP Line-averaged_electron_number_density_(/m3)______________________________ (nd_electron_line)_________________________ 9.41418090734859715e+19 OP - Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_on_axis)___________________________ 9.10652371594472905e+05 OP + Plasma_pressure_on_axis_(Pa)_____________________________________________ (pres_plasma_thermal_on_axis)___________________________ 9.10652371594472905e+05 OP Volume_averaged_plasma_pressure_(Pa)_____________________________________ (pres_plasma_thermal_vol_avg)_____________ 2.63957209157818230e+05 OP Line-averaged_electron_density_/_Greenwald_density_______________________ (dnla_gw)______________________ 1.19999999999999951e+00 OP Total_Ion_number_density_(/m3)___________________________________________ (nd_ions_total)________________ 7.34892735750798868e+19 OP diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index fa823fb8da..bcf5ca6dae 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -822,7 +822,7 @@ """margin to vertical stability""" -pres_plasma_on_axis: float = None +pres_plasma_thermal_on_axis: float = None """central total plasma pressure (Pa)""" pres_plasma_total_profile: list[float] = None @@ -1514,7 +1514,7 @@ def init_physics_variables(): global nd_plasma_electron_on_axis global nd_plasma_ions_on_axis global m_s_limit - global pres_plasma_on_axis + global pres_plasma_thermal_on_axis global pres_plasma_total_profile global pres_plasma_electron_profile global pres_plasma_ion_total_profile @@ -1773,7 +1773,7 @@ def init_physics_variables(): nd_plasma_electron_on_axis = 0.0 nd_plasma_ions_on_axis = 0.0 m_s_limit = 0.3 - pres_plasma_on_axis = 0.0 + pres_plasma_thermal_on_axis = 0.0 pres_plasma_total_profile = [] pres_plasma_electron_profile = [] pres_plasma_ion_total_profile = [] diff --git a/process/physics.py b/process/physics.py index 1282b65308..fed67bd64a 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1628,7 +1628,7 @@ def physics(self): physics_variables.i_plasma_current, physics_variables.kappa, physics_variables.kappa95, - physics_variables.pres_plasma_on_axis, + physics_variables.pres_plasma_thermal_on_axis, physics_variables.len_plasma_poloidal, physics_variables.q95, physics_variables.rmajor, @@ -2074,7 +2074,7 @@ def physics(self): current_drive_variables.cboot * self.bootstrap_fraction_andrade( beta_poloidal=physics_variables.beta_poloidal, - core_pressure=physics_variables.pres_plasma_on_axis, + core_pressure=physics_variables.pres_plasma_thermal_on_axis, average_pressure=physics_variables.pres_plasma_thermal_vol_avg, inverse_aspect=physics_variables.eps, ) @@ -4767,8 +4767,8 @@ def outplas(self): po.ovarre( self.outfile, "Plasma pressure on axis (Pa)", - "(pres_plasma_on_axis)", - physics_variables.pres_plasma_on_axis, + "(pres_plasma_thermal_on_axis)", + physics_variables.pres_plasma_thermal_on_axis, "OP ", ) po.ovarre( diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index e18ce0da37..ddd6a798e9 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -27,7 +27,7 @@ class PlasmaProfile: parameterise_plasma(): Initializes the density and temperature profile averages and peak values. parabolic_paramterisation(): Parameterizes plasma profiles in the case where i_plasma_pedestal=0. pedestal_parameterisation(): Instance temperature and density profiles then integrate them, setting physics variables temp_plasma_electron_density_weighted_kev and temp_plasma_ion_density_weighted_kev. - calculate_profile_factors(): Calculate and set the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap). + calculate_profile_factors(): Calculate and set the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap). calculate_parabolic_profile_factors(): The gradient information for i_plasma_pedestal = 0. """ @@ -255,10 +255,10 @@ def pedestal_parameterisation(self) -> None: def calculate_profile_factors(self) -> None: """ - Calculate and set the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap). + Calculate and set the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap). - This method calculates the central pressure (pres_plasma_on_axis) using the ideal gas law and the pressure profile index (alphap). - It sets the value of the physics variable `pres_plasma_on_axis`. + This method calculates the central pressure (pres_plasma_thermal_on_axis) using the ideal gas law and the pressure profile index (alphap). + It sets the value of the physics variable `pres_plasma_thermal_on_axis`. Args: None @@ -269,7 +269,7 @@ def calculate_profile_factors(self) -> None: # Central pressure (Pa), from ideal gas law : p = nkT - physics_variables.pres_plasma_on_axis = ( + physics_variables.pres_plasma_thermal_on_axis = ( ( physics_variables.nd_plasma_electron_on_axis * physics_variables.temp_plasma_electron_on_axis_kev @@ -312,7 +312,7 @@ def calculate_profile_factors(self) -> None: ) # Pressure profile index (N.B. no pedestal effects included here) - # N.B. pres_plasma_on_axis is NOT equal to

* (1 + alphap), but p(rho) = n(rho)*T(rho) + # N.B. pres_plasma_thermal_on_axis is NOT equal to

* (1 + alphap), but p(rho) = n(rho)*T(rho) # and

= .T_n where <...> denotes volume-averages and T_n is the # density-weighted temperature @@ -321,7 +321,8 @@ def calculate_profile_factors(self) -> None: # Shall assume that the pressure profile is parabolic. Can find volume average from # profile index and core value the same as for density and temperature physics_variables.pres_plasma_thermal_vol_avg = ( - physics_variables.pres_plasma_on_axis / (physics_variables.alphap + 1) + physics_variables.pres_plasma_thermal_on_axis + / (physics_variables.alphap + 1) ) # Central plasma current density (A/m^2) diff --git a/tests/unit/test_plasma_profiles.py b/tests/unit/test_plasma_profiles.py index 4c5f633116..37dbe8e27a 100644 --- a/tests/unit/test_plasma_profiles.py +++ b/tests/unit/test_plasma_profiles.py @@ -141,7 +141,7 @@ class PlasmaProfilesParam(NamedTuple): temp_plasma_electron_on_axis_kev: float = 0.0 - pres_plasma_on_axis: float = 0.0 + pres_plasma_thermal_on_axis: float = 0.0 nd_plasma_separatrix_electron: float = 0.0 @@ -225,7 +225,7 @@ class PlasmaProfilesParam(NamedTuple): alphap=0.0, tbeta=2, temp_plasma_electron_on_axis_kev=0.0, - pres_plasma_on_axis=0.0, + pres_plasma_thermal_on_axis=0.0, nd_plasma_separatrix_electron=3.6421334486704804e19, temp_plasma_separatrix_kev=0.10000000000000001, pcoef=0.0, @@ -270,7 +270,7 @@ class PlasmaProfilesParam(NamedTuple): alphap=2.4500000000000002, tbeta=2, temp_plasma_electron_on_axis_kev=27.369013322953624, - pres_plasma_on_axis=868071.46874220832, + pres_plasma_thermal_on_axis=868071.46874220832, nd_plasma_separatrix_electron=3.6421334486704804e19, temp_plasma_separatrix_kev=0.10000000000000001, pcoef=1.1110842637642833, @@ -355,8 +355,8 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch): monkeypatch.setattr( physics_variables, - "pres_plasma_on_axis", - plasmaprofilesparam.pres_plasma_on_axis, + "pres_plasma_thermal_on_axis", + plasmaprofilesparam.pres_plasma_thermal_on_axis, ) monkeypatch.setattr( @@ -482,7 +482,7 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch): plasmaprofilesparam.expected_te0 ) - assert physics_variables.pres_plasma_on_axis == pytest.approx( + assert physics_variables.pres_plasma_thermal_on_axis == pytest.approx( plasmaprofilesparam.expected_p0 ) From bdfee7772dd797e8a67b85a4e403082089d6ed27 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:03:43 +0100 Subject: [PATCH 20/29] :sparkle: Add visualization for plasma outer circular flux surface in magnetic field plot --- process/io/plot_proc.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index f36456fa0e..ba5a03964a 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10743,6 +10743,19 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): color="green", label="Total B-field [T]", ) + + # Plot plasma minor radius as a circle (ensure it's visible and not covered) + circle = plt.Circle( + (rmajor, 0), + rminor, + color="purple", + fill=False, + linewidth=2, + linestyle="--", + label="Plasma Outer Circlular Flux Surface", + zorder=10, + ) + axis.add_patch(circle) # Plot plasma on top of magnetic field, displaced vertically by bt plot_plasma(axis, mfile_data, scan, colour_scheme=1) From 7f740bd2978ebfe03461a2e38d0b70a8d11ddaac Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:06:35 +0100 Subject: [PATCH 21/29] :sparkle: Add c_plasma_circular variable for circular plasma current and initialize in init_physics_variables function --- 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 bcf5ca6dae..3bea28006f 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -970,6 +970,9 @@ plasma_current: float = None """plasma current (A)""" +c_plasma_circular: float = None +"""Plasma current for circular plasma (A)""" + p_plasma_neutron_mw: float = None """Neutron fusion power from just the plasma [MW]""" @@ -1552,6 +1555,7 @@ def init_physics_variables(): global pflux_fw_rad_mw global pden_ion_electron_equilibration_mw global plasma_current + global c_plasma_circular global p_plasma_neutron_mw global p_neutron_total_mw global pden_neutron_total_mw @@ -1811,6 +1815,7 @@ def init_physics_variables(): pflux_fw_rad_mw = 0.0 pden_ion_electron_equilibration_mw = 0.0 plasma_current = 0.0 + c_plasma_circular = 0.0 p_plasma_neutron_mw = 0.0 p_neutron_total_mw = 0.0 pden_neutron_total_mw = 0.0 From 4d73c88bd6643a060ce1dd39236854e2156a1707 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:14:54 +0100 Subject: [PATCH 22/29] =?UTF-8?q?=F0=9F=94=84=20Update=20plasma=20pressure?= =?UTF-8?q?=20variable=20to=20thermal=20plasma=20pressure=20and=20add=20ca?= =?UTF-8?q?lculation=20for=20circular=20plasma=20current?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/io/plot_proc.py | 2 +- process/physics.py | 15 ++++++++++++++- 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index ba5a03964a..ceee730e74 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10743,7 +10743,7 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): color="green", label="Total B-field [T]", ) - + # Plot plasma minor radius as a circle (ensure it's visible and not covered) circle = plt.Circle( (rmajor, 0), diff --git a/process/physics.py b/process/physics.py index fed67bd64a..92d5fd3726 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2665,7 +2665,7 @@ def physics(self): physics_variables.beta_norm_max_thloreus = ( self.calculate_beta_norm_max_thloreus( c_beta=physics_variables.c_beta, - pres_plasma_on_axis=physics_variables.pres_plasma_on_axis, + pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, ) ) @@ -3857,6 +3857,12 @@ def calculate_plasma_current( * b_plasma_toroidal_on_axis ) # i_plasma_current == 2 case covered above + physics_variables.c_plasma_circular = ( + (constants.TWOPI / constants.RMU0) + * rminor**2 + / (rmajor * q95) + * b_plasma_toroidal_on_axis + ) # Calculate cyclindrical safety factor from IPDG89 qstar = ( @@ -4286,6 +4292,13 @@ def outplas(self): physics_variables.plasma_current / 1.0e6, "OP ", ) + po.ovarrf( + self.outfile, + "Plasma current for a circular plasma (A)", + "(c_plasma_circular)", + physics_variables.c_plasma_circular, + "OP ", + ) if physics_variables.i_alphaj == 1: po.ovarrf( From e5e3b70904c743e2c97b799598bf7db52d282a4b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:28:34 +0100 Subject: [PATCH 23/29] :sparkle: Implement calculation for on-axis circular plasma current density and update plotting function to visualize it --- process/io/plot_proc.py | 34 +++++++++++++++++++++------------- process/physics.py | 7 +++++++ process/plasma_profiles.py | 11 +++++++++++ 3 files changed, 39 insertions(+), 13 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index ceee730e74..cc8a091261 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -3706,40 +3706,48 @@ def plot_n_profiles(prof, demo_ranges, mfile_data, scan): # --- -def plot_jprofile(prof): +def plot_jprofile(axis, mfile_data, scan): """Function to plot density profile Arguments: prof --> axis object to add plot to """ - prof.set_xlabel(r"$\rho \quad [r/a]$") - prof.set_ylabel(r"Current density $[kA/m^2]$") - prof.set_title("$J$ profile") - prof.minorticks_on() - prof.set_xlim([0, 1.0]) + axis.set_xlabel(r"$\rho \quad [r/a]$") + axis.set_ylabel(r"Current density $[kA/m^2]$") + axis.set_title("$J$ profile") + axis.minorticks_on() + axis.set_xlim([0, 1.0]) + + j_plasma_circular_on_axis = mfile_data.data["j_plasma_circular_on_axis"].get_scan( + scan + ) rho = np.linspace(0, 1) y2 = (j_plasma_0 * (1 - rho**2) ** alphaj) / 1e3 + y_circular = (j_plasma_circular_on_axis * (1 - rho**2) ** alphaj) / 1e3 - prof.plot(rho, y2, label="$n_{i}$", color="red") + axis.plot(rho, y2, label="Real", color="red") + axis.plot(rho, y_circular, label="Circular", color="blue", linestyle="--") + axis.legend() textstr_j = "\n".join(( - r"$j_0$: " + f"{y2[0]:.3f} kA m$^{{-2}}$\n", + r"$j_0$: " + f"{j_plasma_0 / 1000:.3f} kA m$^{{-2}}$\n", + r"$j_0,circular$: " + f"{j_plasma_circular_on_axis / 1000:.3f} kA m$^{{-2}}$\n", r"$\alpha_J$: " + f"{alphaj:.3f}", )) props_j = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5} - prof.text( + axis.text( 1.1, 0.75, textstr_j, - transform=prof.transAxes, + transform=axis.transAxes, fontsize=9, verticalalignment="top", bbox=props_j, ) - prof.text( + axis.text( 0.05, 0.04, "*Current profile is assumed to be parabolic", @@ -3747,7 +3755,7 @@ def plot_jprofile(prof): ha="left", transform=plt.gcf().transFigure, ) - prof.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2) + axis.grid(True, which="both", linestyle="--", linewidth=0.5, alpha=0.2) def plot_t_profiles(prof, demo_ranges, mfile_data, scan): @@ -10954,7 +10962,7 @@ def main_plot( # Plot current density profile ax12 = fig4.add_subplot(4, 3, 10) ax12.set_position([0.075, 0.125, 0.25, 0.15]) - plot_jprofile(ax12) + plot_jprofile(ax12, mfile_data=m_file_data, scan=scan) # Plot q profile ax13 = fig4.add_subplot(4, 3, 12) diff --git a/process/physics.py b/process/physics.py index 92d5fd3726..5450964a45 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4333,6 +4333,13 @@ def outplas(self): physics_variables.j_plasma_on_axis, "OP ", ) + po.ovarrf( + self.outfile, + "On-axis circular plasma current density (A/m2)", + "(j_plasma_circular_on_axis)", + physics_variables.j_plasma_circular_on_axis, + "OP ", + ) po.ovarrf( self.outfile, "Vertical field at plasma (T)", diff --git a/process/plasma_profiles.py b/process/plasma_profiles.py index ddd6a798e9..aded19f186 100644 --- a/process/plasma_profiles.py +++ b/process/plasma_profiles.py @@ -336,6 +336,17 @@ def calculate_profile_factors(self) -> None: ) ) + physics_variables.j_plasma_circular_on_axis = ( + (physics_variables.c_plasma_circular) + * 2 + / ( + sp.special.beta(0.5, physics_variables.alphaj + 1) + * 2 + * np.pi + * physics_variables.rminor + ) + ) + @staticmethod def calculate_parabolic_profile_factors() -> None: """ From d294a9c5c7f169868d2cf78bb457e13ed3e31241 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:34:10 +0100 Subject: [PATCH 24/29] :sparkle: Add j_plasma_circular_on_axis variable for central plasma current density in circular plasma and initialize in init_physics_variables function --- 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 3bea28006f..d3aad2d49e 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -840,6 +840,9 @@ j_plasma_on_axis: float = None """Central plasma current density (A/m2)""" +j_plasma_circular_on_axis: float = None +"""Central plasma current density for a circular plasma (A/m2)""" + n_plasma_profile_elements: int = None """Number of elements in plasma profile""" @@ -1523,6 +1526,7 @@ def init_physics_variables(): global pres_plasma_ion_total_profile global pres_plasma_fuel_profile global j_plasma_on_axis + global j_plasma_circular_on_axis global n_plasma_profile_elements global f_dd_branching_trit global pden_plasma_alpha_mw @@ -1783,6 +1787,7 @@ def init_physics_variables(): pres_plasma_ion_total_profile = [] pres_plasma_fuel_profile = [] j_plasma_on_axis = 0.0 + j_plasma_circular_on_axis = 0.0 n_plasma_profile_elements = 500 f_dd_branching_trit = 0.0 pden_plasma_alpha_mw = 0.0 From 39df5ffa7c6f1702ec06a9be1d0e07f8898f5de5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:35:19 +0100 Subject: [PATCH 25/29] Add b_plasma_circular_poloidal_profile variable and initialize in init_physics_variables function --- 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 d3aad2d49e..00eb6ef4d1 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -304,6 +304,9 @@ b_plasma_poloidal_profile: list[float] = None """poloidal field profile in plasma (T)""" +b_plasma_circular_poloidal_profile: list[float] = None +"""poloidal field profile in circular plasma (T)""" + b_plasma_total_profile: list[float] = None """total field profile in plasma (T)""" @@ -1422,6 +1425,7 @@ def init_physics_variables(): global b_plasma_toroidal_on_axis global b_plasma_toroidal_profile global b_plasma_poloidal_profile + global b_plasma_circular_poloidal_profile global b_plasma_total_profile global b_plasma_total global burnup @@ -1683,6 +1687,7 @@ def init_physics_variables(): b_plasma_toroidal_on_axis = 5.68 b_plasma_toroidal_profile = [] b_plasma_poloidal_profile = [] + b_plasma_circular_poloidal_profile = [] b_plasma_total_profile = [] b_plasma_total = 0.0 burnup = 0.0 From eca507febfcfaaa36ca2c1c4b17392bd40efb117 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 18:58:41 +0100 Subject: [PATCH 26/29] :sparkle: Add b_plasma_circular_poloidal_profile calculation and update plotting function for circular plasma visualization --- process/data_structure/physics_variables.py | 5 +++++ process/io/plot_proc.py | 14 ++++++++++++++ process/physics.py | 16 ++++++++++++++++ 3 files changed, 35 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 00eb6ef4d1..2a811aa93d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1125,6 +1125,9 @@ qstar: float = None """cylindrical safety factor""" +q_95_circular: float = None +"""Safety factor at 95% flux surface for circular plasma""" + rad_fraction_sol: float = None """SoL radiation fraction""" @@ -1594,6 +1597,7 @@ def init_physics_variables(): global tauratio global q95_min global qstar + global q_95_circular global rad_fraction_sol global rad_fraction_total global f_nd_alpha_electron @@ -1856,6 +1860,7 @@ def init_physics_variables(): tauratio = 1.0 q95_min = 0.0 qstar = 0.0 + q_95_circular = 0.0 rad_fraction_sol = 0.8 rad_fraction_total = 0.0 f_nd_alpha_electron = 0.1 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index cc8a091261..16d27e7902 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -10721,6 +10721,11 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): for i in range(2 * n_plasma_profile_elements) ] + b_plasma_circular_poloidal_profile = [ + mfile_data.data[f"b_plasma_circular_poloidal_profile{i}"].get_scan(scan) + for i in range(2 * n_plasma_profile_elements) + ] + b_plasma_total_profile = [ mfile_data.data[f"b_plasma_total_profile{i}"].get_scan(scan) for i in range(2 * n_plasma_profile_elements) @@ -10745,6 +10750,15 @@ def plot_magnetic_fields_in_plasma(axis, mfile_data, scan): label="Poloidal B-field [T]", ) + axis.plot( + np.linspace( + rmajor - rminor, rmajor + rminor, len(b_plasma_circular_poloidal_profile) + ), + b_plasma_circular_poloidal_profile, + color="orange", + label="Circular Poloidal B-field [T]", + ) + axis.plot( np.linspace(rmajor - rminor, rmajor + rminor, len(b_plasma_total_profile)), b_plasma_total_profile, diff --git a/process/physics.py b/process/physics.py index 5450964a45..f7213e8287 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1762,6 +1762,15 @@ def physics(self): physics_variables.b_plasma_poloidal_profile ) + physics_variables.b_plasma_circular_poloidal_profile = ( + self.calculate_poloidal_field_profile( + j_plasma_on_axis=physics_variables.j_plasma_circular_on_axis, + alphaj=physics_variables.alphaj, + n_profile_elements=physics_variables.n_plasma_profile_elements * 2, + rminor=physics_variables.rminor, + ) + ) + # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components physics_variables.b_plasma_total_profile = np.sqrt( np.square(physics_variables.b_plasma_toroidal_profile) @@ -4368,6 +4377,13 @@ def outplas(self): f"b_plasma_poloidal_profile{i}", physics_variables.b_plasma_poloidal_profile[i], ) + for i in range(len(physics_variables.b_plasma_circular_poloidal_profile)): + po.ovarre( + self.mfile, + f"Poloidal field in circular plasma at point {i}", + f"b_plasma_circular_poloidal_profile{i}", + physics_variables.b_plasma_circular_poloidal_profile[i], + ) for i in range(len(physics_variables.b_plasma_total_profile)): po.ovarre( self.mfile, From 861f496036d1eeb162e61d688ee21fe0e4690268 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 19:25:40 +0100 Subject: [PATCH 27/29] :sparkle: Add q_circular_profile variable and update init_physics_variables function; calculate q_95_circular in Physics class --- process/data_structure/physics_variables.py | 4 ++++ process/io/plot_proc.py | 5 +++-- process/physics.py | 14 ++++++++++++++ 3 files changed, 21 insertions(+), 2 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 2a811aa93d..a2ab4440a6 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1128,6 +1128,8 @@ q_95_circular: float = None """Safety factor at 95% flux surface for circular plasma""" +q_circular_profile: list[float] = None +"""safety factor profile for circular plasma""" rad_fraction_sol: float = None """SoL radiation fraction""" @@ -1598,6 +1600,7 @@ def init_physics_variables(): global q95_min global qstar global q_95_circular + global q_circular_profile global rad_fraction_sol global rad_fraction_total global f_nd_alpha_electron @@ -1861,6 +1864,7 @@ def init_physics_variables(): q95_min = 0.0 qstar = 0.0 q_95_circular = 0.0 + q_circular_profile = [] rad_fraction_sol = 0.8 rad_fraction_total = 0.0 f_nd_alpha_electron = 0.1 diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 16d27e7902..b0ab9c96b0 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -3897,12 +3897,13 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan): textstr_q = "\n".join(( r"$q_0$: " + f"{q0:.3f}\n", r"$q_{95}$: " + f"{q95:.3f}\n", - r"$q_{\text{cyl}}$: " + f"{mfile_data.data['qstar'].get_scan(scan):.3f}", + r"$q_{\text{cyl}}$: " + f"{mfile_data.data['qstar'].get_scan(scan):.3f}\n", + r"$q_{95,\text{circular}}$: " + f"{mfile_data.data['q_95_circular'].get_scan(scan):.3f}", )) props_q = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5} prof.text( - -0.4, + -0.5, 0.75, textstr_q, transform=prof.transAxes, diff --git a/process/physics.py b/process/physics.py index f7213e8287..aab4cf0754 100644 --- a/process/physics.py +++ b/process/physics.py @@ -3882,6 +3882,14 @@ def calculate_plasma_current( * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) ) + physics_variables.q_95_circular = ( + 2 + * np.pi + * rminor**2 + * b_plasma_toroidal_on_axis + / (constants.RMU0 * rmajor * plasma_current) + ) + # Normalised beta from Troyon beta limit physics_variables.beta_norm_total = ( 1.0e8 @@ -4434,6 +4442,12 @@ def outplas(self): physics_variables.qstar, "OP ", ) + po.ovarrf( + self.outfile, + "Edge safety factor for a circular plasma (q_circular)", + "(q_95_circular)", + physics_variables.q_95_circular, + ) if physics_variables.i_plasma_geometry == 1: po.ovarrf( From c0170390502aeb1c7dc2687cb1d3038d1c03bf2a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 19:49:16 +0100 Subject: [PATCH 28/29] :sparkle: Implement calculation and plotting for circular plasma safety factor profile --- process/io/plot_proc.py | 25 ++++++++++++++++++++++--- process/physics.py | 24 ++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 3 deletions(-) diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index b0ab9c96b0..920fd79fd1 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -3868,9 +3868,27 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan): q_r_nevin = q0 + (q95 - q0) * (rho + rho * rho + rho**3) / (3.0) q_r_sauter = q0 + (q95 - q0) * (rho * rho) + n_plasma_profile_elements = int( + mfile_data.data["n_plasma_profile_elements"].get_scan(scan) + ) + + q_circular_profile = [ + mfile_data.data[f"q_circular_profile{i}"].get_scan(scan) + for i in range(n_plasma_profile_elements) + ] + + prof.plot( + np.linspace(0, 1, n_plasma_profile_elements), + q_circular_profile, + label="Circular", + linestyle="--", + color="black", + ) + prof.plot(rho, q_r_nevin, label="Nevins") prof.plot(rho, q_r_sauter, label="Sauter") prof.legend() + prof.figure.tight_layout() # Ranges # --- @@ -3884,9 +3902,9 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan): prof.set_ylim([0, q95 * 1.2]) prof.text( - 0.6, + 0.45, 0.04, - "*Profile is not calculated, only $q_0$ and $q_{95}$ are known.", + "*Only circular profile is calculated, for others only $q_0$ and $q_{95}$ are known.", fontsize=10, ha="left", transform=plt.gcf().transFigure, @@ -3898,7 +3916,8 @@ def plot_qprofile(prof, demo_ranges, mfile_data, scan): r"$q_0$: " + f"{q0:.3f}\n", r"$q_{95}$: " + f"{q95:.3f}\n", r"$q_{\text{cyl}}$: " + f"{mfile_data.data['qstar'].get_scan(scan):.3f}\n", - r"$q_{95,\text{circular}}$: " + f"{mfile_data.data['q_95_circular'].get_scan(scan):.3f}", + r"$q_{95,\text{circular}}$: " + + f"{mfile_data.data['q_95_circular'].get_scan(scan):.3f}", )) props_q = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5} diff --git a/process/physics.py b/process/physics.py index aab4cf0754..0bd584eb27 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1771,6 +1771,23 @@ def physics(self): ) ) + physics_variables.q_circular_profile = ( + ( + np.linspace( + 0, + physics_variables.rminor, + physics_variables.n_plasma_profile_elements, + ) + / physics_variables.rmajor + ) + * physics_variables.b_plasma_toroidal_profile[ + physics_variables.n_plasma_profile_elements : + ] + / physics_variables.b_plasma_circular_poloidal_profile[ + physics_variables.n_plasma_profile_elements : + ] + ) + # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components physics_variables.b_plasma_total_profile = np.sqrt( np.square(physics_variables.b_plasma_toroidal_profile) @@ -4399,6 +4416,13 @@ def outplas(self): f"b_plasma_total_profile{i}", physics_variables.b_plasma_total_profile[i], ) + for i in range(len(physics_variables.q_circular_profile)): + po.ovarre( + self.mfile, + f"Safety factor for circular plasma at point {i}", + f"q_circular_profile{i}", + physics_variables.q_circular_profile[i], + ) po.ovarrf( self.outfile, "Average poloidal field (T)", From d719ac81bf186b29969c0fecbd34bdc4d7f14282 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 10 Oct 2025 23:03:08 +0100 Subject: [PATCH 29/29] Enhance b_plasma_circular_poloidal_profile calculation and update q_circular_profile formula in Physics class --- process/physics.py | 67 ++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 56 insertions(+), 11 deletions(-) diff --git a/process/physics.py b/process/physics.py index 0bd584eb27..97bac88907 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1771,21 +1771,66 @@ def physics(self): ) ) + physics_variables.b_plasma_circular_poloidal_profile = ( + ( + constants.RMU0 + * physics_variables.j_plasma_circular_on_axis + * physics_variables.rminor**2 + * ( + 1 + - ( + 1 + - ( + np.linspace( + -1, 1, 2 * physics_variables.n_plasma_profile_elements + ) + ) + ** 2 + ) + ** (physics_variables.alphaj + 1) + ) + ) + / (2 * (physics_variables.alphaj + 1)) + * np.linspace(-1, 1, 2 * physics_variables.n_plasma_profile_elements) + * physics_variables.rminor + ) + + # physics_variables.q_circular_profile = ( + # ( + # np.linspace( + # 0, + # physics_variables.rminor, + # physics_variables.n_plasma_profile_elements, + # ) + # / physics_variables.rmajor + # ) + # * physics_variables.b_plasma_toroidal_profile[ + # physics_variables.n_plasma_profile_elements : + # ] + # / physics_variables.b_plasma_circular_poloidal_profile[ + # physics_variables.n_plasma_profile_elements : + # ] + # ) + rho = np.linspace( + 0, + 1.0, + physics_variables.n_plasma_profile_elements, + ) physics_variables.q_circular_profile = ( ( - np.linspace( - 0, - physics_variables.rminor, - physics_variables.n_plasma_profile_elements, + ( + 2 + * (physics_variables.alphaj + 1) + * physics_variables.b_plasma_toroidal_on_axis + ) + / ( + constants.RMU0 + * physics_variables.rmajor + * physics_variables.j_plasma_circular_on_axis ) - / physics_variables.rmajor ) - * physics_variables.b_plasma_toroidal_profile[ - physics_variables.n_plasma_profile_elements : - ] - / physics_variables.b_plasma_circular_poloidal_profile[ - physics_variables.n_plasma_profile_elements : - ] + * rho**2 + / (1 - (1 - rho**2) ** (physics_variables.alphaj + 1)) ) # Calculate the total magnetic field profile at each point by summing the squares of the toroidal and poloidal components