From e731b7f48f2e82a5956e3e1b982fecc58b8ca0e9 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 10:43:25 +0000 Subject: [PATCH 01/36] Add new plasma physics variables and update iteration variable count --- process/core/solver/iteration_variables.py | 12 ++++++++++++ process/data_structure/numerics.py | 2 +- process/data_structure/physics_variables.py | 14 ++++++++++++++ 3 files changed, 27 insertions(+), 1 deletion(-) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index edae1a4092..57496ca7f9 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -291,6 +291,18 @@ class IterationVariable: 176: IterationVariable( "f_st_coil_aspect", data_structure.stellarator_variables, 0.70, 1.30 ), + 177: IterationVariable( + "f_plasma_particles_lcfs_recycled", data_structure.physics_variables, 0.01, 1.0 + ), + 178: IterationVariable( + "eta_plasma_fuelling", data_structure.physics_variables, 0.01, 1.0 + ), + 179: IterationVariable( + "molflow_plasma_fuelling_vv_injected", + data_structure.physics_variables, + 0.0, + 1e24, + ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 6eafc3aa99..1cc1842fce 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -1,6 +1,6 @@ import numpy as np -ipnvars: int = 177 +ipnvars: int = 180 """total number of variables available for iteration""" ipeqns: int = 92 diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index fccd7a8294..b7ad3626f8 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1065,6 +1065,14 @@ molflow_plasma_fuelling_required: float = None """plasma fuelling rate (nucleus-pairs/s)""" +f_plasma_particles_lcfs_recycled: float = None +"""fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)""" + +eta_plasma_fuelling: float = None +"""fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)""" + +molflow_plasma_fuelling_vv_injected: float = None +"""plasma fuelling rate into the vacuum vessel (nucleus-pairs/s)""" tauratio: float = None """tauratio /1.0/ : ratio of He and pellet particle confinement times""" @@ -1617,6 +1625,9 @@ def init_physics_variables(): q0, \ q95, \ molflow_plasma_fuelling_required, \ + f_plasma_particles_lcfs_recycled, \ + eta_plasma_fuelling, \ + molflow_plasma_fuelling_vv_injected, \ tauratio, \ q95_min, \ qstar, \ @@ -1908,6 +1919,9 @@ def init_physics_variables(): q0 = 1.0 q95 = 0.0 molflow_plasma_fuelling_required = 0.0 + f_plasma_particles_lcfs_recycled = 0.0 + eta_plasma_fuelling = 0.0 + molflow_plasma_fuelling_vv_injected = 0.0 tauratio = 1.0 q95_min = 0.0 qstar = 0.0 From 02d6623d0d9260cf35d4b502fc425a0e102beac6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 10:59:44 +0000 Subject: [PATCH 02/36] Add particle balance constraint equation and update constraint count --- process/core/solver/constraints.py | 36 ++++++++++++++++++++++++++++++ process/data_structure/numerics.py | 2 +- 2 files changed, 37 insertions(+), 1 deletion(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 949b4439bf..cdaa11f070 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1855,6 +1855,42 @@ def constraint_equation_92(constraint_registration): ) +@ConstraintManager.register_constraint(93, "particles/s", "=") +def constraint_equation_93(): + """ + Particle balance + """ + # pscaling: total transport power per volume (MW/m3) + # NEED TO ADD PROPER FUSION RATES FOR EACH REACTION + + # Numerator shall be the tritium particle balance + numerator = ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + - data_structure.physics_variables.fusrat_total + ) - ( + ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_tritium + ) + / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + ) + ) + + # Alpha particle balance + denominator = data_structure.physics_variables.fusrat_total - ( + data_structure.physics_variables.nd_plasma_alphas_vol_avg + * data_structure.physics_variables.vol_plasma + ) / (data_structure.physics_variables.t_alpha_confinement) + + cc = 1.0 - numerator / numerator + + return ConstraintResult(cc, denominator * (1.0 - cc), denominator * cc) + + def constraint_eqns(m: int, ieqn: int): """Evaluates the constraints given the current state of PROCESS. diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 1cc1842fce..35ac0a04e0 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -3,7 +3,7 @@ ipnvars: int = 180 """total number of variables available for iteration""" -ipeqns: int = 92 +ipeqns: int = 93 """number of constraint equations available""" ipnfoms: int = 19 From 1678fc17e0139596716bf4fdafe6b88d3114d164 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 11:03:21 +0000 Subject: [PATCH 03/36] Add new plasma physics input variables for recycling and fuelling --- process/core/input.py | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/process/core/input.py b/process/core/input.py index 7cb7b76577..405d55ace4 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1208,6 +1208,15 @@ def __post_init__(self): "molflow_vac_pumps": InputVariable( data_structure.vacuum_variables, float, range=(0.0, 1e30) ), + "f_plasma_particles_lcfs_recycled": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "eta_plasma_fuelling": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "molflow_plasma_fuelling_vv_injected": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1e24) + ), "pflux_plant_floor_electric": InputVariable( data_structure.heat_transport_variables, float, range=(0.0, 1000.0) ), From 762276e70f22eaa5235cca0276662dc6db60ff3f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 11:12:39 +0000 Subject: [PATCH 04/36] Refactor particle balance constraint equation and update iteration variable for consistency --- process/core/solver/constraints.py | 6 ++---- process/core/solver/iteration_variables.py | 2 +- process/data_structure/numerics.py | 3 +++ 3 files changed, 6 insertions(+), 5 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index cdaa11f070..3044f3e99a 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1856,7 +1856,7 @@ def constraint_equation_92(constraint_registration): @ConstraintManager.register_constraint(93, "particles/s", "=") -def constraint_equation_93(): +def constraint_equation_93(constraint_registration): """ Particle balance """ @@ -1886,9 +1886,7 @@ def constraint_equation_93(): * data_structure.physics_variables.vol_plasma ) / (data_structure.physics_variables.t_alpha_confinement) - cc = 1.0 - numerator / numerator - - return ConstraintResult(cc, denominator * (1.0 - cc), denominator * cc) + return eq(numerator, denominator, constraint_registration) def constraint_eqns(m: int, ieqn: int): diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 57496ca7f9..34eea42820 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -300,7 +300,7 @@ class IterationVariable: 179: IterationVariable( "molflow_plasma_fuelling_vv_injected", data_structure.physics_variables, - 0.0, + 1.0, 1e24, ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 35ac0a04e0..132943b939 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -188,6 +188,7 @@ * (90) Lower Limit on number of stress load cycles for CS * (91) Checking if the design point is ECRH ignitable * (92) D/T/He3 ratio in fuel sums to 1 +* (93) Particle balance consistency constraint """ ixc: list[int] = None @@ -314,6 +315,7 @@ * (115) NOT USED * (116) NOT USED * (117) NOT USED +* (118) NOT USED * (119) temp_plasma_separatrix_kev: separatrix temperature calculated by the Kallenbach divertor model * (120) ttarget: Plasma temperature adjacent to divertor sheath [eV] * (121) neratio: ratio of mean SOL density at OMP to separatrix density at OMP @@ -621,6 +623,7 @@ def init_numerics(): "CS stress load cycles ", "ECRH ignitability ", "Fuel composition consistency ", + "Particle balance consistency ", ] ixc = np.array([0] * ipnvars) From 60932ab7d29cc58d1bcdedbbe393786173749dc4 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 14:56:21 +0000 Subject: [PATCH 05/36] Update particle balance constraints and increment equation count --- process/core/solver/constraints.py | 42 +++++++++++++++++++++++++++--- process/data_structure/numerics.py | 6 +++-- 2 files changed, 42 insertions(+), 6 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 3044f3e99a..2bcd0f2d38 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1855,10 +1855,10 @@ def constraint_equation_92(constraint_registration): ) -@ConstraintManager.register_constraint(93, "particles/s", "=") +@ConstraintManager.register_constraint(93, "particles/s", "<=") def constraint_equation_93(constraint_registration): """ - Particle balance + Tritium particle balance """ # pscaling: total transport power per volume (MW/m3) # NEED TO ADD PROPER FUSION RATES FOR EACH REACTION @@ -1880,13 +1880,47 @@ def constraint_equation_93(constraint_registration): ) ) + return leq(numerator, 1e10, constraint_registration) + + +@ConstraintManager.register_constraint(94, "particles/s", "<=") +def constraint_equation_94(constraint_registration): + """ + Deuterium particle balance + """ + + numerator = ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + - data_structure.physics_variables.fusrat_total + ) - ( + ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_deuterium + ) + / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + ) + ) + + return leq(numerator, 1e10, constraint_registration) + + +@ConstraintManager.register_constraint(95, "particles/s", "<=") +def constraint_equation_95(constraint_registration): + """ + Alpha particle balance + """ + # Alpha particle balance - denominator = data_structure.physics_variables.fusrat_total - ( + numerator = data_structure.physics_variables.fusrat_total - ( data_structure.physics_variables.nd_plasma_alphas_vol_avg * data_structure.physics_variables.vol_plasma ) / (data_structure.physics_variables.t_alpha_confinement) - return eq(numerator, denominator, constraint_registration) + return leq(numerator, 1e10, constraint_registration) def constraint_eqns(m: int, ieqn: int): diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 132943b939..18fa48ee97 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -3,7 +3,7 @@ ipnvars: int = 180 """total number of variables available for iteration""" -ipeqns: int = 93 +ipeqns: int = 95 """number of constraint equations available""" ipnfoms: int = 19 @@ -623,7 +623,9 @@ def init_numerics(): "CS stress load cycles ", "ECRH ignitability ", "Fuel composition consistency ", - "Particle balance consistency ", + "Tritium particle balance consistency ", + "Deuterium particle balance consistency ", + "Alpha particle balance consistency ", ] ixc = np.array([0] * ipnvars) From 4245b7947860f86681401823ad709b94cf58ac57 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 15:41:44 +0000 Subject: [PATCH 06/36] Update constraint limits --- process/core/solver/constraints.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 2bcd0f2d38..fc37afe593 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1855,7 +1855,7 @@ def constraint_equation_92(constraint_registration): ) -@ConstraintManager.register_constraint(93, "particles/s", "<=") +@ConstraintManager.register_constraint(93, "particles/s", "=") def constraint_equation_93(constraint_registration): """ Tritium particle balance @@ -1880,10 +1880,10 @@ def constraint_equation_93(constraint_registration): ) ) - return leq(numerator, 1e10, constraint_registration) + return eq(numerator / 1e20, 1e-8, constraint_registration) -@ConstraintManager.register_constraint(94, "particles/s", "<=") +@ConstraintManager.register_constraint(94, "particles/s", "=") def constraint_equation_94(constraint_registration): """ Deuterium particle balance @@ -1905,10 +1905,10 @@ def constraint_equation_94(constraint_registration): ) ) - return leq(numerator, 1e10, constraint_registration) + return eq(numerator / 1e20, 1e-8, constraint_registration) -@ConstraintManager.register_constraint(95, "particles/s", "<=") +@ConstraintManager.register_constraint(95, "particles/s", "=") def constraint_equation_95(constraint_registration): """ Alpha particle balance @@ -1918,9 +1918,12 @@ def constraint_equation_95(constraint_registration): numerator = data_structure.physics_variables.fusrat_total - ( data_structure.physics_variables.nd_plasma_alphas_vol_avg * data_structure.physics_variables.vol_plasma - ) / (data_structure.physics_variables.t_alpha_confinement) + ) / ( + (data_structure.physics_variables.t_energy_confinement * 4.0) + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + ) - return leq(numerator, 1e10, constraint_registration) + return eq(numerator / 1e20, 1e-8, constraint_registration) def constraint_eqns(m: int, ieqn: int): From 5cd3fd518d55a4220139d773e0cc342e826a79f5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 16:49:07 +0000 Subject: [PATCH 07/36] Add new total D-T fusion rate variables --- process/core/io/plot_proc.py | 10 ++++++---- process/core/solver/constraints.py | 6 +++--- process/data_structure/physics_variables.py | 10 ++++++++++ process/models/physics/physics.py | 22 +++++++++++++++++++++ 4 files changed, 41 insertions(+), 7 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 9165eb4840..a985a0b064 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -11322,9 +11322,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ============================================================================ textstr_dt = ( - f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.2f} MW \n" - f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.2f} MW\n" + f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.4f} MW\n" + f"Total fusion rate: {mfile.get('fusrat_dt_total', scan=scan):.4e} reactions/s\n" + f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.4f} MW \n" + f"Plasma fusion rate: {mfile.get('fusrat_plasma_dt', scan=scan):.4e} reactions/s\n" + f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.4f} MW" ) axis.text( @@ -11344,7 +11346,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.24, + 0.275, 0.8, "$\\text{D - T}$", fontsize=20, diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index fc37afe593..9dbdae6887 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1867,7 +1867,7 @@ def constraint_equation_93(constraint_registration): numerator = ( data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - - data_structure.physics_variables.fusrat_total + - data_structure.physics_variables.fusrat_dt_total ) - ( ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg @@ -1892,7 +1892,7 @@ def constraint_equation_94(constraint_registration): numerator = ( data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - - data_structure.physics_variables.fusrat_total + - data_structure.physics_variables.fusrat_dt_total ) - ( ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg @@ -1915,7 +1915,7 @@ def constraint_equation_95(constraint_registration): """ # Alpha particle balance - numerator = data_structure.physics_variables.fusrat_total - ( + numerator = data_structure.physics_variables.fusrat_dt_total - ( data_structure.physics_variables.nd_plasma_alphas_vol_avg * data_structure.physics_variables.vol_plasma ) / ( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index b7ad3626f8..35e46e296b 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -458,6 +458,12 @@ fusrat_total: float = None """fusion reaction rate, from beams and plasma (reactions/sec)""" +fusrat_plasma_dt: float = None +""" D-T fusion reaction rate in plasma, (reactions/sec)""" + +fusrat_dt_total: float = None +""" Total D-T fusion reaction rate from beams and plasma, (reactions/sec)""" + fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -1502,6 +1508,8 @@ def init_physics_variables(): f_plasma_fuel_tritium, \ fusden_total, \ fusrat_total, \ + fusrat_plasma_dt, \ + fusrat_dt_total, \ fusrat_plasma_dt_profile, \ fusrat_plasma_dd_triton_profile, \ fusrat_plasma_dd_helion_profile, \ @@ -1794,6 +1802,8 @@ def init_physics_variables(): f_plasma_fuel_tritium = 0.5 fusden_total = 0.0 fusrat_total = 0.0 + fusrat_plasma_dt = 0.0 + fusrat_dt_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] fusrat_plasma_dd_helion_profile = [] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 34f987e425..22f4564b7a 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -655,6 +655,7 @@ def run(self): + (1.0 / (1.0 - constants.DT_NEUTRON_ENERGY_FRACTION)) * physics_variables.p_beam_alpha_mw ) + physics_variables.p_beam_neutron_mw = physics_variables.p_beam_alpha_mw * ( constants.DT_NEUTRON_ENERGY_FRACTION / (1 - constants.DT_NEUTRON_ENERGY_FRACTION) @@ -671,6 +672,13 @@ def run(self): physics_variables.fusrat_total = ( physics_variables.fusden_total * physics_variables.vol_plasma ) + physics_variables.fusrat_plasma_dt = (physics_variables.p_plasma_dt_mw * 1e6) / ( + constants.D_T_ENERGY + ) + + physics_variables.fusrat_dt_total = ( + physics_variables.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) + ) # Create some derived values and add beam contribution to fusion power ( @@ -2641,6 +2649,20 @@ def outplas(self): physics_variables.fusrat_total, "OP ", ) + po.ovarre( + self.outfile, + "D-T Fusion rate: total (reactions/sec)", + "(fusrat_dt_total)", + physics_variables.fusrat_dt_total, + "OP ", + ) + po.ovarre( + self.outfile, + "D-T Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dt)", + physics_variables.fusrat_plasma_dt, + "OP ", + ) po.ovarre( self.outfile, "Fusion rate density: total (reactions/m3/sec)", From ac439129fb978ad28dd26baf5312ec6a63708f8f Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 18 Mar 2026 16:30:05 +0000 Subject: [PATCH 08/36] Update consistency equations for particle consistency --- process/core/scan.py | 6 ++++ process/core/solver/constraints.py | 45 ++++++++++++++---------------- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/process/core/scan.py b/process/core/scan.py index b2a7035ad1..0ac7c495f7 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -615,6 +615,12 @@ def post_optimise(self, ifail: int): f"(ineq_con{numerics.icc[i]:03d})", -constraint.normalised_residual, ) + process_output.ovarre( + constants.MFILE, + f"{name} error", + f"(ineq_err{numerics.icc[i]:03d})", + -constraint.residual, + ) process_output.ovarre( constants.MFILE, f"{name} physical value", diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 9dbdae6887..3f48581128 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1868,19 +1868,17 @@ def constraint_equation_93(constraint_registration): data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - data_structure.physics_variables.fusrat_dt_total - ) - ( - ( - data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg - * data_structure.physics_variables.vol_plasma - * data_structure.physics_variables.f_plasma_fuel_tritium - ) - / ( - data_structure.physics_variables.t_energy_confinement - / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) - ) + ) + denominator = ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_tritium + ) / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) ) - return eq(numerator / 1e20, 1e-8, constraint_registration) + return eq(numerator, denominator, constraint_registration) @ConstraintManager.register_constraint(94, "particles/s", "=") @@ -1893,19 +1891,17 @@ def constraint_equation_94(constraint_registration): data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - data_structure.physics_variables.fusrat_dt_total - ) - ( - ( - data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg - * data_structure.physics_variables.vol_plasma - * data_structure.physics_variables.f_plasma_fuel_deuterium - ) - / ( - data_structure.physics_variables.t_energy_confinement - / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) - ) + ) + denominator = ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_deuterium + ) / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) ) - return eq(numerator / 1e20, 1e-8, constraint_registration) + return eq(numerator, denominator, constraint_registration) @ConstraintManager.register_constraint(95, "particles/s", "=") @@ -1915,7 +1911,8 @@ def constraint_equation_95(constraint_registration): """ # Alpha particle balance - numerator = data_structure.physics_variables.fusrat_dt_total - ( + numerator = data_structure.physics_variables.fusrat_dt_total + denominator = ( data_structure.physics_variables.nd_plasma_alphas_vol_avg * data_structure.physics_variables.vol_plasma ) / ( @@ -1923,7 +1920,7 @@ def constraint_equation_95(constraint_registration): / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) ) - return eq(numerator / 1e20, 1e-8, constraint_registration) + return eq(numerator, denominator, constraint_registration) def constraint_eqns(m: int, ieqn: int): From 5fb3687984f31ea60ed6f9a144777bb358b44d7e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:32:06 +0000 Subject: [PATCH 09/36] Update molflow plasma fuelling variable ranges for consistency --- process/core/input.py | 2 +- process/core/solver/iteration_variables.py | 2 +- process/data_structure/physics_variables.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/process/core/input.py b/process/core/input.py index 405d55ace4..2b55e7b3a3 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1215,7 +1215,7 @@ def __post_init__(self): data_structure.physics_variables, float, range=(0.0, 1.0) ), "molflow_plasma_fuelling_vv_injected": InputVariable( - data_structure.physics_variables, float, range=(0.0, 1e24) + data_structure.physics_variables, float, range=(1e18, 1e24) ), "pflux_plant_floor_electric": InputVariable( data_structure.heat_transport_variables, float, range=(0.0, 1000.0) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 34eea42820..efbc96a559 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -300,7 +300,7 @@ class IterationVariable: 179: IterationVariable( "molflow_plasma_fuelling_vv_injected", data_structure.physics_variables, - 1.0, + 1e18, 1e24, ), } diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 35e46e296b..ed3c1190d0 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1931,7 +1931,7 @@ def init_physics_variables(): molflow_plasma_fuelling_required = 0.0 f_plasma_particles_lcfs_recycled = 0.0 eta_plasma_fuelling = 0.0 - molflow_plasma_fuelling_vv_injected = 0.0 + molflow_plasma_fuelling_vv_injected = 1e20 tauratio = 1.0 q95_min = 0.0 qstar = 0.0 From 25c5c6836db5713b6074bb823267c3afb4743b08 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:39:04 +0000 Subject: [PATCH 10/36] Refine fusion power output formatting for improved precision in plot displays --- process/core/io/plot_proc.py | 40 ++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index a985a0b064..1170f4112c 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -2978,7 +2978,7 @@ def plot_main_plasma_information( f"| D: {mfile.get('f_plasma_fuel_deuterium', scan=scan):.2f} | T: {mfile.get('f_plasma_fuel_tritium', scan=scan):.2f} | 3He: {mfile.get('f_plasma_fuel_helium3', scan=scan):.2f} |\n\n" f"Fusion Power, $P_{{\\text{{fus}}}}:$ {mfile.get('p_fusion_total_mw', scan=scan):,.2f} MW\n" f"D-T Power, $P_{{\\text{{fus,DT}}}}:$ {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" + f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" f"D-3He Power, $P_{{\\text{{fus,D3He}}}}:$ {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW\n" f"Alpha Power, $P_{{\\alpha}}:$ {mfile.get('p_alpha_total_mw', scan=scan):,.2f} MW" ) @@ -11299,8 +11299,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # Add plasma volume, areas and shaping information textstr_general = ( f"Total fusion rate: {mfile.get('fusrat_total', scan=scan):.4e} reactions/s\n" - f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m3/s\n" - f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m3/s\n" + f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m$^3$/s\n" + f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m$^3$/s" ) axis.text( @@ -11346,7 +11346,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.275, + 0.285, 0.8, "$\\text{D - T}$", fontsize=20, @@ -11357,7 +11357,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ================================================= textstr_dd = ( - f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" + f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" ) @@ -11388,7 +11388,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ================================================= - textstr_dhe3 = f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW \n\n" + textstr_dhe3 = f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.4f} MW \n\n" axis.text( 0.05, @@ -11418,15 +11418,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ================================================= textstr_alpha = ( - f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.2f} MW\n" - f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.2f} MW\n" - f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.2f} MW\n\n" - f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m3/sec\n" - f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m3/sec\n\n" - f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n" - f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n" - f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n" - f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n" + f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.4f} MW\n" + f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.4f} MW\n" + f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.4f} MW\n\n" + f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m$^3$/sec\n" + f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m$^3$/sec\n\n" + f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m$^3$\n" + f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m$^3$\n\n" + f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m$^3$\n" + f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m$^3$" ) axis.text( @@ -11457,11 +11457,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ================================================= textstr_neutron = ( - f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.2f} MW\n" - f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.2f} MW\n" - f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.2f} MW\n\n" - f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m3\n" - f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m3\n" + f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.4f} MW\n" + f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.4f} MW\n" + f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.4f} MW\n\n" + f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m$^3$\n" + f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m$^3$" ) axis.text( From e5ab85b1898090fe827a9729f0757224454c0bfb Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:42:10 +0000 Subject: [PATCH 11/36] Add D-D fusion reaction rate variables for helium and tritium branches --- process/data_structure/physics_variables.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ed3c1190d0..b0256d6eb1 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -464,6 +464,12 @@ fusrat_dt_total: float = None """ Total D-T fusion reaction rate from beams and plasma, (reactions/sec)""" +fusrat_plasma_dd_helion: float = None +"""D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" + +fusrat_plasma_dd_triton: float = None +"""D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)""" + fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -1510,6 +1516,8 @@ def init_physics_variables(): fusrat_total, \ fusrat_plasma_dt, \ fusrat_dt_total, \ + fusrat_plasma_dd_helion, \ + fusrat_plasma_dd_triton, \ fusrat_plasma_dt_profile, \ fusrat_plasma_dd_triton_profile, \ fusrat_plasma_dd_helion_profile, \ @@ -1803,6 +1811,8 @@ def init_physics_variables(): fusden_total = 0.0 fusrat_total = 0.0 fusrat_plasma_dt = 0.0 + fusrat_plasma_dd_helion = 0.0 + fusrat_plasma_dd_triton = 0.0 fusrat_dt_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] From 873d1ebd333e8679b8cd074134274516b0c90b9b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:08:30 +0000 Subject: [PATCH 12/36] Add D-D fusion rate calculations for tritium and helium branches --- process/core/io/plot_proc.py | 8 +++++--- process/models/physics/fusion_reactions.py | 8 ++++++++ process/models/physics/physics.py | 15 +++++++++++++++ 3 files changed, 28 insertions(+), 3 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 1170f4112c..fdbce5477f 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -11358,7 +11358,9 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): textstr_dd = ( f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" - f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" + f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" + f"D+D -> T fusion rate: {mfile.get('fusrat_plasma_dd_triton', scan=scan):.4e} reactions/s \n" + f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s " ) axis.text( @@ -11378,8 +11380,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.22, - 0.685, + 0.31, + 0.69, "$\\text{D - D}$", fontsize=20, verticalalignment="top", diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 21c111dc23..7a5513c23e 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -427,6 +427,10 @@ def dd_helion_reaction(self): alpha_rate_density = 0.0 proton_rate_density = 0.0 + physics_variables.fusrat_plasma_dd_helion = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-D power density self.dd_power_density += fusion_power_density @@ -521,6 +525,10 @@ def dd_triton_reaction(self): alpha_rate_density = 0.0 proton_rate_density = fusion_rate_density # Proton production rate [m^3/second] + physics_variables.fusrat_plasma_dd_triton = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-D power density self.dd_power_density += fusion_power_density diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 22f4564b7a..00607c24ab 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2663,6 +2663,21 @@ def outplas(self): physics_variables.fusrat_plasma_dt, "OP ", ) + + po.ovarre( + self.outfile, + "D-D -> 3He Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dd_helion)", + physics_variables.fusrat_plasma_dd_helion, + "OP ", + ) + po.ovarre( + self.outfile, + "D-D -> T Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dd_triton)", + physics_variables.fusrat_plasma_dd_triton, + "OP ", + ) po.ovarre( self.outfile, "Fusion rate density: total (reactions/m3/sec)", From 9ef2fc9df0ead6c00db6d40d563c93b9007878f5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:53:07 +0000 Subject: [PATCH 13/36] Add plasma tritium flow rate calculations and contour plotting --- process/core/io/plot_proc.py | 5 ++ process/models/physics/exhaust.py | 96 +++++++++++++++++++++++++++++++ process/models/physics/physics.py | 7 +++ 3 files changed, 108 insertions(+) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index fdbce5477f..eb8cb9697b 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -67,6 +67,7 @@ ElectronBernstein, ElectronCyclotron, ) +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.impurity_radiation import read_impurity_file from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel @@ -13864,6 +13865,10 @@ def main_plot( ax24.set_position([0.08, 0.35, 0.84, 0.57]) plot_system_power_profiles_over_time(ax24, m_file, scan, figs[33]) + + ax25 = figs[31].add_subplot(221) + PlasmaExhaust().plot_tritium_flow_contour(ax25, m_file, scan) + def create_thickness_builds(m_file, scan: int): # Build the dictionaries of radial and vertical build values and cumulative values diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index eb0021fb8e..159703ae01 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -1,5 +1,9 @@ import logging +import matplotlib.pyplot as plt +import numpy as np + +import process.core.io.mfile as mf from process.core import constants logger = logging.getLogger(__name__) @@ -118,3 +122,95 @@ def calculate_eu_demo_re_attachment_metric( return (p_plasma_separatrix_mw * b_plasma_toroidal_on_axis) / ( q95 * aspect * rmajor ) + + @staticmethod + def calculate_plasma_tritium_flow_rate( + eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected, + fusrat_dt_total, + fusrat_plasma_dd, + t_energy_confinement, + f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg, + vol_plasma, + f_plasma_fuel_tritium, + ): + """Calculate the tritium flow rate in the plasma exhaust.""" + + return ( + eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected + - fusrat_dt_total + + fusrat_plasma_dd + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_deuterium_flow_rate( + eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected, + fusrat_dt_total, + fusrat_plasma_dd, + t_energy_confinement, + f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg, + vol_plasma, + f_plasma_fuel_tritium, + ): + """Calculate the tritium flow rate in the plasma exhaust.""" + + return ( + eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected + - fusrat_dt_total + + fusrat_plasma_dd + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of tritium flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + tritium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + tritium_flow[i, j] = self.calculate_plasma_tritium_flow_rate( + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dd=mfile.get("fusrat_plasma_dd_helion", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_tritium=mfile.get("f_plasma_fuel_tritium", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, recycling_range, tritium_flow, levels=15, cmap="RdBu_r" + ) + axis.contour( + fuelling_range, + recycling_range, + tritium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Tritium Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 00607c24ab..18361532d1 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -3293,6 +3293,13 @@ def outplas(self): physics_variables.molflow_plasma_fuelling_required, "OP ", ) + po.ovarre( + self.outfile, + "Fuelling rate (nucleus-pairs/s)", + "(molflow_plasma_fuelling_vv_injected)", + physics_variables.molflow_plasma_fuelling_vv_injected, + "OP ", + ) po.ovarre( self.outfile, "Fuel burn-up rate (reactions/s)", From f0229670530cc7f765290e6ff03cc4f296baa26e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:57:34 +0000 Subject: [PATCH 14/36] Add total D-D fusion rate calculation and update output formatting --- process/core/io/plot_proc.py | 3 ++- process/data_structure/physics_variables.py | 5 +++++ process/models/physics/physics.py | 11 +++++++++++ 3 files changed, 18 insertions(+), 1 deletion(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index eb8cb9697b..f6a2debbb3 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -11361,7 +11361,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" f"D+D -> T fusion rate: {mfile.get('fusrat_plasma_dd_triton', scan=scan):.4e} reactions/s \n" - f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s " + f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s \n" + f"Total D-D fusion rate: {mfile.get('fusrat_plasma_dd_total', scan=scan):.4e} reactions/s" ) axis.text( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index b0256d6eb1..7b3b27c1be 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -470,6 +470,9 @@ fusrat_plasma_dd_triton: float = None """D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)""" +fusrat_plasma_dd_total: float = None +"""Total D-D fusion reaction rate in plasma, (reactions/sec)""" + fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -1518,6 +1521,7 @@ def init_physics_variables(): fusrat_dt_total, \ fusrat_plasma_dd_helion, \ fusrat_plasma_dd_triton, \ + fusrat_plasma_dd_total, \ fusrat_plasma_dt_profile, \ fusrat_plasma_dd_triton_profile, \ fusrat_plasma_dd_helion_profile, \ @@ -1813,6 +1817,7 @@ def init_physics_variables(): fusrat_plasma_dt = 0.0 fusrat_plasma_dd_helion = 0.0 fusrat_plasma_dd_triton = 0.0 + fusrat_plasma_dd_total = 0.0 fusrat_dt_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 18361532d1..f231faacb1 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -675,6 +675,10 @@ def run(self): physics_variables.fusrat_plasma_dt = (physics_variables.p_plasma_dt_mw * 1e6) / ( constants.D_T_ENERGY ) + physics_variables.fusrat_plasma_dd_total = ( + physics_variables.fusrat_plasma_dd_helion + + physics_variables.fusrat_plasma_dd_triton + ) physics_variables.fusrat_dt_total = ( physics_variables.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) @@ -2678,6 +2682,13 @@ def outplas(self): physics_variables.fusrat_plasma_dd_triton, "OP ", ) + po.ovarre( + self.outfile, + "D-D Fusion rate: total (reactions/sec)", + "(fusrat_plasma_dd_total)", + physics_variables.fusrat_plasma_dd_total, + "OP ", + ) po.ovarre( self.outfile, "Fusion rate density: total (reactions/m3/sec)", From db568e3f79ca1235b62b4ad8f2991cada5166ba5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 16:35:54 +0000 Subject: [PATCH 15/36] Add deuterium and alpha particle flow rate calculations and contour plotting --- process/core/io/plot_proc.py | 4 + process/models/physics/exhaust.py | 142 +++++++++++++++++++++++++++--- 2 files changed, 134 insertions(+), 12 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index f6a2debbb3..360ac89d3b 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -13869,6 +13869,10 @@ def main_plot( ax25 = figs[31].add_subplot(221) PlasmaExhaust().plot_tritium_flow_contour(ax25, m_file, scan) + ax26 = figs[31].add_subplot(222) + PlasmaExhaust().plot_deuterium_flow_contour(ax26, m_file, scan) + ax27 = figs[31].add_subplot(223) + PlasmaExhaust().plot_alpha_flow_contour(ax27, m_file, scan) def create_thickness_builds(m_file, scan: int): diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index 159703ae01..83c59ae64f 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -128,7 +128,7 @@ def calculate_plasma_tritium_flow_rate( eta_plasma_fuelling, molflow_plasma_fuelling_vv_injected, fusrat_dt_total, - fusrat_plasma_dd, + fusrat_plasma_dd_triton, t_energy_confinement, f_plasma_particles_lcfs_recycled, nd_plasma_fuel_ions_vol_avg, @@ -136,11 +136,11 @@ def calculate_plasma_tritium_flow_rate( f_plasma_fuel_tritium, ): """Calculate the tritium flow rate in the plasma exhaust.""" - + # Assuming 50/50 D-T fuel mix, the tritium fuelling rate is half the total fuelling rate, and the tritium loss includes both D and T contributions from fusion reactions. return ( - eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected + (eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected / 2) - fusrat_dt_total - + fusrat_plasma_dd + + fusrat_plasma_dd_triton - ( (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) @@ -152,25 +152,43 @@ def calculate_plasma_deuterium_flow_rate( eta_plasma_fuelling, molflow_plasma_fuelling_vv_injected, fusrat_dt_total, - fusrat_plasma_dd, + fusrat_plasma_dd_total, t_energy_confinement, f_plasma_particles_lcfs_recycled, nd_plasma_fuel_ions_vol_avg, vol_plasma, - f_plasma_fuel_tritium, + f_plasma_fuel_deuterium, ): - """Calculate the tritium flow rate in the plasma exhaust.""" - + """Calculate the deuterium flow rate in the plasma exhaust.""" + # Assuming 50/50 D-T fuel mix, the deuterium fuelling rate is half the total fuelling rate, and the deuterium loss includes both D and T contributions from fusion reactions. return ( - eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected + (eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected / 2) - fusrat_dt_total - + fusrat_plasma_dd + - 2 * fusrat_plasma_dd_total - ( - (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium) / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) ) ) + @staticmethod + def calculate_plasma_alphas_flow_rate( + fusrat_dt_total, + t_energy_confinement, + f_t_alpha_energy_confinement, + f_plasma_particles_lcfs_recycled, + nd_plasma_alphas_vol_avg, + vol_plasma, + ): + """Calculate the alpha particle flow rate in the plasma exhaust.""" + + # Alpha particle balance + + return fusrat_dt_total - (nd_plasma_alphas_vol_avg * vol_plasma) / ( + (t_energy_confinement * f_t_alpha_energy_confinement) + / (1 - f_plasma_particles_lcfs_recycled) + ) + def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): """Plot contour of tritium flow rate vs recycling and fuelling rate.""" @@ -186,7 +204,9 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): "molflow_plasma_fuelling_vv_injected", scan=scan ), fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), - fusrat_plasma_dd=mfile.get("fusrat_plasma_dd_helion", scan=scan), + fusrat_plasma_dd_triton=mfile.get( + "fusrat_plasma_dd_triton", scan=scan + ), t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), f_plasma_particles_lcfs_recycled=recycling, nd_plasma_fuel_ions_vol_avg=mfile.get( @@ -214,3 +234,101 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): axis.grid(True, which="major", linestyle="-", alpha=0.7) axis.grid(True, which="minor", linestyle=":", alpha=0.4) plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") + + def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of deuterium flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + deuterium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + deuterium_flow[i, j] = self.calculate_plasma_deuterium_flow_rate( + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dd_total=mfile.get( + "fusrat_plasma_dd_total", scan=scan + ), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_deuterium=mfile.get( + "f_plasma_fuel_deuterium", scan=scan + ), + ) + + contour = axis.contourf( + fuelling_range, recycling_range, deuterium_flow, levels=15, cmap="RdBu_r" + ) + axis.contour( + fuelling_range, + recycling_range, + deuterium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Deuterium Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") + + def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + f_t_alpha_energy_confinement_range = np.linspace(4.0, 8.0, 20) + alpha_flow = np.zeros(( + len(recycling_range), + len(f_t_alpha_energy_confinement_range), + )) + + for i, recycling in enumerate(recycling_range): + for j, f_t_alpha_energy_confinement in enumerate( + f_t_alpha_energy_confinement_range + ): + alpha_flow[i, j] = self.calculate_plasma_alphas_flow_rate( + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_alphas_vol_avg=mfile.get( + "nd_plasma_alphas_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_t_alpha_energy_confinement=f_t_alpha_energy_confinement, + ) + + contour = axis.contourf( + f_t_alpha_energy_confinement_range, + recycling_range, + alpha_flow, + levels=15, + cmap="RdBu_r", + ) + axis.contour( + f_t_alpha_energy_confinement_range, + recycling_range, + alpha_flow, + levels=[0], + colors="black", + linewidths=2, + ) + axis.set_xlabel( + "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" + ) + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") From 0949e73755810151259cff7bb9860b9f0d645381 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 17:23:48 +0000 Subject: [PATCH 16/36] Add plotting of fuelling efficiency and recycling fraction in plasma flow contours --- process/models/physics/exhaust.py | 44 ++++++++++++++++++++++++++++++- process/models/physics/physics.py | 14 ++++++++++ 2 files changed, 57 insertions(+), 1 deletion(-) diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index 83c59ae64f..a1004f55d8 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -227,6 +227,20 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): colors="black", linewidths=2, ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") axis.set_ylabel("Recycling Fraction [$R$]") axis.set_title("Plasma Tritium Flow Rate (particles/s)") @@ -275,6 +289,20 @@ def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int colors="black", linewidths=2, ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") axis.set_ylabel("Recycling Fraction [$R$]") axis.set_title("Plasma Deuterium Flow Rate (particles/s)") @@ -287,7 +315,7 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" recycling_range = np.linspace(0.01, 0.99, 20) - f_t_alpha_energy_confinement_range = np.linspace(4.0, 8.0, 20) + f_t_alpha_energy_confinement_range = np.linspace(4.0, 10.0, 20) alpha_flow = np.zeros(( len(recycling_range), len(f_t_alpha_energy_confinement_range), @@ -323,6 +351,20 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): colors="black", linewidths=2, ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan) + axis.plot( + f_t_alpha_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + axis.set_xlabel( "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" ) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index f231faacb1..6e05af8932 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -3311,6 +3311,20 @@ def outplas(self): physics_variables.molflow_plasma_fuelling_vv_injected, "OP ", ) + po.ovarre( + self.outfile, + "Fuelling efficiency", + "(eta_plasma_fuelling)", + physics_variables.eta_plasma_fuelling, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma particles recycled at the LCFS", + "(f_plasma_particles_lcfs_recycled)", + physics_variables.f_plasma_particles_lcfs_recycled, + "OP ", + ) po.ovarre( self.outfile, "Fuel burn-up rate (reactions/s)", From 5a223063b92e88a5e8b5896452b96dd249ac5858 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 15:08:52 +0000 Subject: [PATCH 17/36] Add plasma fuelling equations and documentation --- .../source/physics-models/plasma_fuelling.md | 44 +++++++++++++++++++ mkdocs.yml | 1 + 2 files changed, 45 insertions(+) create mode 100644 documentation/source/physics-models/plasma_fuelling.md diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md new file mode 100644 index 0000000000..c467232a1d --- /dev/null +++ b/documentation/source/physics-models/plasma_fuelling.md @@ -0,0 +1,44 @@ +# Plasma Fuelling + +The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$. + +$$ +\frac{dn_{\text{T}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{T}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} +$$ + +$$ +\frac{dn_{\text{D}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{T}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} +$$ + +$$ +\frac{dn_{\text{3He}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{3He}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +$$ + +$$ +\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*} +$$ + +In a steady state equilibrium all 4 of these equations should balance, therefore: + +$$ +\frac{dn_{\text{D}}}{dt} = \frac{dn_{\text{T}}}{dt} = \frac{dn_{\text{3He}}}{dt} = \frac{dn_{\alpha}}{dt} = 0 +$$ + +Here $\eta_{\text{fuelling}}$ is the fuelling efficiecny which represents the method of injecting fuel into the plasma. Gas puffing on the low field side is probably around 0.01-0.1, supersonic gas is 0.1 and 0.2 and using pellets can get you close to unity with 0.5-0.9. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel, so $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ together presents the fraction of injected fuel that actually makes it into the plasma core to fuse. + + - $N$ is the total amount of ions in the plasma. + + - $\tau_{\text{fuel}}^*$ is the recycling corrected fuel particle confinement time given by: + + $\tau_{\text{fuel}}^* = (\tau_p) / (1-R)$ + + +Where $\tau_p$ is the particle confinement time which we can assume is approximately equal to the energy confinement time ($\tau_p = \tau_E$). + +The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1]. + + + + +[^1]: G. L. Jackson, V. S. Chan, and R. D. Stambaugh, “An Analytic Expression for the Tritium Burnup Fraction in Burning-Plasma Devices,” Fusion Science and Technology, vol. 64, no. 1, pp. 8–12, Jul. 2013, doi: https://doi.org/10.13182/fst13-a17042. +‌ \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index 6ccdc0a653..bd75f39961 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -52,6 +52,7 @@ nav: - Overview: physics-models/plasma_beta/plasma_beta.md - Fast Alpha: physics-models/plasma_beta/plasma_alpha_beta_contribution.md - Density Limit: physics-models/plasma_density.md + - Fuelling: physics-models/plasma_fuelling.md - Composition & Impurities: physics-models/plasma_composition.md - Radiation: physics-models/plasma_radiation.md - Plasma Current: From e9558aaa1c7e1b53f9e7d25b502f018540f14e40 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 15:35:12 +0000 Subject: [PATCH 18/36] Add D-3He fusion rate and neutron production rate calculations --- process/core/io/plot_proc.py | 22 ++++++++++++--------- process/data_structure/physics_variables.py | 10 ++++++++++ process/models/physics/fusion_reactions.py | 4 ++++ process/models/physics/physics.py | 18 +++++++++++++++++ 4 files changed, 45 insertions(+), 9 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 360ac89d3b..8776cc4b7c 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -11359,7 +11359,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): textstr_dd = ( f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" - f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" + f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n\n" f"D+D -> T fusion rate: {mfile.get('fusrat_plasma_dd_triton', scan=scan):.4e} reactions/s \n" f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s \n" f"Total D-D fusion rate: {mfile.get('fusrat_plasma_dd_total', scan=scan):.4e} reactions/s" @@ -11367,7 +11367,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): axis.text( 0.05, - 0.65, + 0.625, textstr_dd, fontsize=9, verticalalignment="bottom", @@ -11392,11 +11392,14 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): # ================================================= - textstr_dhe3 = f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.4f} MW \n\n" + textstr_dhe3 = ( + f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.4f} MW \n" + f"D+3He fusion rate: {mfile.get('fusrat_plasma_dhe3', scan=scan):.4e} reactions/s \n" + ) axis.text( 0.05, - 0.55, + 0.525, textstr_dhe3, fontsize=9, verticalalignment="bottom", @@ -11411,8 +11414,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.21, - 0.59, + 0.285, + 0.56, "$\\text{D - 3He}$", fontsize=20, verticalalignment="top", @@ -11465,7 +11468,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.4f} MW\n" f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.4f} MW\n\n" f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m$^3$\n" - f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m$^3$" + f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m$^3$\n\n" + f"Neutron production rate: {mfile.get('fusrat_neutron_production_total', scan=scan):.4e} particles/s" ) axis.text( @@ -11485,8 +11489,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.25, - 0.2, + 0.26, + 0.21, "$n$", fontsize=20, verticalalignment="top", diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 7b3b27c1be..7e18c1890d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -473,6 +473,12 @@ fusrat_plasma_dd_total: float = None """Total D-D fusion reaction rate in plasma, (reactions/sec)""" +fusrat_plasma_dhe3: float = None +"""D-3He fusion reaction rate in plasma, (reactions/sec)""" + +fusrat_neutron_production_total: float = None +"""Total neutron production rate from plasma and beams (neutrons/sec)""" + fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -1522,6 +1528,8 @@ def init_physics_variables(): fusrat_plasma_dd_helion, \ fusrat_plasma_dd_triton, \ fusrat_plasma_dd_total, \ + fusrat_plasma_dhe3, \ + fusrat_neutron_production_total, \ fusrat_plasma_dt_profile, \ fusrat_plasma_dd_triton_profile, \ fusrat_plasma_dd_helion_profile, \ @@ -1818,6 +1826,8 @@ def init_physics_variables(): fusrat_plasma_dd_helion = 0.0 fusrat_plasma_dd_triton = 0.0 fusrat_plasma_dd_total = 0.0 + fusrat_plasma_dhe3 = 0.0 + fusrat_neutron_production_total = 0.0 fusrat_dt_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 7a5513c23e..6c9be402af 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -328,6 +328,10 @@ def dhe3_reaction(self): alpha_rate_density = fusion_rate_density proton_rate_density = fusion_rate_density # Proton production rate [m^3/second] + physics_variables.fusrat_plasma_dhe3 = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-3He power density self.dhe3_power_density = fusion_power_density diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 6e05af8932..3f4364cb4d 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -680,6 +680,10 @@ def run(self): + physics_variables.fusrat_plasma_dd_triton ) + physics_variables.fusrat_neutron_production_total = ( + physics_variables.fusrat_plasma_dd_helion + physics_variables.fusrat_dt_total + ) + physics_variables.fusrat_dt_total = ( physics_variables.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) ) @@ -2689,6 +2693,20 @@ def outplas(self): physics_variables.fusrat_plasma_dd_total, "OP ", ) + po.ovarre( + self.outfile, + "D-3He Fusion rate: total (reactions/sec)", + "(fusrat_plasma_dhe3)", + physics_variables.fusrat_plasma_dhe3, + "OP ", + ) + po.ovarre( + self.outfile, + "Neutron production rate: total (particles/sec)", + "(fusrat_neutron_production_total)", + physics_variables.fusrat_neutron_production_total, + "OP ", + ) po.ovarre( self.outfile, "Fusion rate density: total (reactions/m3/sec)", From f1667a022ca2783e2c2a53d9bb580af09fa7a403 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 16:42:25 +0000 Subject: [PATCH 19/36] Add fuelling composition variables and constraints for deuterium, tritium, and helium-3 --- .../source/physics-models/plasma_fuelling.md | 9 +++-- process/core/input.py | 9 +++++ process/core/solver/constraints.py | 36 ++++++++++++++----- process/core/solver/iteration_variables.py | 18 ++++++++++ process/data_structure/numerics.py | 5 +-- process/data_structure/physics_variables.py | 15 ++++++++ process/models/physics/physics.py | 21 +++++++++++ 7 files changed, 100 insertions(+), 13 deletions(-) diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md index c467232a1d..e43a965667 100644 --- a/documentation/source/physics-models/plasma_fuelling.md +++ b/documentation/source/physics-models/plasma_fuelling.md @@ -3,15 +3,15 @@ The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$. $$ -\frac{dn_{\text{T}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{T}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} +\frac{dn_{\text{T}}}{dt} = f_{\text{fuel,T}}\eta_{\text{fuelling}}\Gamma_{\text{T}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} $$ $$ -\frac{dn_{\text{D}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{T}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} +\frac{dn_{\text{D}}}{dt} = f_{\text{fuel,D}}\eta_{\text{fuelling}}\Gamma_{\text{T}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} $$ $$ -\frac{dn_{\text{3He}}}{dt} = \eta_{\text{fuelling}}\Gamma_{\text{3He}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +\frac{dn_{\text{3He}}}{dt} = f_{\text{fuel,3He}}\eta_{\text{fuelling}}\Gamma_{\text{3He}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} $$ $$ @@ -26,6 +26,9 @@ $$ Here $\eta_{\text{fuelling}}$ is the fuelling efficiecny which represents the method of injecting fuel into the plasma. Gas puffing on the low field side is probably around 0.01-0.1, supersonic gas is 0.1 and 0.2 and using pellets can get you close to unity with 0.5-0.9. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel, so $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ together presents the fraction of injected fuel that actually makes it into the plasma core to fuse. + +The fuelling fractional compositions is given by $f$ + - $N$ is the total amount of ions in the plasma. - $\tau_{\text{fuel}}^*$ is the recycling corrected fuel particle confinement time given by: diff --git a/process/core/input.py b/process/core/input.py index 2b55e7b3a3..476c8c4395 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1217,6 +1217,15 @@ def __post_init__(self): "molflow_plasma_fuelling_vv_injected": InputVariable( data_structure.physics_variables, float, range=(1e18, 1e24) ), + "f_molflow_plasma_fuelling_deuterium": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "f_molflow_plasma_fuelling_tritium": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "f_molflow_plasma_fuelling_helium3": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), "pflux_plant_floor_electric": InputVariable( data_structure.heat_transport_variables, float, range=(0.0, 1000.0) ), diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 3f48581128..407ffe5183 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1867,8 +1867,8 @@ def constraint_equation_93(constraint_registration): numerator = ( data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - - data_structure.physics_variables.fusrat_dt_total - ) + * data_structure.physics_variables.f_molflow_plasma_fuelling_tritium + ) + data_structure.physics_variables.fusrat_plasma_dd_triton denominator = ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg * data_structure.physics_variables.vol_plasma @@ -1876,7 +1876,7 @@ def constraint_equation_93(constraint_registration): ) / ( data_structure.physics_variables.t_energy_confinement / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) - ) + ) - data_structure.physics_variables.fusrat_dt_total return eq(numerator, denominator, constraint_registration) @@ -1890,8 +1890,8 @@ def constraint_equation_94(constraint_registration): numerator = ( data_structure.physics_variables.eta_plasma_fuelling * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - - data_structure.physics_variables.fusrat_dt_total - ) + * data_structure.physics_variables.f_molflow_plasma_fuelling_deuterium + ) - data_structure.physics_variables.fusrat_dt_total denominator = ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg * data_structure.physics_variables.vol_plasma @@ -1911,18 +1911,38 @@ def constraint_equation_95(constraint_registration): """ # Alpha particle balance - numerator = data_structure.physics_variables.fusrat_dt_total + numerator = ( + data_structure.physics_variables.fusrat_dt_total + + data_structure.physics_variables.fusrat_plasma_dhe3 + ) denominator = ( data_structure.physics_variables.nd_plasma_alphas_vol_avg * data_structure.physics_variables.vol_plasma ) / ( - (data_structure.physics_variables.t_energy_confinement * 4.0) - / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + data_structure.physics_variables.t_energy_confinement + * data_structure.physics_variables.f_alpha_energy_confinement ) return eq(numerator, denominator, constraint_registration) +@ConstraintManager.register_constraint(96, "", "=") +def constraint_equation_96(constraint_registration): + """Equation for checking the fuelling composition is consistent. + + f_molflow_plasma_fuelling_deuterium: fraction of deuterium ions + f_molflow_plasma_fuelling_tritium: fraction of tritium ions + f_molflow_plasma_fuelling_helium3: fraction of helium-3 ions + """ + return eq( + data_structure.physics_variables.f_molflow_plasma_fuelling_deuterium + + data_structure.physics_variables.f_molflow_plasma_fuelling_tritium + + data_structure.physics_variables.f_molflow_plasma_fuelling_helium3, + 1.0, + constraint_registration, + ) + + def constraint_eqns(m: int, ieqn: int): """Evaluates the constraints given the current state of PROCESS. diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index efbc96a559..fd761b2495 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -303,6 +303,24 @@ class IterationVariable: 1e18, 1e24, ), + 180: IterationVariable( + "f_molflow_plasma_fuelling_deuterium", + data_structure.physics_variables, + 0.0, + 1.0, + ), + 181: IterationVariable( + "f_molflow_plasma_fuelling_tritium", + data_structure.physics_variables, + 0.0, + 1.0, + ), + 182: IterationVariable( + "f_molflow_plasma_fuelling_helium3", + data_structure.physics_variables, + 0.0, + 1.0, + ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 18fa48ee97..2a388e3665 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -1,9 +1,9 @@ import numpy as np -ipnvars: int = 180 +ipnvars: int = 183 """total number of variables available for iteration""" -ipeqns: int = 95 +ipeqns: int = 96 """number of constraint equations available""" ipnfoms: int = 19 @@ -626,6 +626,7 @@ def init_numerics(): "Tritium particle balance consistency ", "Deuterium particle balance consistency ", "Alpha particle balance consistency ", + "Fuelling composition consistency ", ] ixc = np.array([0] * ipnvars) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 7e18c1890d..a21fcf98b1 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1095,6 +1095,15 @@ molflow_plasma_fuelling_vv_injected: float = None """plasma fuelling rate into the vacuum vessel (nucleus-pairs/s)""" +f_molflow_plasma_fuelling_deuterium: float = None +"""fraction of plasma fuelling that is deuterium""" + +f_molflow_plasma_fuelling_tritium: float = None +"""fraction of plasma fuelling that is tritium""" + +f_molflow_plasma_fuelling_helium3: float = None +"""fraction of plasma fuelling that is helium-3""" + tauratio: float = None """tauratio /1.0/ : ratio of He and pellet particle confinement times""" @@ -1656,6 +1665,9 @@ def init_physics_variables(): f_plasma_particles_lcfs_recycled, \ eta_plasma_fuelling, \ molflow_plasma_fuelling_vv_injected, \ + f_molflow_plasma_fuelling_deuterium, \ + f_molflow_plasma_fuelling_tritium, \ + f_molflow_plasma_fuelling_helium3, \ tauratio, \ q95_min, \ qstar, \ @@ -1957,6 +1969,9 @@ def init_physics_variables(): f_plasma_particles_lcfs_recycled = 0.0 eta_plasma_fuelling = 0.0 molflow_plasma_fuelling_vv_injected = 1e20 + f_molflow_plasma_fuelling_deuterium = 0.5 + f_molflow_plasma_fuelling_tritium = 0.5 + f_molflow_plasma_fuelling_helium3 = 0.0 tauratio = 1.0 q95_min = 0.0 qstar = 0.0 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 3f4364cb4d..b9c8b9aaac 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -3329,6 +3329,27 @@ def outplas(self): physics_variables.molflow_plasma_fuelling_vv_injected, "OP ", ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is deuterium", + "(f_molflow_plasma_fuelling_deuterium)", + physics_variables.f_molflow_plasma_fuelling_deuterium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is tritium", + "(f_molflow_plasma_fuelling_tritium)", + physics_variables.f_molflow_plasma_fuelling_tritium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is helium-3", + "(f_molflow_plasma_fuelling_helium3)", + physics_variables.f_molflow_plasma_fuelling_helium3, + "OP ", + ) po.ovarre( self.outfile, "Fuelling efficiency", From 605462de76f01afae1ad894a69582209cfbc70c2 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 17:31:12 +0000 Subject: [PATCH 20/36] Enhance plasma flow calculations by adding tritium and deuterium flow rate parameters, and update plotting functions for alpha particle flow rate and fusion DT rate. --- process/core/solver/constraints.py | 27 ++++++---- process/core/solver/iteration_variables.py | 2 +- process/models/physics/exhaust.py | 57 ++++++++++++++-------- 3 files changed, 57 insertions(+), 29 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 407ffe5183..25629ed7b6 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1865,10 +1865,14 @@ def constraint_equation_93(constraint_registration): # Numerator shall be the tritium particle balance numerator = ( - data_structure.physics_variables.eta_plasma_fuelling - * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - * data_structure.physics_variables.f_molflow_plasma_fuelling_tritium - ) + data_structure.physics_variables.fusrat_plasma_dd_triton + ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + * data_structure.physics_variables.f_molflow_plasma_fuelling_tritium + ) + + data_structure.physics_variables.fusrat_plasma_dd_triton + - data_structure.physics_variables.fusrat_dt_total + ) denominator = ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg * data_structure.physics_variables.vol_plasma @@ -1876,7 +1880,7 @@ def constraint_equation_93(constraint_registration): ) / ( data_structure.physics_variables.t_energy_confinement / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) - ) - data_structure.physics_variables.fusrat_dt_total + ) return eq(numerator, denominator, constraint_registration) @@ -1888,10 +1892,15 @@ def constraint_equation_94(constraint_registration): """ numerator = ( - data_structure.physics_variables.eta_plasma_fuelling - * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected - * data_structure.physics_variables.f_molflow_plasma_fuelling_deuterium - ) - data_structure.physics_variables.fusrat_dt_total + ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + * data_structure.physics_variables.f_molflow_plasma_fuelling_deuterium + ) + - data_structure.physics_variables.fusrat_dt_total + - data_structure.physics_variables.fusrat_plasma_dhe3 + - 2.0 * data_structure.physics_variables.fusrat_plasma_dd_total + ) denominator = ( data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg * data_structure.physics_variables.vol_plasma diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index fd761b2495..870f5837cd 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -301,7 +301,7 @@ class IterationVariable: "molflow_plasma_fuelling_vv_injected", data_structure.physics_variables, 1e18, - 1e24, + 1e22, ), 180: IterationVariable( "f_molflow_plasma_fuelling_deuterium", diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index a1004f55d8..8de99277a6 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -125,6 +125,7 @@ def calculate_eu_demo_re_attachment_metric( @staticmethod def calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium, eta_plasma_fuelling, molflow_plasma_fuelling_vv_injected, fusrat_dt_total, @@ -136,9 +137,12 @@ def calculate_plasma_tritium_flow_rate( f_plasma_fuel_tritium, ): """Calculate the tritium flow rate in the plasma exhaust.""" - # Assuming 50/50 D-T fuel mix, the tritium fuelling rate is half the total fuelling rate, and the tritium loss includes both D and T contributions from fusion reactions. return ( - (eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected / 2) + ( + f_molflow_plasma_fuelling_tritium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) - fusrat_dt_total + fusrat_plasma_dd_triton - ( @@ -149,9 +153,11 @@ def calculate_plasma_tritium_flow_rate( @staticmethod def calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium, eta_plasma_fuelling, molflow_plasma_fuelling_vv_injected, fusrat_dt_total, + fusrat_plasma_dhe3, fusrat_plasma_dd_total, t_energy_confinement, f_plasma_particles_lcfs_recycled, @@ -160,11 +166,15 @@ def calculate_plasma_deuterium_flow_rate( f_plasma_fuel_deuterium, ): """Calculate the deuterium flow rate in the plasma exhaust.""" - # Assuming 50/50 D-T fuel mix, the deuterium fuelling rate is half the total fuelling rate, and the deuterium loss includes both D and T contributions from fusion reactions. return ( - (eta_plasma_fuelling * molflow_plasma_fuelling_vv_injected / 2) + ( + f_molflow_plasma_fuelling_deuterium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) - fusrat_dt_total - 2 * fusrat_plasma_dd_total + - fusrat_plasma_dhe3 - ( (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium) / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) @@ -174,9 +184,9 @@ def calculate_plasma_deuterium_flow_rate( @staticmethod def calculate_plasma_alphas_flow_rate( fusrat_dt_total, + fusrat_plasma_dhe3, t_energy_confinement, f_t_alpha_energy_confinement, - f_plasma_particles_lcfs_recycled, nd_plasma_alphas_vol_avg, vol_plasma, ): @@ -184,9 +194,11 @@ def calculate_plasma_alphas_flow_rate( # Alpha particle balance - return fusrat_dt_total - (nd_plasma_alphas_vol_avg * vol_plasma) / ( - (t_energy_confinement * f_t_alpha_energy_confinement) - / (1 - f_plasma_particles_lcfs_recycled) + return ( + fusrat_dt_total + + fusrat_plasma_dhe3 + - (nd_plasma_alphas_vol_avg * vol_plasma) + / (t_energy_confinement * f_t_alpha_energy_confinement) ) def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): @@ -199,6 +211,9 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): for i, recycling in enumerate(recycling_range): for j, fuelling in enumerate(fuelling_range): tritium_flow[i, j] = self.calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium=mfile.get( + "f_molflow_plasma_fuelling_tritium", scan=scan + ), eta_plasma_fuelling=fuelling, molflow_plasma_fuelling_vv_injected=mfile.get( "molflow_plasma_fuelling_vv_injected", scan=scan @@ -259,11 +274,15 @@ def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int for i, recycling in enumerate(recycling_range): for j, fuelling in enumerate(fuelling_range): deuterium_flow[i, j] = self.calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium=mfile.get( + "f_molflow_plasma_fuelling_deuterium", scan=scan + ), eta_plasma_fuelling=fuelling, molflow_plasma_fuelling_vv_injected=mfile.get( "molflow_plasma_fuelling_vv_injected", scan=scan ), fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), fusrat_plasma_dd_total=mfile.get( "fusrat_plasma_dd_total", scan=scan ), @@ -314,21 +333,21 @@ def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" - recycling_range = np.linspace(0.01, 0.99, 20) - f_t_alpha_energy_confinement_range = np.linspace(4.0, 10.0, 20) + fusion_dt_range = np.linspace(1e19, 5e21, 20) + f_t_alpha_energy_confinement_range = np.linspace(2.0, 10.0, 20) alpha_flow = np.zeros(( - len(recycling_range), + len(fusion_dt_range), len(f_t_alpha_energy_confinement_range), )) - for i, recycling in enumerate(recycling_range): + for i, fusion_dt in enumerate(fusion_dt_range): for j, f_t_alpha_energy_confinement in enumerate( f_t_alpha_energy_confinement_range ): alpha_flow[i, j] = self.calculate_plasma_alphas_flow_rate( - fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_dt_total=fusion_dt, + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), - f_plasma_particles_lcfs_recycled=recycling, nd_plasma_alphas_vol_avg=mfile.get( "nd_plasma_alphas_vol_avg", scan=scan ), @@ -338,14 +357,14 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): contour = axis.contourf( f_t_alpha_energy_confinement_range, - recycling_range, + fusion_dt_range, alpha_flow, levels=15, cmap="RdBu_r", ) axis.contour( f_t_alpha_energy_confinement_range, - recycling_range, + fusion_dt_range, alpha_flow, levels=[0], colors="black", @@ -353,11 +372,11 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): ) # Plot star for mfile values - recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan) f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan) axis.plot( f_t_alpha_mfile, - recycling_mfile, + fusion_dt_mfile, marker="*", markersize=15, color="yellow", @@ -368,7 +387,7 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): axis.set_xlabel( "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" ) - axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]") axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)") axis.minorticks_on() axis.grid(True, which="major", linestyle="-", alpha=0.7) From 39fdfb8988d368f477ccb37124d591c7191320e1 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 21:01:12 +0000 Subject: [PATCH 21/36] Rename burn-up fraction variable to f_plasma_fuel_burnup across multiple files and update related calculations and tests --- process/core/io/plot_proc.py | 2 +- process/data_structure/physics_variables.py | 6 +++--- process/models/physics/physics.py | 16 ++++++++-------- process/models/stellarator/stellarator.py | 2 +- tests/unit/test_physics.py | 4 ++-- 5 files changed, 15 insertions(+), 15 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 8776cc4b7c..b9e0e39f44 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -3005,7 +3005,7 @@ def plot_main_plasma_information( f" - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n" f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n" f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n" - f"Burn-up fraction: {mfile.get('burnup', scan=scan):.4f} \n" + f"Burn-up fraction: {mfile.get('f_plasma_fuel_burnup', scan=scan):.4f} \n" ) axis.text( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index a21fcf98b1..c373098d2d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -253,7 +253,7 @@ """Plasma stored magnetic energy (J)""" -burnup: float = None +f_plasma_fuel_burnup: float = None """fractional plasma burnup""" @@ -1486,7 +1486,7 @@ def init_physics_variables(): b_plasma_toroidal_profile, \ b_plasma_total, \ e_plasma_magnetic_stored, \ - burnup, \ + f_plasma_fuel_burnup, \ burnup_in, \ b_plasma_vertical_required, \ c_beta, \ @@ -1788,7 +1788,7 @@ def init_physics_variables(): b_plasma_toroidal_profile = [] b_plasma_total = 0.0 e_plasma_magnetic_stored = 0.0 - burnup = 0.0 + f_plasma_fuel_burnup = 0.0 burnup_in = 0.0 b_plasma_vertical_required = 0.0 c_beta = 0.5 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index b9c8b9aaac..3abc500e8e 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -942,7 +942,7 @@ def run(self): # Calculate auxiliary physics related information sbar = 1.0e0 ( - physics_variables.burnup, + physics_variables.f_plasma_fuel_burnup, physics_variables.figmer, physics_variables.fusrat, physics_variables.molflow_plasma_fuelling_required, @@ -1473,7 +1473,7 @@ def phyaux( ------- tuple A tuple containing: - - burnup (float): Fractional plasma burnup. + - f_plasma_fuel_burnup (float): Fractional plasma burnup. - figmer (float): Physics figure of merit. - fusrat (float): Number of fusion reactions per second. - molflow_plasma_fuelling_required (float): Fuelling rate for D-T (nucleus-pairs/sec). @@ -1501,12 +1501,12 @@ def phyaux( # The ratio of ash to fuel particle confinement times is given by # tauratio # Possible logic... - # burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3; + # f_plasma_fuel_burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3; # fuel ion-pairs burned/m3 = alpha particles/m3 (for both D-T and D-He3 reactions) # initial fuel ion-pairs/m3 = burnt fuel ion-pairs/m3 + unburnt fuel-ion pairs/m3 # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3 if physics_variables.burnup_in <= 1.0e-9: - burnup = ( + f_plasma_fuel_burnup = ( nd_plasma_alphas_vol_avg / (nd_plasma_alphas_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg) / physics_variables.tauratio @@ -1518,12 +1518,12 @@ def phyaux( rndfuel = fusrat # Required fuelling rate (fuel ion pairs/second) (previously Amps) - molflow_plasma_fuelling_required = rndfuel / burnup + molflow_plasma_fuelling_required = rndfuel / f_plasma_fuel_burnup f_alpha_energy_confinement = t_alpha_confinement / t_energy_confinement return ( - burnup, + f_plasma_fuel_burnup, figmer, fusrat, molflow_plasma_fuelling_required, @@ -3374,8 +3374,8 @@ def outplas(self): po.ovarrf( self.outfile, "Burn-up fraction", - "(burnup)", - physics_variables.burnup, + "(f_plasma_fuel_burnup)", + physics_variables.f_plasma_fuel_burnup, "OP ", ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 1df7adc0b8..b154449ba4 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2357,7 +2357,7 @@ def st_phys(self, output): sbar = 1.0e0 ( - physics_variables.burnup, + physics_variables.f_plasma_fuel_burnup, physics_variables.figmer, _fusrat, physics_variables.molflow_plasma_fuelling_required, diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 9e4b7157a7..3cd17ec6bc 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2251,7 +2251,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): monkeypatch.setattr(physics_variables, "burnup_in", phyauxparam.burnup_in) ( - burnup, + f_plasma_fuel_burnup, figmer, fusrat, molflow_plasma_fuelling_required, @@ -2270,7 +2270,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): vol_plasma=phyauxparam.vol_plasma, ) - assert burnup == pytest.approx(phyauxparam.expected_burnup) + assert f_plasma_fuel_burnup == pytest.approx(phyauxparam.expected_burnup) assert figmer == pytest.approx(phyauxparam.expected_figmer) From 05f62afb11cf8a2176fef9558a13adcfe1738c50 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 21:22:38 +0000 Subject: [PATCH 22/36] Refactor fuel burn-up rate variable to fusrat_total across multiple files and update related calculations and tests --- process/core/io/plot_proc.py | 4 ++-- process/data_structure/physics_variables.py | 6 ------ process/models/costs/costs.py | 2 +- process/models/physics/physics.py | 6 +++--- process/models/stellarator/stellarator.py | 2 +- tests/unit/test_costs_1990.py | 10 +++++----- 6 files changed, 12 insertions(+), 18 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index b9e0e39f44..454cc0fee4 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -3003,8 +3003,8 @@ def plot_main_plasma_information( f" - Average mass of all plasma ions: {mfile.get('m_ions_total_amu', scan=scan):.3f} amu\n" f"Fuel mass: {mfile.get('m_plasma_fuel_ions', scan=scan) * 1000:.4f} g\n" f" - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n" - f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n" - f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n" + f"Fueling rate: {mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.3e} nucleus-pairs/s\n" + f"Fuel burn-up rate: {mfile.get('fusrat_total', scan=scan):.3e} reactions/s \n" f"Burn-up fraction: {mfile.get('f_plasma_fuel_burnup', scan=scan):.4f} \n" ) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index c373098d2d..096e9ea211 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1167,10 +1167,6 @@ """n_carbon / n_e""" -rndfuel: float = None -"""fuel burnup rate (reactions/second)""" - - f_nd_plasma_iron_argon_electron: float = None """n_highZ / n_e""" @@ -1685,7 +1681,6 @@ def init_physics_variables(): rminor, \ f_nd_beam_electron, \ f_nd_plasma_carbon_electron, \ - rndfuel, \ f_nd_plasma_iron_argon_electron, \ f_nd_plasma_oxygen_electron, \ f_res_plasma_neo, \ @@ -1989,7 +1984,6 @@ def init_physics_variables(): rminor = 0.0 f_nd_beam_electron = 0.005 f_nd_plasma_carbon_electron = 0.0 - rndfuel = 0.0 f_nd_plasma_iron_argon_electron = 0.0 f_nd_plasma_oxygen_electron = 0.0 f_res_plasma_neo = 0.0 diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py index be7f2198e3..05e520411b 100644 --- a/process/models/costs/costs.py +++ b/process/models/costs/costs.py @@ -2362,7 +2362,7 @@ def acc2272(self): # New calculation: 2 nuclei * reactions/sec * kg/nucleus * g/kg * sec/day physics_variables.wtgpd = ( 2.0e0 - * physics_variables.rndfuel + * physics_variables.fusrat_total * physics_variables.m_fuel_amu * constants.UMASS * 1000.0e0 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 3abc500e8e..69ad306afd 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -946,7 +946,7 @@ def run(self): physics_variables.figmer, physics_variables.fusrat, physics_variables.molflow_plasma_fuelling_required, - physics_variables.rndfuel, + _, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.phyaux( @@ -3367,8 +3367,8 @@ def outplas(self): po.ovarre( self.outfile, "Fuel burn-up rate (reactions/s)", - "(rndfuel)", - physics_variables.rndfuel, + "(fusrat_total)", + physics_variables.fusrat_total, "OP ", ) po.ovarrf( diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index b154449ba4..d75ca10b16 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2361,7 +2361,7 @@ def st_phys(self, output): physics_variables.figmer, _fusrat, physics_variables.molflow_plasma_fuelling_required, - physics_variables.rndfuel, + _, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.physics.phyaux( diff --git a/tests/unit/test_costs_1990.py b/tests/unit/test_costs_1990.py index 1843269b92..ebeb7c3f8b 100644 --- a/tests/unit/test_costs_1990.py +++ b/tests/unit/test_costs_1990.py @@ -164,7 +164,7 @@ def test_acc2272(monkeypatch, costs): :param monkeypatch: Mock fixture :type monkeypatch: object """ - monkeypatch.setattr(physics_variables, "rndfuel", 7.158e20) + monkeypatch.setattr(physics_variables, "fusrat_total", 7.158e20) monkeypatch.setattr(physics_variables, "m_fuel_amu", 2.5) monkeypatch.setattr(cost_variables, "fkind", 1) monkeypatch.setattr(cost_variables, "c2271", 0) @@ -4380,7 +4380,7 @@ class Acc2272Param(NamedTuple): wtgpd: Any = None - rndfuel: Any = None + fusrat_total: Any = None m_fuel_amu: Any = None @@ -4406,7 +4406,7 @@ class Acc2272Param(NamedTuple): gain=0, edrive=5000000, wtgpd=0, - rndfuel=7.0799717510383796e20, + fusrat_total=7.0799717510383796e20, m_fuel_amu=2.5, c227=0, c2272=0, @@ -4422,7 +4422,7 @@ class Acc2272Param(NamedTuple): gain=0, edrive=5000000, wtgpd=507.88376577416528, - rndfuel=7.0777619721108953e20, + fusrat_total=7.0777619721108953e20, m_fuel_amu=2.5, c227=284.96904049038437, c2272=114.02873340990777, @@ -4459,7 +4459,7 @@ def test_acc2272_rut(acc2272param, monkeypatch, costs): monkeypatch.setattr(physics_variables, "wtgpd", acc2272param.wtgpd) - monkeypatch.setattr(physics_variables, "rndfuel", acc2272param.rndfuel) + monkeypatch.setattr(physics_variables, "fusrat_total", acc2272param.fusrat_total) monkeypatch.setattr(physics_variables, "m_fuel_amu", acc2272param.m_fuel_amu) From 0baf1c1a7e07f6eb0167107808f6d6b943b068ac Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 11:21:44 +0000 Subject: [PATCH 23/36] Refactor plasma exhaust calculations to plasma fuelling, updating related plotting functions for tritium, deuterium, and alpha particle flow rates --- process/core/io/plot_proc.py | 2 +- process/models/physics/exhaust.py | 275 --------------------------- process/models/physics/fuelling.py | 288 +++++++++++++++++++++++++++++ 3 files changed, 289 insertions(+), 276 deletions(-) create mode 100644 process/models/physics/fuelling.py diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 454cc0fee4..83350b9712 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -67,7 +67,7 @@ ElectronBernstein, ElectronCyclotron, ) -from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import read_impurity_file from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel from process.models.physics.plasma_current import PlasmaCurrent, PlasmaCurrentModel diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index 8de99277a6..eb0021fb8e 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -1,9 +1,5 @@ import logging -import matplotlib.pyplot as plt -import numpy as np - -import process.core.io.mfile as mf from process.core import constants logger = logging.getLogger(__name__) @@ -122,274 +118,3 @@ def calculate_eu_demo_re_attachment_metric( return (p_plasma_separatrix_mw * b_plasma_toroidal_on_axis) / ( q95 * aspect * rmajor ) - - @staticmethod - def calculate_plasma_tritium_flow_rate( - f_molflow_plasma_fuelling_tritium, - eta_plasma_fuelling, - molflow_plasma_fuelling_vv_injected, - fusrat_dt_total, - fusrat_plasma_dd_triton, - t_energy_confinement, - f_plasma_particles_lcfs_recycled, - nd_plasma_fuel_ions_vol_avg, - vol_plasma, - f_plasma_fuel_tritium, - ): - """Calculate the tritium flow rate in the plasma exhaust.""" - return ( - ( - f_molflow_plasma_fuelling_tritium - * eta_plasma_fuelling - * molflow_plasma_fuelling_vv_injected - ) - - fusrat_dt_total - + fusrat_plasma_dd_triton - - ( - (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) - / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) - ) - ) - - @staticmethod - def calculate_plasma_deuterium_flow_rate( - f_molflow_plasma_fuelling_deuterium, - eta_plasma_fuelling, - molflow_plasma_fuelling_vv_injected, - fusrat_dt_total, - fusrat_plasma_dhe3, - fusrat_plasma_dd_total, - t_energy_confinement, - f_plasma_particles_lcfs_recycled, - nd_plasma_fuel_ions_vol_avg, - vol_plasma, - f_plasma_fuel_deuterium, - ): - """Calculate the deuterium flow rate in the plasma exhaust.""" - return ( - ( - f_molflow_plasma_fuelling_deuterium - * eta_plasma_fuelling - * molflow_plasma_fuelling_vv_injected - ) - - fusrat_dt_total - - 2 * fusrat_plasma_dd_total - - fusrat_plasma_dhe3 - - ( - (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium) - / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) - ) - ) - - @staticmethod - def calculate_plasma_alphas_flow_rate( - fusrat_dt_total, - fusrat_plasma_dhe3, - t_energy_confinement, - f_t_alpha_energy_confinement, - nd_plasma_alphas_vol_avg, - vol_plasma, - ): - """Calculate the alpha particle flow rate in the plasma exhaust.""" - - # Alpha particle balance - - return ( - fusrat_dt_total - + fusrat_plasma_dhe3 - - (nd_plasma_alphas_vol_avg * vol_plasma) - / (t_energy_confinement * f_t_alpha_energy_confinement) - ) - - def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): - """Plot contour of tritium flow rate vs recycling and fuelling rate.""" - - recycling_range = np.linspace(0.01, 0.99, 20) - fuelling_range = np.linspace(0.01, 1.0, 20) - tritium_flow = np.zeros((len(recycling_range), len(fuelling_range))) - - for i, recycling in enumerate(recycling_range): - for j, fuelling in enumerate(fuelling_range): - tritium_flow[i, j] = self.calculate_plasma_tritium_flow_rate( - f_molflow_plasma_fuelling_tritium=mfile.get( - "f_molflow_plasma_fuelling_tritium", scan=scan - ), - eta_plasma_fuelling=fuelling, - molflow_plasma_fuelling_vv_injected=mfile.get( - "molflow_plasma_fuelling_vv_injected", scan=scan - ), - fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), - fusrat_plasma_dd_triton=mfile.get( - "fusrat_plasma_dd_triton", scan=scan - ), - t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), - f_plasma_particles_lcfs_recycled=recycling, - nd_plasma_fuel_ions_vol_avg=mfile.get( - "nd_plasma_fuel_ions_vol_avg", scan=scan - ), - vol_plasma=mfile.get("vol_plasma", scan=scan), - f_plasma_fuel_tritium=mfile.get("f_plasma_fuel_tritium", scan=scan), - ) - - contour = axis.contourf( - fuelling_range, recycling_range, tritium_flow, levels=15, cmap="RdBu_r" - ) - axis.contour( - fuelling_range, - recycling_range, - tritium_flow, - levels=[0], - colors="black", - linewidths=2, - ) - - # Plot star for mfile values - recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) - fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) - axis.plot( - fuelling_mfile, - recycling_mfile, - marker="*", - markersize=15, - color="yellow", - markeredgecolor="black", - markeredgewidth=1.5, - ) - - axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") - axis.set_ylabel("Recycling Fraction [$R$]") - axis.set_title("Plasma Tritium Flow Rate (particles/s)") - axis.minorticks_on() - axis.grid(True, which="major", linestyle="-", alpha=0.7) - axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") - - def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): - """Plot contour of deuterium flow rate vs recycling and fuelling rate.""" - - recycling_range = np.linspace(0.01, 0.99, 20) - fuelling_range = np.linspace(0.01, 1.0, 20) - deuterium_flow = np.zeros((len(recycling_range), len(fuelling_range))) - - for i, recycling in enumerate(recycling_range): - for j, fuelling in enumerate(fuelling_range): - deuterium_flow[i, j] = self.calculate_plasma_deuterium_flow_rate( - f_molflow_plasma_fuelling_deuterium=mfile.get( - "f_molflow_plasma_fuelling_deuterium", scan=scan - ), - eta_plasma_fuelling=fuelling, - molflow_plasma_fuelling_vv_injected=mfile.get( - "molflow_plasma_fuelling_vv_injected", scan=scan - ), - fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), - fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), - fusrat_plasma_dd_total=mfile.get( - "fusrat_plasma_dd_total", scan=scan - ), - t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), - f_plasma_particles_lcfs_recycled=recycling, - nd_plasma_fuel_ions_vol_avg=mfile.get( - "nd_plasma_fuel_ions_vol_avg", scan=scan - ), - vol_plasma=mfile.get("vol_plasma", scan=scan), - f_plasma_fuel_deuterium=mfile.get( - "f_plasma_fuel_deuterium", scan=scan - ), - ) - - contour = axis.contourf( - fuelling_range, recycling_range, deuterium_flow, levels=15, cmap="RdBu_r" - ) - axis.contour( - fuelling_range, - recycling_range, - deuterium_flow, - levels=[0], - colors="black", - linewidths=2, - ) - - # Plot star for mfile values - recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) - fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) - axis.plot( - fuelling_mfile, - recycling_mfile, - marker="*", - markersize=15, - color="yellow", - markeredgecolor="black", - markeredgewidth=1.5, - ) - - axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") - axis.set_ylabel("Recycling Fraction [$R$]") - axis.set_title("Plasma Deuterium Flow Rate (particles/s)") - axis.minorticks_on() - axis.grid(True, which="major", linestyle="-", alpha=0.7) - axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") - - def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): - """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" - - fusion_dt_range = np.linspace(1e19, 5e21, 20) - f_t_alpha_energy_confinement_range = np.linspace(2.0, 10.0, 20) - alpha_flow = np.zeros(( - len(fusion_dt_range), - len(f_t_alpha_energy_confinement_range), - )) - - for i, fusion_dt in enumerate(fusion_dt_range): - for j, f_t_alpha_energy_confinement in enumerate( - f_t_alpha_energy_confinement_range - ): - alpha_flow[i, j] = self.calculate_plasma_alphas_flow_rate( - fusrat_dt_total=fusion_dt, - fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), - t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), - nd_plasma_alphas_vol_avg=mfile.get( - "nd_plasma_alphas_vol_avg", scan=scan - ), - vol_plasma=mfile.get("vol_plasma", scan=scan), - f_t_alpha_energy_confinement=f_t_alpha_energy_confinement, - ) - - contour = axis.contourf( - f_t_alpha_energy_confinement_range, - fusion_dt_range, - alpha_flow, - levels=15, - cmap="RdBu_r", - ) - axis.contour( - f_t_alpha_energy_confinement_range, - fusion_dt_range, - alpha_flow, - levels=[0], - colors="black", - linewidths=2, - ) - - # Plot star for mfile values - fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan) - f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan) - axis.plot( - f_t_alpha_mfile, - fusion_dt_mfile, - marker="*", - markersize=15, - color="yellow", - markeredgecolor="black", - markeredgewidth=1.5, - ) - - axis.set_xlabel( - "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" - ) - axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]") - axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)") - axis.minorticks_on() - axis.grid(True, which="major", linestyle="-", alpha=0.7) - axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py new file mode 100644 index 0000000000..704fe503a2 --- /dev/null +++ b/process/models/physics/fuelling.py @@ -0,0 +1,288 @@ +import logging + +import matplotlib.pyplot as plt +import numpy as np + +import process.core.io.mfile as mf +from process.core import constants + +logger = logging.getLogger(__name__) + + +class PlasmaFuelling: + """Class to hold plasma fuelling calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + @staticmethod + def calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium, + eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected, + fusrat_dt_total, + fusrat_plasma_dd_triton, + t_energy_confinement, + f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg, + vol_plasma, + f_plasma_fuel_tritium, + ): + """Calculate the tritium flow rate in the plasma exhaust.""" + return ( + ( + f_molflow_plasma_fuelling_tritium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + - fusrat_dt_total + + fusrat_plasma_dd_triton + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium, + eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected, + fusrat_dt_total, + fusrat_plasma_dhe3, + fusrat_plasma_dd_total, + t_energy_confinement, + f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg, + vol_plasma, + f_plasma_fuel_deuterium, + ): + """Calculate the deuterium flow rate in the plasma exhaust.""" + return ( + ( + f_molflow_plasma_fuelling_deuterium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + - fusrat_dt_total + - 2 * fusrat_plasma_dd_total + - fusrat_plasma_dhe3 + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_alphas_flow_rate( + fusrat_dt_total, + fusrat_plasma_dhe3, + t_energy_confinement, + f_t_alpha_energy_confinement, + nd_plasma_alphas_vol_avg, + vol_plasma, + ): + """Calculate the alpha particle flow rate in the plasma exhaust.""" + + # Alpha particle balance + + return ( + fusrat_dt_total + + fusrat_plasma_dhe3 + - (nd_plasma_alphas_vol_avg * vol_plasma) + / (t_energy_confinement * f_t_alpha_energy_confinement) + ) + + def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of tritium flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + tritium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + tritium_flow[i, j] = self.calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium=mfile.get( + "f_molflow_plasma_fuelling_tritium", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dd_triton=mfile.get( + "fusrat_plasma_dd_triton", scan=scan + ), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_tritium=mfile.get("f_plasma_fuel_tritium", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, recycling_range, tritium_flow, levels=15, cmap="RdBu_r" + ) + axis.contour( + fuelling_range, + recycling_range, + tritium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Tritium Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") + + def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of deuterium flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + deuterium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + deuterium_flow[i, j] = self.calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium=mfile.get( + "f_molflow_plasma_fuelling_deuterium", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + fusrat_plasma_dd_total=mfile.get( + "fusrat_plasma_dd_total", scan=scan + ), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_deuterium=mfile.get( + "f_plasma_fuel_deuterium", scan=scan + ), + ) + + contour = axis.contourf( + fuelling_range, recycling_range, deuterium_flow, levels=15, cmap="RdBu_r" + ) + axis.contour( + fuelling_range, + recycling_range, + deuterium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Deuterium Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") + + def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" + + fusion_dt_range = np.linspace(1e19, 5e21, 20) + f_t_alpha_energy_confinement_range = np.linspace(2.0, 10.0, 20) + alpha_flow = np.zeros(( + len(fusion_dt_range), + len(f_t_alpha_energy_confinement_range), + )) + + for i, fusion_dt in enumerate(fusion_dt_range): + for j, f_t_alpha_energy_confinement in enumerate( + f_t_alpha_energy_confinement_range + ): + alpha_flow[i, j] = self.calculate_plasma_alphas_flow_rate( + fusrat_dt_total=fusion_dt, + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + nd_plasma_alphas_vol_avg=mfile.get( + "nd_plasma_alphas_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_t_alpha_energy_confinement=f_t_alpha_energy_confinement, + ) + + contour = axis.contourf( + f_t_alpha_energy_confinement_range, + fusion_dt_range, + alpha_flow, + levels=15, + cmap="RdBu_r", + ) + axis.contour( + f_t_alpha_energy_confinement_range, + fusion_dt_range, + alpha_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan) + f_t_alpha_mfile = mfile.get("f_alpha_energy_confinement", scan=scan) + axis.plot( + f_t_alpha_mfile, + fusion_dt_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel( + "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" + ) + axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]") + axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)") + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") From 17b4b687b8be859edb8871200ac9c0b0e09ffda0 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 11:59:43 +0000 Subject: [PATCH 24/36] Refactor plasma fuelling plotting functions to improve organization and add fuelling information display --- process/models/physics/fuelling.py | 218 +++++++++++++++++++++++++---- 1 file changed, 193 insertions(+), 25 deletions(-) diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 704fe503a2..7a5c82d64f 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -16,20 +16,72 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + @staticmethod + def calculate_fuel_burnup_fraction( + fusrat_total: float, molflow_plasma_fuelling_vv_injected: float + ) -> float: + """Calculate the fuel burnup fraction + + Parameters + ---------- + fusrat_total : float + Total fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + + Returns + ------- + float Fuel burnup fraction (dimensionless). + + """ + + return fusrat_total / molflow_plasma_fuelling_vv_injected + @staticmethod def calculate_plasma_tritium_flow_rate( - f_molflow_plasma_fuelling_tritium, - eta_plasma_fuelling, - molflow_plasma_fuelling_vv_injected, - fusrat_dt_total, - fusrat_plasma_dd_triton, - t_energy_confinement, - f_plasma_particles_lcfs_recycled, - nd_plasma_fuel_ions_vol_avg, - vol_plasma, - f_plasma_fuel_tritium, - ): - """Calculate the tritium flow rate in the plasma exhaust.""" + f_molflow_plasma_fuelling_tritium: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_dt_total: float, + fusrat_plasma_dd_triton: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_tritium: float, + ) -> float: + """Calculate the tritium flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_tritium : float + Fraction of tritium in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dd_triton : float + Tritium production rate from DD fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_tritium : float + Fraction of tritium in the plasma fuel. + + Returns + ------- + float + Tritium flow rate in the plasma exhaust (particles/s). + + """ + return ( ( f_molflow_plasma_fuelling_tritium @@ -46,19 +98,52 @@ def calculate_plasma_tritium_flow_rate( @staticmethod def calculate_plasma_deuterium_flow_rate( - f_molflow_plasma_fuelling_deuterium, - eta_plasma_fuelling, - molflow_plasma_fuelling_vv_injected, - fusrat_dt_total, - fusrat_plasma_dhe3, - fusrat_plasma_dd_total, - t_energy_confinement, - f_plasma_particles_lcfs_recycled, - nd_plasma_fuel_ions_vol_avg, - vol_plasma, - f_plasma_fuel_deuterium, - ): - """Calculate the deuterium flow rate in the plasma exhaust.""" + f_molflow_plasma_fuelling_deuterium: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_dt_total: float, + fusrat_plasma_dhe3: float, + fusrat_plasma_dd_total: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_deuterium: float, + ) -> float: + """Calculate the deuterium flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_deuterium : float + Fraction of deuterium in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + fusrat_plasma_dd_total : float + Total deuterium consumption rate from DD fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_deuterium : float + Fraction of deuterium in the plasma fuel. + + Returns + ------- + float + Deuterium flow rate in the plasma exhaust (particles/s). + + + """ return ( ( f_molflow_plasma_fuelling_deuterium @@ -74,6 +159,61 @@ def calculate_plasma_deuterium_flow_rate( ) ) + @staticmethod + def calculate_plasma_helium3_flow_rate( + f_molflow_plasma_fuelling_helium3: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_plasma_dhe3: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_helium3: float, + ) -> float: + """Calculate the helium-3 flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_helium3 : float + Fraction of helium-3 in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_helium3 : float + Fraction of helium-3 in the plasma fuel. + + Returns + ------- + float + Helium-3 flow rate in the plasma exhaust (particles/s). + + """ + + return ( + ( + f_molflow_plasma_fuelling_helium3 + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + + fusrat_plasma_dhe3 + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_helium3) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + @staticmethod def calculate_plasma_alphas_flow_rate( fusrat_dt_total, @@ -286,3 +426,31 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): axis.grid(True, which="major", linestyle="-", alpha=0.7) axis.grid(True, which="minor", linestyle=":", alpha=0.4) plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") + + def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): + """Plot fuelling information.""" + msg = ( + f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n" + f"Total fuelling rate:" + f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n" + f"Fuelling Rate Efficiency ($\\eta_{{\\text{{fuelling}}}}$): " + f"{mfile.get('eta_plasma_fuelling', scan=scan):.4f}\n" + f"Recycling Fraction ($R$): " + f"{mfile.get('f_plasma_particles_lcfs_recycled', scan=scan):.4f}\n\n" + f"Fraction of Tritium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_tritium', scan=scan):.4f}\n" + f"Fraction of Deuterium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_deuterium', scan=scan):.4f}\n" + f"Fraction of 3-Helium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}" + ) + fig.text( + 0.75, + 0.25, + msg, + ha="center", + va="center", + transform=fig.transFigure, + fontsize=9, + bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0}, + ) From dd7fa2bc2e5a172101149cd427ada9daf6ea0ce6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 12:24:48 +0000 Subject: [PATCH 25/36] Integrate PlasmaFuelling class into physics calculations and update fuel burnup fraction computation across models --- process/main.py | 3 + process/models/physics/fuelling.py | 4 +- process/models/physics/physics.py | 68 +++-------------------- process/models/stellarator/stellarator.py | 17 ++---- 4 files changed, 20 insertions(+), 72 deletions(-) diff --git a/process/main.py b/process/main.py index c37e24dadb..32b1fa6f02 100644 --- a/process/main.py +++ b/process/main.py @@ -97,6 +97,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( @@ -709,6 +710,7 @@ def __init__(self): self.plasma_inductance = PlasmaInductance() self.plasma_density_limit = PlasmaDensityLimit() self.plasma_exhaust = PlasmaExhaust() + self.plasma_fuelling = PlasmaFuelling() self.plasma_bootstrap_current = PlasmaBootstrapCurrent( plasma_profile=self.plasma_profile ) @@ -726,6 +728,7 @@ def __init__(self): plasma_confinement=self.plasma_confinement, plasma_transition=self.plasma_transition, plasma_current=self.plasma_current, + plasma_fuelling=self.plasma_fuelling, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 7a5c82d64f..9df1ea57cb 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -442,7 +442,9 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): f"Fraction of Deuterium Fuelling: " f"{mfile.get('f_molflow_plasma_fuelling_deuterium', scan=scan):.4f}\n" f"Fraction of 3-Helium Fuelling: " - f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}" + f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}\n\n" + f"Fuel Burnup Fraction: " + f"{mfile.get('f_plasma_fuel_burnup', scan=scan):.4f}\n" ) fig.text( 0.75, diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 69ad306afd..c742b3994f 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -30,6 +30,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.plasma_current import PlasmaCurrent @@ -214,6 +215,7 @@ def __init__( plasma_confinement: PlasmaConfinementTime, plasma_transition: PlasmaConfinementTransition, plasma_current: PlasmaCurrent, + plasma_fuelling: PlasmaFuelling, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -227,6 +229,7 @@ def __init__( self.confinement = plasma_confinement self.plasma_transition = plasma_transition self.current = plasma_current + self.fuelling = plasma_fuelling def output(self): self.calculate_effective_charge_ionisation_profiles() @@ -939,26 +942,18 @@ def run(self): 0.5e0 * physics_variables.ind_plasma * physics_variables.plasma_current**2 ) - # Calculate auxiliary physics related information - sbar = 1.0e0 ( - physics_variables.f_plasma_fuel_burnup, - physics_variables.figmer, - physics_variables.fusrat, - physics_variables.molflow_plasma_fuelling_required, - _, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.phyaux( - physics_variables.aspect, - physics_variables.nd_plasma_fuel_ions_vol_avg, - physics_variables.fusden_total, physics_variables.fusden_alpha_total, - physics_variables.plasma_current, - sbar, physics_variables.nd_plasma_alphas_vol_avg, physics_variables.t_energy_confinement, - physics_variables.vol_plasma, + ) + + physics_variables.f_plasma_fuel_burnup = self.fuelling.calculate_fuel_burnup_fraction( + fusrat_total=physics_variables.fusrat_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, ) physics_variables.ntau, physics_variables.nTtau = ( @@ -1436,16 +1431,10 @@ def plasma_composition(): @staticmethod @nb.njit(cache=True) def phyaux( - aspect: float, - nd_plasma_fuel_ions_vol_avg: float, - fusden_total: float, fusden_alpha_total: float, - plasma_current: float, - sbar: float, nd_plasma_alphas_vol_avg: float, t_energy_confinement: float, - vol_plasma: float, - ) -> tuple[float, float, float, float, float, float, float, float]: + ) -> tuple[float, float]: """Auxiliary physics quantities Parameters @@ -1473,20 +1462,11 @@ def phyaux( ------- tuple A tuple containing: - - f_plasma_fuel_burnup (float): Fractional plasma burnup. - - figmer (float): Physics figure of merit. - - fusrat (float): Number of fusion reactions per second. - - molflow_plasma_fuelling_required (float): Fuelling rate for D-T (nucleus-pairs/sec). - - rndfuel (float): Fuel burnup rate (reactions/s). - t_alpha_confinement (float): Alpha particle confinement time (s). - f_alpha_energy_confinement (float): Fraction of alpha energy confinement. This subroutine calculates extra physics related items needed by other parts of the code. """ - figmer = 1e-6 * plasma_current * aspect**sbar - - # Fusion reactions per second - fusrat = fusden_total * vol_plasma # Alpha particle confinement time (s) # Number of alphas / alpha production rate @@ -1495,39 +1475,9 @@ def phyaux( else: # only likely if DD is only active fusion reaction t_alpha_confinement = 0.0 - # Fractional burnup - # (Consider detailed model in: G. L. Jackson, V. S. Chan, R. D. Stambaugh, - # Fusion Science and Technology, vol.64, no.1, July 2013, pp.8-12) - # The ratio of ash to fuel particle confinement times is given by - # tauratio - # Possible logic... - # f_plasma_fuel_burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3; - # fuel ion-pairs burned/m3 = alpha particles/m3 (for both D-T and D-He3 reactions) - # initial fuel ion-pairs/m3 = burnt fuel ion-pairs/m3 + unburnt fuel-ion pairs/m3 - # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3 - if physics_variables.burnup_in <= 1.0e-9: - f_plasma_fuel_burnup = ( - nd_plasma_alphas_vol_avg - / (nd_plasma_alphas_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg) - / physics_variables.tauratio - ) - else: - burnup = physics_variables.burnup_in - - # Fuel burnup rate (reactions/second) (previously Amps) - rndfuel = fusrat - - # Required fuelling rate (fuel ion pairs/second) (previously Amps) - molflow_plasma_fuelling_required = rndfuel / f_plasma_fuel_burnup - f_alpha_energy_confinement = t_alpha_confinement / t_energy_confinement return ( - f_plasma_fuel_burnup, - figmer, - fusrat, - molflow_plasma_fuelling_required, - rndfuel, t_alpha_confinement, f_alpha_energy_confinement, ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index d75ca10b16..ff59e1d1c2 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2355,25 +2355,18 @@ def st_phys(self, output): # Calculate auxiliary physics related information # for the rest of the code - sbar = 1.0e0 ( - physics_variables.f_plasma_fuel_burnup, - physics_variables.figmer, - _fusrat, - physics_variables.molflow_plasma_fuelling_required, - _, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.physics.phyaux( - physics_variables.aspect, - physics_variables.nd_plasma_fuel_ions_vol_avg, - physics_variables.fusden_total, physics_variables.fusden_alpha_total, - physics_variables.plasma_current, - sbar, physics_variables.nd_plasma_alphas_vol_avg, physics_variables.t_energy_confinement, - physics_variables.vol_plasma, + ) + + physics_variables.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction( + fusrat_total=physics_variables.fusrat_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, ) # Calculate the neoclassical sanity check with PROCESS parameters From f32e570f22c89ec49b53e2c310f870f1d11a5d21 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 13:40:30 +0000 Subject: [PATCH 26/36] Enhance plasma fuelling documentation with flow rate equations and key constraints for consistency models --- .../source/physics-models/plasma_fuelling.md | 90 ++++++++++++++++++- 1 file changed, 86 insertions(+), 4 deletions(-) diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md index e43a965667..ededb6966c 100644 --- a/documentation/source/physics-models/plasma_fuelling.md +++ b/documentation/source/physics-models/plasma_fuelling.md @@ -1,17 +1,17 @@ -# Plasma Fuelling +# Plasma Fuelling | `PlasmaFuelling()` The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$. $$ -\frac{dn_{\text{T}}}{dt} = f_{\text{fuel,T}}\eta_{\text{fuelling}}\Gamma_{\text{T}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} +\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} $$ $$ -\frac{dn_{\text{D}}}{dt} = f_{\text{fuel,D}}\eta_{\text{fuelling}}\Gamma_{\text{T}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} +\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} $$ $$ -\frac{dn_{\text{3He}}}{dt} = f_{\text{fuel,3He}}\eta_{\text{fuelling}}\Gamma_{\text{3He}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} $$ $$ @@ -42,6 +42,88 @@ The recycling coefficient $R$, defined as the fraction of particles crossing the +-------------- + +## Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()` + +$$ +\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*} +$$ + +--------------- + +## Deuterium Flow Rate | `calculate_plasma_deuterium_flow_rate()` + +$$ +\frac{dn_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*} +$$ + +--------------- + +## Helium-3 Flow Rate | `calculate_plasma_helium3_flow_rate()` + +$$ +\frac{dn_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +$$ + +--------------- + +## Alpha Particle Flow Rate | `calculate_plasma_alphas_flow_rate()` + +$$ +\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*} +$$ + +----------------- + +## Key Constraints + +### Deuterium Flow Consistency + +This constraint can be activated by stating `icc = 94` in the input file. + +This constraint ensures that the change in deuterium particles as a function of time is zero. It ensures the output of `calculate_plasma_deuterium_flow_rate()` is zero + +**It is recommended to have this constraint on as it is a plasma consistency model** + +---------------- + +### Tritium Flow Consistency + +This constraint can be activated by stating `icc = 93` in the input file. + +This constraint ensures that the change in tritium particles as a function of time is zero. It ensures the output of `calculate_plasma_tritium_flow_rate()` is zero + +**It is recommended to have this constraint on as it is a plasma consistency model** + +----------------- + +### Helium-3 Flow Consistency + +### Alpha Particle Flow Consistency + +This constraint can be activated by stating `icc = 95` in the input file. + +This constraint ensures that the change in alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_flow_rate(()` is zero + +**It is recommended to have this constraint on as it is a plasma consistency model** + +------------------ + +### Fuelling Proportion Consistency + +This constraint can be activated by stating `icc = 96` in the input file. + +This ensures that all 3 injected fuelling fractions sum up to 1: + +$$ +f_{\text{fuelling,D}} + f_{\text{fuelling,T}} + f_{\text{fuelling,3He}} = 1.0 +$$ + +**It is recommended to have this constraint on as it is a plasma consistency model** + +----------------- + [^1]: G. L. Jackson, V. S. Chan, and R. D. Stambaugh, “An Analytic Expression for the Tritium Burnup Fraction in Burning-Plasma Devices,” Fusion Science and Technology, vol. 64, no. 1, pp. 8–12, Jul. 2013, doi: https://doi.org/10.13182/fst13-a17042. ‌ \ No newline at end of file From 21c17dd9462459cc6299e38cf3202b2ec11131f8 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 14:56:41 +0000 Subject: [PATCH 27/36] Refactor fuelling output logic to use dedicated method in PlasmaFuelling class and update plotting function calls for consistency --- process/models/physics/fuelling.py | 101 +++++++++++++++++++++++++++++ process/models/physics/physics.py | 93 +------------------------- 2 files changed, 102 insertions(+), 92 deletions(-) diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 9df1ea57cb..ea8879e21a 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -5,6 +5,12 @@ import process.core.io.mfile as mf from process.core import constants +from process.core import process_output as po +from process.data_structure import ( + numerics, + physics_variables, + reinke_variables, +) logger = logging.getLogger(__name__) @@ -456,3 +462,98 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): fontsize=9, bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0}, ) + + def output_fuelling_info(self): + """Output fuelling information to mfile.""" + po.oheadr(self.outfile, "Plasma Fuelling") + po.ovarre( + self.outfile, + "Ratio of He and pellet particle confinement times", + "(tauratio)", + physics_variables.tauratio, + ) + po.ovarre( + self.outfile, + "Fuelling rate (nucleus-pairs/s)", + "(molflow_plasma_fuelling_required)", + physics_variables.molflow_plasma_fuelling_required, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling rate (nucleus-pairs/s)", + "(molflow_plasma_fuelling_vv_injected)", + physics_variables.molflow_plasma_fuelling_vv_injected, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is deuterium", + "(f_molflow_plasma_fuelling_deuterium)", + physics_variables.f_molflow_plasma_fuelling_deuterium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is tritium", + "(f_molflow_plasma_fuelling_tritium)", + physics_variables.f_molflow_plasma_fuelling_tritium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is helium-3", + "(f_molflow_plasma_fuelling_helium3)", + physics_variables.f_molflow_plasma_fuelling_helium3, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling efficiency", + "(eta_plasma_fuelling)", + physics_variables.eta_plasma_fuelling, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma particles recycled at the LCFS", + "(f_plasma_particles_lcfs_recycled)", + physics_variables.f_plasma_particles_lcfs_recycled, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuel burn-up rate (reactions/s)", + "(fusrat_total)", + physics_variables.fusrat_total, + "OP ", + ) + po.ovarrf( + self.outfile, + "Burn-up fraction", + "(f_plasma_fuel_burnup)", + physics_variables.f_plasma_fuel_burnup, + "OP ", + ) + + if 78 in numerics.icc: + po.osubhd(self.outfile, "Reinke Criterion :") + po.ovarin( + self.outfile, + "index of impurity to be iterated for divertor detachment", + "(impvardiv)", + reinke_variables.impvardiv, + ) + po.ovarre( + self.outfile, + "Minimum Impurity fraction from Reinke", + "(fzmin)", + reinke_variables.fzmin, + "OP ", + ) + po.ovarre( + self.outfile, + "Actual Impurity fraction", + "(fzactual)", + reinke_variables.fzactual, + ) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index c742b3994f..8abb705116 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -3258,98 +3258,7 @@ def outplas(self): if stellarator_variables.istell == 0: self.plasma_bootstrap_current.output() - po.osubhd(self.outfile, "Fuelling :") - po.ovarre( - self.outfile, - "Ratio of He and pellet particle confinement times", - "(tauratio)", - physics_variables.tauratio, - ) - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - physics_variables.molflow_plasma_fuelling_required, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_vv_injected)", - physics_variables.molflow_plasma_fuelling_vv_injected, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of plasma fuelling that is deuterium", - "(f_molflow_plasma_fuelling_deuterium)", - physics_variables.f_molflow_plasma_fuelling_deuterium, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of plasma fuelling that is tritium", - "(f_molflow_plasma_fuelling_tritium)", - physics_variables.f_molflow_plasma_fuelling_tritium, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of plasma fuelling that is helium-3", - "(f_molflow_plasma_fuelling_helium3)", - physics_variables.f_molflow_plasma_fuelling_helium3, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuelling efficiency", - "(eta_plasma_fuelling)", - physics_variables.eta_plasma_fuelling, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of plasma particles recycled at the LCFS", - "(f_plasma_particles_lcfs_recycled)", - physics_variables.f_plasma_particles_lcfs_recycled, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuel burn-up rate (reactions/s)", - "(fusrat_total)", - physics_variables.fusrat_total, - "OP ", - ) - po.ovarrf( - self.outfile, - "Burn-up fraction", - "(f_plasma_fuel_burnup)", - physics_variables.f_plasma_fuel_burnup, - "OP ", - ) - - if 78 in numerics.icc: - po.osubhd(self.outfile, "Reinke Criterion :") - po.ovarin( - self.outfile, - "index of impurity to be iterated for divertor detachment", - "(impvardiv)", - reinke_variables.impvardiv, - ) - po.ovarre( - self.outfile, - "Minimum Impurity fraction from Reinke", - "(fzmin)", - reinke_variables.fzmin, - "OP ", - ) - po.ovarre( - self.outfile, - "Actual Impurity fraction", - "(fzactual)", - reinke_variables.fzactual, - ) + self.fuelling.output_fuelling_info() @staticmethod @nb.njit(cache=True) From 146524699b9beba3e49007faf63bdfbc2c9a2a5a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 14:59:02 +0000 Subject: [PATCH 28/36] :fire: Remove tauratio variable from input variables and related tests for cleanup --- process/core/input.py | 3 --- process/data_structure/physics_variables.py | 5 ----- process/models/physics/fuelling.py | 6 ------ tests/unit/test_physics.py | 4 ---- 4 files changed, 18 deletions(-) diff --git a/process/core/input.py b/process/core/input.py index 476c8c4395..94dfb526d2 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1482,9 +1482,6 @@ def __post_init__(self): "t_plasma_energy_confinement_max": InputVariable( data_structure.physics_variables, float, range=(0.1, 100.0) ), - "tauratio": InputVariable( - data_structure.physics_variables, float, range=(0.1, 100.0) - ), "n_beam_decay_lengths_core_required": InputVariable( data_structure.current_drive_variables, float, range=(0.0, 10.0) ), diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 096e9ea211..9d3a1a04af 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1104,9 +1104,6 @@ f_molflow_plasma_fuelling_helium3: float = None """fraction of plasma fuelling that is helium-3""" -tauratio: float = None -"""tauratio /1.0/ : ratio of He and pellet particle confinement times""" - q95_min: float = None """lower limit for edge safety factor""" @@ -1664,7 +1661,6 @@ def init_physics_variables(): f_molflow_plasma_fuelling_deuterium, \ f_molflow_plasma_fuelling_tritium, \ f_molflow_plasma_fuelling_helium3, \ - tauratio, \ q95_min, \ qstar, \ rad_fraction_sol, \ @@ -1967,7 +1963,6 @@ def init_physics_variables(): f_molflow_plasma_fuelling_deuterium = 0.5 f_molflow_plasma_fuelling_tritium = 0.5 f_molflow_plasma_fuelling_helium3 = 0.0 - tauratio = 1.0 q95_min = 0.0 qstar = 0.0 rad_fraction_sol = 0.8 diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index ea8879e21a..338b565907 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -466,12 +466,6 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): def output_fuelling_info(self): """Output fuelling information to mfile.""" po.oheadr(self.outfile, "Plasma Fuelling") - po.ovarre( - self.outfile, - "Ratio of He and pellet particle confinement times", - "(tauratio)", - physics_variables.tauratio, - ) po.ovarre( self.outfile, "Fuelling rate (nucleus-pairs/s)", diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 3cd17ec6bc..423ccaca12 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2147,7 +2147,6 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - tauratio: Any = None burnup_in: Any = None @@ -2194,7 +2193,6 @@ class PhyauxParam(NamedTuple): "phyauxparam", ( PhyauxParam( - tauratio=1, burnup_in=0, aspect=3, nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19, @@ -2213,7 +2211,6 @@ class PhyauxParam(NamedTuple): expected_t_alpha_confinement=37.993985551650177, ), PhyauxParam( - tauratio=1, burnup_in=0, aspect=3, nd_plasma_fuel_ions_vol_avg=5.8576156204039725e19, @@ -2246,7 +2243,6 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr(physics_variables, "tauratio", phyauxparam.tauratio) monkeypatch.setattr(physics_variables, "burnup_in", phyauxparam.burnup_in) From 4ab2c30c3ecf5e1d7b212a2135840a8bb91c53fa Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 15:11:36 +0000 Subject: [PATCH 29/36] Add tritium and deuterium burnup calculations and update physics variables --- process/data_structure/physics_variables.py | 13 ++- process/models/physics/fuelling.py | 97 ++++++++++++++++++++- process/models/physics/physics.py | 13 +++ tests/unit/test_physics.py | 2 - 4 files changed, 119 insertions(+), 6 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 9d3a1a04af..c0e4f25292 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -254,7 +254,14 @@ f_plasma_fuel_burnup: float = None -"""fractional plasma burnup""" +"""Total fuel burnup fraction in plasma""" + + +f_plasma_tritium_burnup: float = None +"""Tritium burnup fraction in plasma""" + +f_plasma_deuterium_burnup: float = None +"""Deuterium burnup fraction in plasma""" burnup_in: float = None @@ -1480,6 +1487,8 @@ def init_physics_variables(): b_plasma_total, \ e_plasma_magnetic_stored, \ f_plasma_fuel_burnup, \ + f_plasma_tritium_burnup, \ + f_plasma_deuterium_burnup, \ burnup_in, \ b_plasma_vertical_required, \ c_beta, \ @@ -1780,6 +1789,8 @@ def init_physics_variables(): b_plasma_total = 0.0 e_plasma_magnetic_stored = 0.0 f_plasma_fuel_burnup = 0.0 + f_plasma_tritium_burnup = 0.0 + f_plasma_deuterium_burnup = 0.0 burnup_in = 0.0 b_plasma_vertical_required = 0.0 c_beta = 0.5 diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 338b565907..8429b727ee 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -39,9 +39,82 @@ def calculate_fuel_burnup_fraction( ------- float Fuel burnup fraction (dimensionless). + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + + """ + + return 2 * fusrat_total / molflow_plasma_fuelling_vv_injected + + @staticmethod + def calculate_tritium_burnup_fraction( + fusrat_dt_total: float, + molflow_plasma_fuelling_vv_injected: float, + f_molflow_plasma_fuelling_tritium: float, + ) -> float: + """Calculate the tritium burnup fraction + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + f_molflow_plasma_fuelling_tritium : float + Fraction of tritium in the plasma fuelling. + + Returns + ------- + float Tritium burnup fraction (dimensionless). + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + + """ + + return fusrat_dt_total / ( + molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_tritium + ) + + @staticmethod + def calculate_deuterium_burnup_fraction( + fusrat_dt_total: float, + molflow_plasma_fuelling_vv_injected: float, + f_molflow_plasma_fuelling_deuterium: float, + fusrat_plasma_dd_total: float, + fusrat_plasma_dhe3: float, + ) -> float: + """Calculate the deuterium burnup fraction + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + f_molflow_plasma_fuelling_deuterium : float + Fraction of deuterium in the plasma fuelling. + fusrat_plasma_dd_total : float + Total deuterium consumption rate from DD fusion (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + + Returns + ------- + float Deuterium burnup fraction (dimensionless). + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + """ - return fusrat_total / molflow_plasma_fuelling_vv_injected + return (fusrat_dt_total + 2 * fusrat_plasma_dd_total + fusrat_plasma_dhe3) / ( + molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium + ) @staticmethod def calculate_plasma_tritium_flow_rate( @@ -449,8 +522,12 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): f"{mfile.get('f_molflow_plasma_fuelling_deuterium', scan=scan):.4f}\n" f"Fraction of 3-Helium Fuelling: " f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}\n\n" - f"Fuel Burnup Fraction: " + f"Total Fuel Burnup Fraction: " f"{mfile.get('f_plasma_fuel_burnup', scan=scan):.4f}\n" + f"Tritium Burnup Fraction: " + f"{mfile.get('f_plasma_tritium_burnup', scan=scan):.4f}\n" + f"Deuterium Burnup Fraction: " + f"{mfile.get('f_plasma_deuterium_burnup', scan=scan):.4f}" ) fig.text( 0.75, @@ -524,11 +601,25 @@ def output_fuelling_info(self): ) po.ovarrf( self.outfile, - "Burn-up fraction", + "Total fuel burn-up fraction", "(f_plasma_fuel_burnup)", physics_variables.f_plasma_fuel_burnup, "OP ", ) + po.ovarrf( + self.outfile, + "Tritium burn-up fraction", + "(f_plasma_tritium_burnup)", + physics_variables.f_plasma_tritium_burnup, + "OP ", + ) + po.ovarrf( + self.outfile, + "Deuterium burn-up fraction", + "(f_plasma_deuterium_burnup)", + physics_variables.f_plasma_deuterium_burnup, + "OP ", + ) if 78 in numerics.icc: po.osubhd(self.outfile, "Reinke Criterion :") diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 8abb705116..45f01e8582 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -955,6 +955,19 @@ def run(self): fusrat_total=physics_variables.fusrat_total, molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, ) + physics_variables.f_plasma_tritium_burnup = self.fuelling.calculate_tritium_burnup_fraction( + fusrat_dt_total=physics_variables.fusrat_dt_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_tritium=physics_variables.f_molflow_plasma_fuelling_tritium, + ) + + physics_variables.f_plasma_deuterium_burnup = self.fuelling.calculate_deuterium_burnup_fraction( + fusrat_plasma_dd_total=physics_variables.fusrat_plasma_dd_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_deuterium=physics_variables.f_molflow_plasma_fuelling_deuterium, + fusrat_dt_total=physics_variables.fusrat_dt_total, + fusrat_plasma_dhe3=physics_variables.fusrat_plasma_dhe3, + ) physics_variables.ntau, physics_variables.nTtau = ( self.confinement.calculate_double_and_triple_product( diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 423ccaca12..6e1fa12571 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2147,7 +2147,6 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - burnup_in: Any = None aspect: Any = None @@ -2243,7 +2242,6 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr(physics_variables, "burnup_in", phyauxparam.burnup_in) ( From f32c5d2304126dbdce1a5cb1c81a5325b8b0f8a5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 15:37:35 +0000 Subject: [PATCH 30/36] :fire: Refactor plasma fuelling variables: remove unused 'molflow_plasma_fuelling_required' and update references to 'molflow_plasma_fuelling_vv_injected' across models and tests --- process/data_structure/physics_variables.py | 7 +------ process/models/physics/fuelling.py | 7 ------- process/models/vacuum.py | 8 ++++---- tests/unit/test_physics.py | 4 ++-- tests/unit/test_vacuum.py | 2 +- 5 files changed, 8 insertions(+), 20 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index c0e4f25292..b8bbd58ae7 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1090,9 +1090,6 @@ """ -molflow_plasma_fuelling_required: float = None -"""plasma fuelling rate (nucleus-pairs/s)""" - f_plasma_particles_lcfs_recycled: float = None """fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)""" @@ -1100,7 +1097,7 @@ """fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)""" molflow_plasma_fuelling_vv_injected: float = None -"""plasma fuelling rate into the vacuum vessel (nucleus-pairs/s)""" +"""plasma fuelling rate into the vacuum vessel (particles/s)""" f_molflow_plasma_fuelling_deuterium: float = None """fraction of plasma fuelling that is deuterium""" @@ -1663,7 +1660,6 @@ def init_physics_variables(): pden_ion_transport_loss_mw, \ q0, \ q95, \ - molflow_plasma_fuelling_required, \ f_plasma_particles_lcfs_recycled, \ eta_plasma_fuelling, \ molflow_plasma_fuelling_vv_injected, \ @@ -1967,7 +1963,6 @@ def init_physics_variables(): pden_ion_transport_loss_mw = 0.0 q0 = 1.0 q95 = 0.0 - molflow_plasma_fuelling_required = 0.0 f_plasma_particles_lcfs_recycled = 0.0 eta_plasma_fuelling = 0.0 molflow_plasma_fuelling_vv_injected = 1e20 diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 8429b727ee..5ba1b7d56d 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -543,13 +543,6 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): def output_fuelling_info(self): """Output fuelling information to mfile.""" po.oheadr(self.outfile, "Plasma Fuelling") - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - physics_variables.molflow_plasma_fuelling_required, - "OP ", - ) po.ovarre( self.outfile, "Fuelling rate (nucleus-pairs/s)", diff --git a/process/models/vacuum.py b/process/models/vacuum.py index 345153a9a7..71497d71db 100644 --- a/process/models/vacuum.py +++ b/process/models/vacuum.py @@ -54,7 +54,7 @@ def run(self, output: bool = False): # MDK Check this!! gasld = ( 2.0e0 - * physics_variables.molflow_plasma_fuelling_required + * physics_variables.molflow_plasma_fuelling_vv_injected * physics_variables.m_fuel_amu * constants.UMASS ) @@ -123,7 +123,7 @@ def vacuum_simple(self, output) -> float: # One ITER torus cryopump has a throughput of 50 Pa m3/s = 1.2155e+22 molecules/s # Issue #304 n_iter_vacuum_pumps = ( - physics_variables.molflow_plasma_fuelling_required + physics_variables.molflow_plasma_fuelling_vv_injected / vacuum_variables.molflow_vac_pumps ) @@ -167,8 +167,8 @@ def vacuum_simple(self, output) -> float: process_output.ovarre( self.outfile, "Plasma fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - physics_variables.molflow_plasma_fuelling_required, + "(molflow_plasma_fuelling_vv_injected)", + physics_variables.molflow_plasma_fuelling_vv_injected, "OP ", ) process_output.ocmmnt( diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 6e1fa12571..bfc14f2cd2 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2248,7 +2248,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): f_plasma_fuel_burnup, figmer, fusrat, - molflow_plasma_fuelling_required, + molflow_plasma_fuelling_vv_injected, rndfuel, t_alpha_confinement, _, @@ -2270,7 +2270,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): assert fusrat == pytest.approx(phyauxparam.expected_fusrat) - assert molflow_plasma_fuelling_required == pytest.approx( + assert molflow_plasma_fuelling_vv_injected == pytest.approx( phyauxparam.expected_molflow_plasma_fuelling_required ) diff --git a/tests/unit/test_vacuum.py b/tests/unit/test_vacuum.py index a20ef5a1e8..3af9ff46eb 100644 --- a/tests/unit/test_vacuum.py +++ b/tests/unit/test_vacuum.py @@ -42,7 +42,7 @@ def test_simple_model(self, monkeypatch, vacuum): :type tfcoil: tests.unit.test_tfcoil.tfcoil (functional fixture) """ monkeypatch.setattr( - pv, "molflow_plasma_fuelling_required", 7.5745668997694112e22 + pv, "molflow_plasma_fuelling_vv_injected", 7.5745668997694112e22 ) monkeypatch.setattr(pv, "a_plasma_surface", 1500.3146527709359) monkeypatch.setattr(tfv, "n_tf_coils", 18) From 5df343d1ae91448ec3bc2e6578cd7fad33908f22 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 16:01:39 +0000 Subject: [PATCH 31/36] Add Avogadro's number constant with reference documentation --- process/core/constants.py | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/process/core/constants.py b/process/core/constants.py index 100fe4d657..abf35cfc03 100644 --- a/process/core/constants.py +++ b/process/core/constants.py @@ -7,6 +7,13 @@ MFILE = 13 """Machine-optimised output file unit""" + +AVOGADRO_NUMBER = 6.02214076e23 +"""Avogadro's number [1/mol] +Reference: National Institute of Standards and Technology (NIST) +https://physics.nist.gov/cgi-bin/cuu/Value?na|search_for=avogadro +""" + ELECTRON_CHARGE = 1.602176634e-19 """Electron / elementary charge [C] Reference: National Institute of Standards and Technology (NIST) From 5525f22ee76c5a603b0a94170ed65ccfa6867ee4 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 24 Mar 2026 08:58:16 +0000 Subject: [PATCH 32/36] Add plasma fuelling loss calculations and update related variables in the PlasmaFuelling class --- process/data_structure/physics_variables.py | 15 ++++++ process/models/physics/fuelling.py | 59 +++++++++++++++++++++ process/models/physics/physics.py | 18 +------ 3 files changed, 75 insertions(+), 17 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index b8bbd58ae7..5a92fe9013 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -1099,6 +1099,15 @@ molflow_plasma_fuelling_vv_injected: float = None """plasma fuelling rate into the vacuum vessel (particles/s)""" +molflow_plasma_fuelling_vv_injected_moles: float = None +"""plasma fuelling rate into the vacuum vessel (moles/s)""" + +molflow_plasma_fuelling_loss: float = None +"""plasma fuelling rate that dosent make it to plasma (particles/s)""" + +molflow_plasma_fuelling_loss_moles: float = None +"""plasma fuelling rate that dosent make it to plasma (moles/s)""" + f_molflow_plasma_fuelling_deuterium: float = None """fraction of plasma fuelling that is deuterium""" @@ -1663,6 +1672,9 @@ def init_physics_variables(): f_plasma_particles_lcfs_recycled, \ eta_plasma_fuelling, \ molflow_plasma_fuelling_vv_injected, \ + molflow_plasma_fuelling_vv_injected_moles, \ + molflow_plasma_fuelling_loss, \ + molflow_plasma_fuelling_loss_moles, \ f_molflow_plasma_fuelling_deuterium, \ f_molflow_plasma_fuelling_tritium, \ f_molflow_plasma_fuelling_helium3, \ @@ -1966,6 +1978,9 @@ def init_physics_variables(): f_plasma_particles_lcfs_recycled = 0.0 eta_plasma_fuelling = 0.0 molflow_plasma_fuelling_vv_injected = 1e20 + molflow_plasma_fuelling_vv_injected_moles = 0.0 + molflow_plasma_fuelling_loss = 0.0 + molflow_plasma_fuelling_loss_moles = 0.0 f_molflow_plasma_fuelling_deuterium = 0.5 f_molflow_plasma_fuelling_tritium = 0.5 f_molflow_plasma_fuelling_helium3 = 0.0 diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 5ba1b7d56d..23dc9fc6a9 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -22,6 +22,38 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + def run(self): + physics_variables.molflow_plasma_fuelling_vv_injected_moles = ( + physics_variables.molflow_plasma_fuelling_vv_injected + / constants.AVOGADRO_NUMBER + ) + + physics_variables.molflow_plasma_fuelling_loss = ( + physics_variables.molflow_plasma_fuelling_vv_injected + * (1 - physics_variables.eta_plasma_fuelling) + ) + physics_variables.molflow_plasma_fuelling_loss_moles = ( + physics_variables.molflow_plasma_fuelling_loss / constants.AVOGADRO_NUMBER + ) + + physics_variables.f_plasma_fuel_burnup = self.calculate_fuel_burnup_fraction( + fusrat_total=physics_variables.fusrat_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + ) + physics_variables.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction( + fusrat_dt_total=physics_variables.fusrat_dt_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_tritium=physics_variables.f_molflow_plasma_fuelling_tritium, + ) + + physics_variables.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction( + fusrat_plasma_dd_total=physics_variables.fusrat_plasma_dd_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_deuterium=physics_variables.f_molflow_plasma_fuelling_deuterium, + fusrat_dt_total=physics_variables.fusrat_dt_total, + fusrat_plasma_dhe3=physics_variables.fusrat_plasma_dhe3, + ) + @staticmethod def calculate_fuel_burnup_fraction( fusrat_total: float, molflow_plasma_fuelling_vv_injected: float @@ -512,6 +544,12 @@ def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n" f"Total fuelling rate:" f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n" + f"Total fuelling rate: " + f"{mfile.get('molflow_plasma_fuelling_vv_injected_moles', scan=scan):.4e} moles/s\n" + f"Total fuelling loss: " + f"{mfile.get('molflow_plasma_fuelling_loss', scan=scan):.4e} particles/s\n" + f"Total fuelling loss: " + f"{mfile.get('molflow_plasma_fuelling_loss_moles', scan=scan):.4e} moles/s\n" f"Fuelling Rate Efficiency ($\\eta_{{\\text{{fuelling}}}}$): " f"{mfile.get('eta_plasma_fuelling', scan=scan):.4f}\n" f"Recycling Fraction ($R$): " @@ -550,6 +588,27 @@ def output_fuelling_info(self): physics_variables.molflow_plasma_fuelling_vv_injected, "OP ", ) + po.ovarre( + self.outfile, + "Fuelling rate (moles/s)", + "(molflow_plasma_fuelling_vv_injected_moles)", + physics_variables.molflow_plasma_fuelling_vv_injected_moles, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling loss (nucleus-pairs/s)", + "(molflow_plasma_fuelling_loss)", + physics_variables.molflow_plasma_fuelling_loss, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling loss (moles/s)", + "(molflow_plasma_fuelling_loss_moles)", + physics_variables.molflow_plasma_fuelling_loss_moles, + "OP ", + ) po.ovarre( self.outfile, "Fraction of plasma fuelling that is deuterium", diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 45f01e8582..8a5c582a2d 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -951,23 +951,7 @@ def run(self): physics_variables.t_energy_confinement, ) - physics_variables.f_plasma_fuel_burnup = self.fuelling.calculate_fuel_burnup_fraction( - fusrat_total=physics_variables.fusrat_total, - molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, - ) - physics_variables.f_plasma_tritium_burnup = self.fuelling.calculate_tritium_burnup_fraction( - fusrat_dt_total=physics_variables.fusrat_dt_total, - molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, - f_molflow_plasma_fuelling_tritium=physics_variables.f_molflow_plasma_fuelling_tritium, - ) - - physics_variables.f_plasma_deuterium_burnup = self.fuelling.calculate_deuterium_burnup_fraction( - fusrat_plasma_dd_total=physics_variables.fusrat_plasma_dd_total, - molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, - f_molflow_plasma_fuelling_deuterium=physics_variables.f_molflow_plasma_fuelling_deuterium, - fusrat_dt_total=physics_variables.fusrat_dt_total, - fusrat_plasma_dhe3=physics_variables.fusrat_plasma_dhe3, - ) + self.fuelling.run() physics_variables.ntau, physics_variables.nTtau = ( self.confinement.calculate_double_and_triple_product( From 25d44c006cc21dd013faf0b52acff8af2304a2df Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 24 Mar 2026 09:38:10 +0000 Subject: [PATCH 33/36] Add helium-3 flow contour plotting and enhance existing contour plots with improved color mapping and titles --- process/models/physics/fuelling.py | 135 +++++++++++++++++++++++++---- 1 file changed, 117 insertions(+), 18 deletions(-) diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 23dc9fc6a9..35e8d4c6e8 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -327,14 +327,36 @@ def calculate_plasma_helium3_flow_rate( @staticmethod def calculate_plasma_alphas_flow_rate( - fusrat_dt_total, - fusrat_plasma_dhe3, - t_energy_confinement, - f_t_alpha_energy_confinement, - nd_plasma_alphas_vol_avg, - vol_plasma, - ): - """Calculate the alpha particle flow rate in the plasma exhaust.""" + fusrat_dt_total: float, + fusrat_plasma_dhe3: float, + t_energy_confinement: float, + f_t_alpha_energy_confinement: float, + nd_plasma_alphas_vol_avg: float, + vol_plasma: float, + ) -> float: + """Calculate the alpha particle flow rate in the plasma exhaust. + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_t_alpha_energy_confinement : float + Ratio of alpha particle confinement time to energy confinement time (dimensionless). + nd_plasma_alphas_vol_avg : float + Volume-averaged density of alpha particles in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + + Returns + ------- + float + Alpha particle flow rate in the plasma exhaust (particles/s). + + """ # Alpha particle balance @@ -376,8 +398,14 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): ) contour = axis.contourf( - fuelling_range, recycling_range, tritium_flow, levels=15, cmap="RdBu_r" + fuelling_range, + recycling_range, + tritium_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), ) + axis.contour( fuelling_range, recycling_range, @@ -402,11 +430,12 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") axis.set_ylabel("Recycling Fraction [$R$]") - axis.set_title("Plasma Tritium Flow Rate (particles/s)") + axis.set_title("Plasma Tritium Flow Rate (particles/s)", pad=20) axis.minorticks_on() axis.grid(True, which="major", linestyle="-", alpha=0.7) axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") + cbar = plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): """Plot contour of deuterium flow rate vs recycling and fuelling rate.""" @@ -442,7 +471,12 @@ def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int ) contour = axis.contourf( - fuelling_range, recycling_range, deuterium_flow, levels=15, cmap="RdBu_r" + fuelling_range, + recycling_range, + deuterium_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), ) axis.contour( fuelling_range, @@ -467,12 +501,12 @@ def plot_deuterium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int ) axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") - axis.set_ylabel("Recycling Fraction [$R$]") - axis.set_title("Plasma Deuterium Flow Rate (particles/s)") + axis.set_title("Plasma Deuterium Flow Rate (particles/s)", pad=20) axis.minorticks_on() axis.grid(True, which="major", linestyle="-", alpha=0.7) axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") + cbar = plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" @@ -504,7 +538,8 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): fusion_dt_range, alpha_flow, levels=15, - cmap="RdBu_r", + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), ) axis.contour( f_t_alpha_energy_confinement_range, @@ -532,11 +567,75 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" ) axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]") - axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)") + axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)", pad=20) + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + cbar = plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + + def plot_helium3_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): + """Plot contour of helium-3 flow rate vs recycling and fuelling rate.""" + + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + helium3_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + helium3_flow[i, j] = self.calculate_plasma_helium3_flow_rate( + f_molflow_plasma_fuelling_helium3=mfile.get( + "f_molflow_plasma_fuelling_helium3", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_helium3=mfile.get("f_plasma_fuel_helium3", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, + recycling_range, + helium3_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), + ) + axis.contour( + fuelling_range, + recycling_range, + helium3_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + cbar = plt.colorbar(contour, ax=axis, label="Helium-3 Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + axis.set_title("Plasma Helium-3 Flow Rate (particles/s)", pad=20) axis.minorticks_on() axis.grid(True, which="major", linestyle="-", alpha=0.7) axis.grid(True, which="minor", linestyle=":", alpha=0.4) - plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") def plot_fuelling_info(self, fig: plt.Figure, mfile: mf.MFile, scan: int): """Plot fuelling information.""" From 12a62cf755a41be3fdf39ce6ed25d01cfa7f6a5b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 25 Mar 2026 13:29:06 +0000 Subject: [PATCH 34/36] Refactor plotting indices in main_plot function to accommodate additional figures --- process/core/io/plot_proc.py | 125 +++++++++++++++++++---------------- 1 file changed, 67 insertions(+), 58 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 83350b9712..4916a35692 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -13716,16 +13716,33 @@ def main_plot( figs[9].text(0.5, 0.5, msg, ha="center", va="center", wrap=True, fontsize=12) figs[10].text(0.5, 0.5, msg, ha="center", va="center", wrap=True, fontsize=12) - plot_plasma_pressure_profiles(figs[11].add_subplot(222), m_file, scan) - plot_plasma_pressure_gradient_profiles(figs[11].add_subplot(224), m_file, scan) + PlasmaFuelling().plot_tritium_flow_contour( + axis=figs[11].add_subplot(231), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_deuterium_flow_contour( + axis=figs[11].add_subplot(232), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_helium3_flow_contour( + axis=figs[11].add_subplot(233), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_alpha_flow_contour( + axis=figs[11].add_subplot(223), mfile=m_file, scan=scan + ) + + PlasmaFuelling().plot_fuelling_info(figs[11], m_file, scan) + + figs[11].subplots_adjust(wspace=0.3, hspace=0.35) + + plot_plasma_pressure_profiles(figs[12].add_subplot(222), m_file, scan) + plot_plasma_pressure_gradient_profiles(figs[12].add_subplot(224), m_file, scan) # Currently only works with Sauter geometry as plasma has a closed surface if i_shape == 1: plot_plasma_poloidal_pressure_contours( - figs[11].add_subplot(121, aspect="equal"), m_file, scan + figs[12].add_subplot(121, aspect="equal"), m_file, scan ) else: - ax = figs[11].add_subplot(131, aspect="equal") + ax = figs[12].add_subplot(131, aspect="equal") msg = ( "Plasma poloidal pressure contours require a closed (Sauter) plasma boundary " "(i_plasma_shape == 1). " @@ -13744,30 +13761,30 @@ def main_plot( ) ax.axis("off") - plot_magnetic_fields_in_plasma(figs[12].add_subplot(122), m_file, scan) - plot_beta_profiles(figs[12].add_subplot(221), m_file, scan) + plot_magnetic_fields_in_plasma(figs[13].add_subplot(122), m_file, scan) + plot_beta_profiles(figs[13].add_subplot(221), m_file, scan) - plot_ebw_ecrh_coupling_graph(figs[13].add_subplot(111), m_file, scan) + plot_ebw_ecrh_coupling_graph(figs[14].add_subplot(111), m_file, scan) - plot_bootstrap_comparison(figs[14].add_subplot(221), m_file, scan) - PlasmaCurrent.plot_plasma_current_comparison(figs[14].add_subplot(224), m_file, scan) - plot_h_threshold_comparison(figs[15].add_subplot(224), m_file, scan) - plot_density_limit_comparison(figs[15].add_subplot(221), m_file, scan) - plot_confinement_time_comparison(figs[16].add_subplot(224), m_file, scan) + plot_bootstrap_comparison(figs[15].add_subplot(221), m_file, scan) + PlasmaCurrent.plot_plasma_current_comparison(figs[15].add_subplot(224), m_file, scan) + plot_h_threshold_comparison(figs[16].add_subplot(224), m_file, scan) + plot_density_limit_comparison(figs[16].add_subplot(221), m_file, scan) + plot_confinement_time_comparison(figs[17].add_subplot(224), m_file, scan) - plot_debye_length_profile(figs[17].add_subplot(232), m_file, scan) - plot_velocity_profile(figs[17].add_subplot(233), m_file, scan) - plot_plasma_coloumb_logarithms(figs[17].add_subplot(231), m_file, scan) + plot_debye_length_profile(figs[18].add_subplot(232), m_file, scan) + plot_velocity_profile(figs[18].add_subplot(233), m_file, scan) + plot_plasma_coloumb_logarithms(figs[18].add_subplot(231), m_file, scan) - plot_electron_frequency_profile(figs[17].add_subplot(212), m_file, scan) + plot_electron_frequency_profile(figs[18].add_subplot(212), m_file, scan) - plot_ion_frequency_profile(figs[18].add_subplot(311), m_file, scan) + plot_ion_frequency_profile(figs[19].add_subplot(311), m_file, scan) - plot_larmor_radius_profile(figs[18].add_subplot(313), m_file, scan) + plot_larmor_radius_profile(figs[19].add_subplot(313), m_file, scan) # Plot poloidal cross-section poloidal_cross_section( - figs[19].add_subplot(121, aspect="equal"), + figs[20].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -13777,7 +13794,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[19].add_subplot(122, aspect="equal"), + figs[20].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -13785,27 +13802,27 @@ def main_plot( ) # Plot color key - ax17 = figs[19].add_subplot(222) + ax17 = figs[20].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) plot_full_machine_poloidal_cross_section( - figs[20].add_subplot(111, aspect="equal"), + figs[21].add_subplot(111, aspect="equal"), m_file, scan, radial_build, colour_scheme, ) - ax18 = figs[21].add_subplot(211) + ax18 = figs[22].add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) plot_radial_build(ax18, m_file, colour_scheme) # Make each axes smaller vertically to leave room for the legend - ax185 = figs[22].add_subplot(211) + ax185 = figs[23].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[22].add_subplot(212) + ax18b = figs[23].add_subplot(212) ax18b.set_position([0.1, 0.13, 0.8, 0.32]) plot_upper_vertical_build(ax185, m_file, colour_scheme) plot_lower_vertical_build(ax18b, m_file, colour_scheme) @@ -13813,70 +13830,62 @@ def main_plot( # Can only plot WP and turn structure if superconducting coil at the moment if m_file.get("i_tf_sup", scan=scan) == 1: # TF coil with WP - ax19 = figs[23].add_subplot(221, aspect="equal") + ax19 = figs[24].add_subplot(221, aspect="equal") ax19.set_position([ 0.025, 0.45, 0.5, 0.5, ]) # Half height, a bit wider, top left - plot_superconducting_tf_wp(ax19, m_file, scan, figs[23]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[24]) # TF coil turn structure - ax20 = figs[24].add_subplot(325, aspect="equal") + ax20 = figs[25].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[24], m_file, scan) - plot_205 = figs[24].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[25], m_file, scan) + plot_205 = figs[25].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) - plot_cable_in_conduit_cable(plot_205, figs[24], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[25], m_file, scan) else: - ax19 = figs[23].add_subplot(211, aspect="equal") + ax19 = figs[24].add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) - plot_resistive_tf_wp(ax19, m_file, scan, figs[23]) - plot_resistive_tf_info(ax19, m_file, scan, figs[23]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[24]) + plot_resistive_tf_info(ax19, m_file, scan, figs[24]) plot_tf_coil_structure( - figs[25].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[26].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[26], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[27], m_file, scan) - plot_tf_stress(figs[27].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[28].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[28].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[29].add_subplot(121, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan ) plot_cs_turn_structure( - figs[29].add_subplot(224, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(224, aspect="equal"), figs[30], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[30].add_subplot(221, aspect="equal"), m_file, scan + figs[31].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[30].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[30].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[31].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[31].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[31], m_file, scan) - ax_blanket = figs[31].add_subplot(122, aspect="equal") - plot_blkt_structure(ax_blanket, figs[31], m_file, scan, radial_build, colour_scheme) + plot_blkt_pipe_bends(figs[32], m_file, scan) + ax_blanket = figs[32].add_subplot(122, aspect="equal") + plot_blkt_structure(ax_blanket, figs[32], m_file, scan, radial_build, colour_scheme) plot_main_power_flow( - figs[32].add_subplot(111, aspect="equal"), m_file, scan, figs[32] + figs[33].add_subplot(111, aspect="equal"), m_file, scan, figs[33] ) - ax24 = figs[33].add_subplot(111) + ax24 = figs[34].add_subplot(111) # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height ax24.set_position([0.08, 0.35, 0.84, 0.57]) - plot_system_power_profiles_over_time(ax24, m_file, scan, figs[33]) - - - ax25 = figs[31].add_subplot(221) - PlasmaExhaust().plot_tritium_flow_contour(ax25, m_file, scan) - ax26 = figs[31].add_subplot(222) - PlasmaExhaust().plot_deuterium_flow_contour(ax26, m_file, scan) - ax27 = figs[31].add_subplot(223) - PlasmaExhaust().plot_alpha_flow_contour(ax27, m_file, scan) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[34]) def create_thickness_builds(m_file, scan: int): @@ -13953,7 +13962,7 @@ def main(args=None): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(34)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(35)] # run main_plot main_plot( From 7687c7bff9787620f119e896c12cd21ca076c0c9 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 27 Mar 2026 09:07:46 +0000 Subject: [PATCH 35/36] Enhance plasma fuelling documentation: add definitions for exhaust efficiency and recycling coefficient, and introduce METIS Alpha Confinement model --- .../source/physics-models/plasma_fuelling.md | 27 ++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md index ededb6966c..86dfd90ea4 100644 --- a/documentation/source/physics-models/plasma_fuelling.md +++ b/documentation/source/physics-models/plasma_fuelling.md @@ -35,12 +35,31 @@ The fuelling fractional compositions is given by $f$ $\tau_{\text{fuel}}^* = (\tau_p) / (1-R)$ +The (effective) exhaust efficiency is is given by , $\eta_{\text{eff}} = 1- R$ + +The factor $\frac{R}{1-R}$ is the mean number of recycling events back into the burning region experieced by a particle before it is pumped away. + +The definition of the recycling coefficient $R = 1- \frac{\Gamma_{\text{pumps}}}{\Gamma_{\text{out}}}$, where $\Gamma_{\text{pumps}}$ is the number of particles exhausted by the pumps per second and $\Gamma_{\text{out}}$ is the number of particles per second transported radially outwards across the separatrix. + + Where $\tau_p$ is the particle confinement time which we can assume is approximately equal to the energy confinement time ($\tau_p = \tau_E$). -The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1]. +!!! note "Quantifying $R$" + + The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1]. + + +## METIS Alpha Confinement + +$$ +\tau_{\alpha} = f_{\alpha}\tau_{\text{E}}\frac{R}{1-R}\tau_{\text{ne}} +$$ + +This is the model currently in METIS[^2] and is found in [^3] + -------------- @@ -126,4 +145,10 @@ $$ [^1]: G. L. Jackson, V. S. Chan, and R. D. Stambaugh, “An Analytic Expression for the Tritium Burnup Fraction in Burning-Plasma Devices,” Fusion Science and Technology, vol. 64, no. 1, pp. 8–12, Jul. 2013, doi: https://doi.org/10.13182/fst13-a17042. + +[^2]: J. F. Artaud et al., “Metis: a fast integrated tokamak modelling tool for scenario design,” Nuclear Fusion, vol. 58, no. 10, pp. 105001–105001, Aug. 2018, doi: https://doi.org/10.1088/1741-4326/aad5b1. + +[^3]: D. Reiter, H. Kever, G. H. Wolf, M. Baelmans, R. Behrisch, and R. Schneider, “Helium removal from tokamaks,” Plasma Physics and Controlled Fusion, vol. 33, no. 13, pp. 1579–1600, Nov. 1991, doi: https://doi.org/10.1088/0741-3335/33/13/008. +‌ +‌ ‌ \ No newline at end of file From cf525c5a6a8d2db850d7ec402960daf40d37a81f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 31 Mar 2026 10:39:05 +0100 Subject: [PATCH 36/36] Refactor Physics model and unit tests to integrate PlasmaFuelling class and update parameter descriptions --- process/models/physics/physics.py | 18 ++----- tests/unit/test_physics.py | 87 ++----------------------------- tests/unit/test_stellarator.py | 2 + 3 files changed, 9 insertions(+), 98 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 8a5c582a2d..3d4f6f0448 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1436,24 +1436,12 @@ def phyaux( Parameters ---------- - aspect : float - Plasma aspect ratio. - nd_plasma_fuel_ions_vol_avg : float - Fuel ion density (/m3). - fusden_total : float - Fusion reaction rate from plasma and beams (/m3/s). fusden_alpha_total : float - Alpha particle production rate (/m3/s). - plasma_current : float - Plasma current (A). - sbar : float - Exponent for aspect ratio (normally 1). + Total alpha particle production rate density (m^-3 s^-1). nd_plasma_alphas_vol_avg : float - Alpha ash density (/m3). + Volume averaged alpha particle density (m^-3). t_energy_confinement : float - Global energy confinement time (s). - vol_plasma : float - Plasma volume (m3). + Energy confinement time (s). Returns ------- diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index bfc14f2cd2..b48959fe5f 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -23,6 +23,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( @@ -66,6 +67,7 @@ def physics(): PlasmaConfinementTime(), PlasmaConfinementTransition(), PlasmaCurrent(), + plasma_fuelling=PlasmaFuelling(), ) @@ -2147,44 +2149,12 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - burnup_in: Any = None - - aspect: Any = None - - nd_plasma_electrons_vol_avg: Any = None - - te: Any = None - - nd_plasma_fuel_ions_vol_avg: Any = None - nd_plasma_alphas_vol_avg: Any = None - fusden_total: Any = None - fusden_alpha_total: Any = None - plasma_current: Any = None - - sbar: Any = None - t_energy_confinement: Any = None - vol_plasma: Any = None - - expected_burnup: Any = None - - expected_ntau: Any = None - - expected_nTtau: Any = None - - expected_figmer: Any = None - - expected_fusrat: Any = None - - expected_molflow_plasma_fuelling_required: Any = None - - expected_rndfuel: Any = None - expected_t_alpha_confinement: Any = None @@ -2192,39 +2162,15 @@ class PhyauxParam(NamedTuple): "phyauxparam", ( PhyauxParam( - burnup_in=0, - aspect=3, - nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19, - nd_plasma_alphas_vol_avg=7.5e18, - fusden_total=1.9852091609123786e17, fusden_alpha_total=1.973996644759543e17, - plasma_current=18398455.678867526, - sbar=1, + nd_plasma_alphas_vol_avg=7.5e18, t_energy_confinement=3.401323521525641, - vol_plasma=1888.1711539956691, - expected_burnup=0.20383432558954095, - expected_figmer=55.195367036602576, - expected_fusrat=3.7484146722826997e20, - expected_molflow_plasma_fuelling_required=1.8389516394951276e21, - expected_rndfuel=3.7484146722826997e20, expected_t_alpha_confinement=37.993985551650177, ), PhyauxParam( - burnup_in=0, - aspect=3, - nd_plasma_fuel_ions_vol_avg=5.8576156204039725e19, - nd_plasma_alphas_vol_avg=7.5e18, - fusden_total=1.9843269653375773e17, fusden_alpha_total=1.9731194318497056e17, - plasma_current=18398455.678867526, - sbar=1, + nd_plasma_alphas_vol_avg=7.5e18, t_energy_confinement=3.402116961408892, - vol_plasma=1888.1711539956691, - expected_burnup=0.20387039462081086, - expected_figmer=55.195367036602576, - expected_fusrat=3.7467489360461772e20, - expected_molflow_plasma_fuelling_required=1.8378092331723546e21, - expected_rndfuel=3.7467489360461772e20, expected_t_alpha_confinement=38.010876984618747, ), ), @@ -2242,40 +2188,15 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr(physics_variables, "burnup_in", phyauxparam.burnup_in) - ( - f_plasma_fuel_burnup, - figmer, - fusrat, - molflow_plasma_fuelling_vv_injected, - rndfuel, t_alpha_confinement, _, ) = physics.phyaux( - aspect=phyauxparam.aspect, - nd_plasma_fuel_ions_vol_avg=phyauxparam.nd_plasma_fuel_ions_vol_avg, nd_plasma_alphas_vol_avg=phyauxparam.nd_plasma_alphas_vol_avg, - fusden_total=phyauxparam.fusden_total, fusden_alpha_total=phyauxparam.fusden_alpha_total, - plasma_current=phyauxparam.plasma_current, - sbar=phyauxparam.sbar, t_energy_confinement=phyauxparam.t_energy_confinement, - vol_plasma=phyauxparam.vol_plasma, - ) - - assert f_plasma_fuel_burnup == pytest.approx(phyauxparam.expected_burnup) - - assert figmer == pytest.approx(phyauxparam.expected_figmer) - - assert fusrat == pytest.approx(phyauxparam.expected_fusrat) - - assert molflow_plasma_fuelling_vv_injected == pytest.approx( - phyauxparam.expected_molflow_plasma_fuelling_required ) - assert rndfuel == pytest.approx(phyauxparam.expected_rndfuel) - assert t_alpha_confinement == pytest.approx(phyauxparam.expected_t_alpha_confinement) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 2c5dfe5951..11b7a49bfd 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -33,6 +33,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( Physics, @@ -98,6 +99,7 @@ def stellarator(): PlasmaConfinementTime(), PlasmaConfinementTransition(), PlasmaCurrent(), + plasma_fuelling=PlasmaFuelling(), ), Neoclassics(), plasma_beta=PlasmaBeta(),