From 73934cabf61d12c8d42b6616404230e84189a528 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 10:43:25 +0000 Subject: [PATCH 01/37] Add new plasma physics variables and update iteration variable count --- process/core/solver/iteration_variables.py | 12 +++ process/data_structure/numerics.py | 97 ++++++++++++++++++++- process/data_structure/physics_variables.py | 8 ++ 3 files changed, 116 insertions(+), 1 deletion(-) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 42abd1dd32..2593cae3dd 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -249,6 +249,18 @@ class IterationVariable: 174: IterationVariable("triang", "physics", 0.00, 1.00), 175: IterationVariable("kappa", "physics", 0.00, 10.00), 176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30), + 177: IterationVariable( + "f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0 + ), + 178: IterationVariable( + "eta_plasma_fuelling", "physics", 0.01, 1.0 + ), + 179: IterationVariable( + "molflow_plasma_fuelling_vv_injected", + "physics", + 0.0, + 1e24, + ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 8499302c80..03a2fda8b8 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -99,7 +99,102 @@ def description(self): return self._description_ -ipnvars: int = 177 +class PROCESSRunMode(IntEnum): + """Enumeration of the available PROCESS run modes, which determine the behaviour + of the code in various places. This is controlled by the `ioptimz` variable + """ + + EVALUATION = (-2, "Evaluation mode (no optimisation)") + """In this mode, the code will not perform any optimisation, and will instead + simply evaluate the constraints for the given input parameters, which is useful + for testing and for evaluating the performance of a given design point without + trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd + root-finding method from `scipy.optimize`) to seek a *consistent* solution by + varying a subset of the iteration variables until the consistency constraints + (equality constraints whose residuals must be driven to zero) are simultaneously + satisfied; no figure-of-merit is optimised, and the solver simply tries to find + a root of the constraint-residual vector. + """ + OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") + """In this mode, the code will perform optimisation using the VMCON solver + (or a custom solver if specified) to try to find a design point that optimises + the figure of merit while satisfying the constraints. This is the default mode + of operation for PROCESS. + """ + + def __new__(cls, value: int, description: str): + """Create a new PROCESSRunMode enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this run mode. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this run mode.""" + return self._description_ + + +class FiguresOfMerit(IntEnum): + """Enumeration of the available figures of merit (FoM) that can be used as + objective functions for optimisation in PROCESS. + """ + + MAJOR_RADIUS = (1, "Plasma major radius (R₀)") + NEUTRON_WALL_LOAD = (3, "Neutron wall load") + P_TF_PLUS_P_PF = (4, "TF & PF coil power") + FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") + COST_OF_ELECTRICITY = (6, "Cost of electricity") + CAPITAL_COST = (7, "Plant capital cost") + ASPECT_RATIO = (8, "Plasma aspect ratio") + DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") + TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") + TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") + PULSE_LENGTH = (14, "Pulse length") + PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") + MIN_R0_MAX_TAU_BURN = ( + 16, + "Linear combination of major radius (minimised) and pulse length (maximised)", + ) + NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") + NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") + MAX_Q_MAX_T_PLANT_PULSE_BURN = ( + 19, + "Linear combination of big Q and pulse length (maximised)", + ) + + def __new__(cls, value: int, description: str): + """Create a new FiguresOfMerit enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this figure of merit. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this figure of merit.""" + return self._description_ + + +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 6ded110158..a5a3f313e8 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -875,6 +875,14 @@ class PhysicsData: molflow_plasma_fuelling_required: float = 0.0 """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 = 1.0 """tauratio /1.0/ : ratio of He and pellet particle confinement times""" From 08959018501061563ff87b980df3780860da41e1 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 10:59:44 +0000 Subject: [PATCH 02/37] Add particle balance constraint equation and update constraint count --- process/core/solver/constraints.py | 73 ++++++++++++++++++++++++++++++ process/data_structure/numerics.py | 2 +- 2 files changed, 74 insertions(+), 1 deletion(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 6d7990a395..71da7ffb58 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1785,6 +1785,79 @@ def constraint_equation_92(constraint_registration, data): ) +@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) + + +@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, data: DataStructure): """Evaluates the constraints given the current state of PROCESS. diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 03a2fda8b8..3bf34f857c 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -197,7 +197,7 @@ def description(self): 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 e4b3864b7b65bb0177ed4ca75148a8f93c0b4a19 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 11:03:21 +0000 Subject: [PATCH 03/37] 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 782e5fa500..c69c835ce8 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -696,6 +696,15 @@ def __post_init__(self): "f_volflow_vac_pumps_impedance": InputVariable("vacuum", float, range=(1e-06, 1.0)), "volflow_vac_pumps_max": InputVariable("vacuum", float, range=(1e-06, 1000.0)), "molflow_vac_pumps": InputVariable("vacuum", 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( "heat_transport", float, range=(0.0, 1000.0) ), From 1513adb8cddfbb18f0a8549f7bc1475977c34650 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 11:12:39 +0000 Subject: [PATCH 04/37] 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 71da7ffb58..dba7e1ef5d 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1786,7 +1786,7 @@ def constraint_equation_92(constraint_registration, data): @ConstraintManager.register_constraint(93, "particles/s", "=") -def constraint_equation_93(): +def constraint_equation_93(constraint_registration): """ Particle balance """ @@ -1816,9 +1816,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) @ConstraintManager.register_constraint(93, "particles/s", "=") diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index 2593cae3dd..ded7806285 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -258,7 +258,7 @@ class IterationVariable: 179: IterationVariable( "molflow_plasma_fuelling_vv_injected", "physics", - 0.0, + 1.0, 1e24, ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 3bf34f857c..96ce1271fd 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -355,6 +355,7 @@ def description(self): * (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 @@ -481,6 +482,7 @@ def description(self): * (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 @@ -763,6 +765,7 @@ def init_numerics(): "CS achievable stress load cycles lower limit ", "ECRH ignitability ", # Stellarator constraint "Fuel composition consistency ", + "Particle balance consistency ", ] ixc = np.array([0] * ipnvars) From 422b518273a318dbda8fb184476b920f248fd424 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 14:56:21 +0000 Subject: [PATCH 05/37] 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 dba7e1ef5d..5983a96e59 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1785,10 +1785,10 @@ def constraint_equation_92(constraint_registration, data): ) -@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 @@ -1810,13 +1810,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) @ConstraintManager.register_constraint(93, "particles/s", "=") diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 96ce1271fd..b1c29f89ab 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -197,7 +197,7 @@ def description(self): ipnvars: int = 180 """total number of variables available for iteration""" -ipeqns: int = 93 +ipeqns: int = 95 """number of constraint equations available""" ipnfoms: int = 19 @@ -765,7 +765,9 @@ def init_numerics(): "CS achievable stress load cycles lower limit ", "ECRH ignitability ", # Stellarator constraint "Fuel composition consistency ", - "Particle balance consistency ", + "Tritium particle balance consistency ", + "Deuterium particle balance consistency ", + "Alpha particle balance consistency ", ] ixc = np.array([0] * ipnvars) From 79ead71add3498448a0bbaa88dc1891140820c87 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 15:41:44 +0000 Subject: [PATCH 06/37] 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 5983a96e59..ddb9eb2002 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1785,7 +1785,7 @@ def constraint_equation_92(constraint_registration, data): ) -@ConstraintManager.register_constraint(93, "particles/s", "<=") +@ConstraintManager.register_constraint(93, "particles/s", "=") def constraint_equation_93(constraint_registration): """ Tritium particle balance @@ -1810,10 +1810,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 @@ -1835,10 +1835,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 @@ -1848,9 +1848,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) @ConstraintManager.register_constraint(93, "particles/s", "=") From 8b17392f6c681116873b3a018a2c2a015e624a52 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 17 Mar 2026 16:49:07 +0000 Subject: [PATCH 07/37] Add new total D-T fusion rate variables --- process/core/io/plot/summary.py | 10 ++++++---- process/core/solver/constraints.py | 6 +++--- process/data_structure/physics_variables.py | 6 ++++++ process/models/physics/physics.py | 21 +++++++++++++++++++++ 4 files changed, 36 insertions(+), 7 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index ba8823a32a..52d0fa99b4 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -11666,9 +11666,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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( @@ -11688,7 +11690,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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 ddb9eb2002..bde8fadd82 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1797,7 +1797,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 @@ -1822,7 +1822,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 @@ -1845,7 +1845,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 a5a3f313e8..539230c1e3 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -366,6 +366,12 @@ class PhysicsData: fusrat_total: float = 0.0 """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] = field(default_factory=list) """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 573ab73acf..66370e64eb 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -645,6 +645,13 @@ def run(self): self.data.physics.fusrat_total = ( self.data.physics.fusden_total * self.data.physics.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 ( @@ -1790,6 +1797,20 @@ def outplas(self): self.data.physics.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/m³/sec)", From 547ec90dc8f9d3f137904b22cb32f174b8801bfd Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 18 Mar 2026 16:30:05 +0000 Subject: [PATCH 08/37] 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 24e170f451..40293ba3fc 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -621,6 +621,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 bde8fadd82..00a419dd25 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1798,19 +1798,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", "=") @@ -1823,19 +1821,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", "=") @@ -1845,7 +1841,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 ) / ( @@ -1853,7 +1850,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) @ConstraintManager.register_constraint(93, "particles/s", "=") From 05fdb3a0fac1b17f0c16bf185b4d8585ccac4e31 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:32:06 +0000 Subject: [PATCH 09/37] Update molflow plasma fuelling variable ranges for consistency --- process/core/input.py | 2 +- process/core/solver/iteration_variables.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/process/core/input.py b/process/core/input.py index c69c835ce8..ad4503f211 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -703,7 +703,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( "heat_transport", float, range=(0.0, 1000.0) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index ded7806285..efc464ec20 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -258,7 +258,7 @@ class IterationVariable: 179: IterationVariable( "molflow_plasma_fuelling_vv_injected", "physics", - 1.0, + 1e18, 1e24, ), } From 50550a25831e84d9546f9e6f1feeb464db157933 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:39:04 +0000 Subject: [PATCH 10/37] Refine fusion power output formatting for improved precision in plot displays --- process/core/io/plot/summary.py | 40 ++++++++++++++++----------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 52d0fa99b4..63f53f8ff7 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -2927,7 +2927,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" ) @@ -11643,8 +11643,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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( @@ -11690,7 +11690,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.275, + 0.285, 0.8, "$\\text{D - T}$", fontsize=20, @@ -11701,7 +11701,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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" ) @@ -11732,7 +11732,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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, @@ -11762,15 +11762,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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( @@ -11801,11 +11801,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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 494e85b6e7a0435fdb32802d7c7a79bb09f0b753 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 14:42:10 +0000 Subject: [PATCH 11/37] Add D-D fusion reaction rate variables for helium and tritium branches --- process/data_structure/physics_variables.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 539230c1e3..af544f8af2 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -366,12 +366,25 @@ class PhysicsData: fusrat_total: float = 0.0 """fusion reaction rate, from beams and plasma (reactions/sec)""" -fusrat_plasma_dt: float = None -""" D-T fusion reaction rate in 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_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_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] = field(default_factory=list) """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" From 17b992366f85706e950e0cc4fed7f2d7ba0a038d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:08:30 +0000 Subject: [PATCH 12/37] Add D-D fusion rate calculations for tritium and helium branches --- process/core/io/plot/summary.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/summary.py b/process/core/io/plot/summary.py index 63f53f8ff7..751b6a9c0a 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -11702,7 +11702,9 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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( @@ -11722,8 +11724,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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 c81e069e50..96165c0ed7 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -451,6 +451,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 @@ -551,6 +555,10 @@ def dd_triton_reaction(self): # Proton production rate [particles/m³/s] proton_rate_density = fusion_rate_density + 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 66370e64eb..1a43156124 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1811,6 +1811,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/m³/sec)", From b0d315e38aa97afbd2440adfc67d18fcee6ea1c2 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:53:07 +0000 Subject: [PATCH 13/37] Add plasma tritium flow rate calculations and contour plotting --- process/core/io/plot/summary.py | 6 ++ process/models/physics/exhaust.py | 96 +++++++++++++++++++++++++++++++ process/models/physics/physics.py | 7 +++ 3 files changed, 109 insertions(+) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 751b6a9c0a..f0ad3bddfe 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -58,6 +58,7 @@ ElectronCyclotron, ) from process.models.physics.density_limit import DensityLimitModel +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.physics import ( @@ -15315,6 +15316,11 @@ 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[35]) + + + ax25 = figs[36].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 2b3a30e321..2d5d01f321 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -2,6 +2,10 @@ import logging +import matplotlib.pyplot as plt +import numpy as np + +import process.core.io.mfile as mf from process.core import constants from process.core import process_output as po from process.core.model import Model @@ -190,3 +194,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 1a43156124..a462f7ed70 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2367,6 +2367,13 @@ def outplas(self): self.data.physics.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 2bdb30c39bf62b6ec30238bc613e5f9466dd9a3b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 15:57:34 +0000 Subject: [PATCH 14/37] Add total D-D fusion rate calculation and update output formatting --- process/core/io/plot/summary.py | 3 ++- process/data_structure/physics_variables.py | 3 +++ process/models/physics/physics.py | 18 ++++++++++++++++++ 3 files changed, 23 insertions(+), 1 deletion(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index f0ad3bddfe..156a2e3e6f 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -11705,7 +11705,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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 af544f8af2..ae042e81e7 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -378,6 +378,9 @@ class PhysicsData: 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_dd_helion: float = None """D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index a462f7ed70..ad38f7d338 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -648,6 +648,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) @@ -1826,6 +1830,20 @@ 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, + "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/m³/sec)", From eeeca04a2672c5b7ec5e3d93264bd390ab233dfd Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 16:35:54 +0000 Subject: [PATCH 15/37] Add deuterium and alpha particle flow rate calculations and contour plotting --- process/core/io/plot/summary.py | 4 + process/models/physics/exhaust.py | 142 +++++++++++++++++++++++++++--- 2 files changed, 134 insertions(+), 12 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 156a2e3e6f..d06f1fe1cd 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -15321,6 +15321,10 @@ def main_plot( ax25 = figs[36].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 2d5d01f321..976b11146b 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -200,7 +200,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, @@ -208,11 +208,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)) @@ -224,25 +224,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.""" @@ -258,7 +276,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( @@ -286,3 +306,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 6ac4a3b3f7fc71b13896e965254a547039fc4b67 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 20 Mar 2026 17:23:48 +0000 Subject: [PATCH 16/37] 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 976b11146b..b1ad8b846e 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -299,6 +299,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)") @@ -347,6 +361,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)") @@ -359,7 +387,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), @@ -395,6 +423,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 ad38f7d338..afcb57fe52 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2392,6 +2392,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 bb233cd9ad6a748ef47779b04806910d9c0c470e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 15:08:52 +0000 Subject: [PATCH 17/37] 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 42595c58bf..ec5ec1d975 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -53,6 +53,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 6391258f336be4536eb99c0bb09c4f41d279a9e1 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 15:35:12 +0000 Subject: [PATCH 18/37] Add D-3He fusion rate and neutron production rate calculations --- process/core/io/plot/summary.py | 22 ++++++++------ process/data_structure/physics_variables.py | 13 +++++++++ process/models/physics/fusion_reactions.py | 4 +++ process/models/physics/physics.py | 32 +++++++++++++++++++++ 4 files changed, 62 insertions(+), 9 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index d06f1fe1cd..d0c8c8857d 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -11703,7 +11703,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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" @@ -11711,7 +11711,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): axis.text( 0.05, - 0.65, + 0.625, textstr_dd, fontsize=9, verticalalignment="bottom", @@ -11736,11 +11736,14 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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", @@ -11755,8 +11758,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.21, - 0.59, + 0.285, + 0.56, "$\\text{D - 3He}$", fontsize=20, verticalalignment="top", @@ -11809,7 +11812,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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( @@ -11829,8 +11833,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: 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 ae042e81e7..e3cee1d7d0 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -381,6 +381,19 @@ class PhysicsData: 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_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_dd_helion: float = None """D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 96165c0ed7..19512707b1 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -347,6 +347,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 afcb57fe52..694e3e2b1b 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -653,6 +653,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) ) @@ -1837,6 +1841,34 @@ 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, + "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, "D-D Fusion rate: total (reactions/sec)", From d111b06ed3702156cad9d4e75d72529676ebb395 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 16:42:25 +0000 Subject: [PATCH 19/37] 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 | 53 ++++- process/core/solver/iteration_variables.py | 18 ++ process/data_structure/numerics.py | 197 +++++++++++++++++- process/data_structure/physics_variables.py | 9 + process/models/physics/physics.py | 21 ++ 7 files changed, 303 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 ad4503f211..989f5ca256 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -705,6 +705,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( "heat_transport", float, range=(0.0, 1000.0) ), diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 00a419dd25..6de2cd2add 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1797,8 +1797,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 @@ -1806,7 +1806,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) @@ -1820,8 +1820,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 @@ -1841,18 +1841,55 @@ 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, + ) + + +@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, + ) + + @ConstraintManager.register_constraint(93, "particles/s", "=") def constraint_equation_93(): """ diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index efc464ec20..1f20501742 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -261,6 +261,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 b1c29f89ab..a7cfb35b26 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -194,10 +194,202 @@ def description(self): return self._description_ -ipnvars: int = 180 + +class PROCESSRunMode(IntEnum): + """Enumeration of the available PROCESS run modes, which determine the behaviour + of the code in various places. This is controlled by the `ioptimz` variable + """ + + EVALUATION = (-2, "Evaluation mode (no optimisation)") + """In this mode, the code will not perform any optimisation, and will instead + simply evaluate the constraints for the given input parameters, which is useful + for testing and for evaluating the performance of a given design point without + trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd + root-finding method from `scipy.optimize`) to seek a *consistent* solution by + varying a subset of the iteration variables until the consistency constraints + (equality constraints whose residuals must be driven to zero) are simultaneously + satisfied; no figure-of-merit is optimised, and the solver simply tries to find + a root of the constraint-residual vector. + """ + OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") + """In this mode, the code will perform optimisation using the VMCON solver + (or a custom solver if specified) to try to find a design point that optimises + the figure of merit while satisfying the constraints. This is the default mode + of operation for PROCESS. + """ + + def __new__(cls, value: int, description: str): + """Create a new PROCESSRunMode enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this run mode. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this run mode.""" + return self._description_ + + +class FiguresOfMerit(IntEnum): + """Enumeration of the available figures of merit (FoM) that can be used as + objective functions for optimisation in PROCESS. + """ + + MAJOR_RADIUS = (1, "Plasma major radius (R₀)") + NEUTRON_WALL_LOAD = (3, "Neutron wall load") + P_TF_PLUS_P_PF = (4, "TF & PF coil power") + FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") + COST_OF_ELECTRICITY = (6, "Cost of electricity") + CAPITAL_COST = (7, "Plant capital cost") + ASPECT_RATIO = (8, "Plasma aspect ratio") + DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") + TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") + TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") + PULSE_LENGTH = (14, "Pulse length") + PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") + MIN_R0_MAX_TAU_BURN = ( + 16, + "Linear combination of major radius (minimised) and pulse length (maximised)", + ) + NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") + NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") + MAX_Q_MAX_T_PLANT_PULSE_BURN = ( + 19, + "Linear combination of big Q and pulse length (maximised)", + ) + + def __new__(cls, value: int, description: str): + """Create a new FiguresOfMerit enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this figure of merit. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this figure of merit.""" + return self._description_ + + +class PROCESSRunMode(IntEnum): + """Enumeration of the available PROCESS run modes, which determine the behaviour + of the code in various places. This is controlled by the `ioptimz` variable + """ + + EVALUATION = (-2, "Evaluation mode (no optimisation)") + """In this mode, the code will not perform any optimisation, and will instead + simply evaluate the constraints for the given input parameters, which is useful + for testing and for evaluating the performance of a given design point without + trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd + root-finding method from `scipy.optimize`) to seek a *consistent* solution by + varying a subset of the iteration variables until the consistency constraints + (equality constraints whose residuals must be driven to zero) are simultaneously + satisfied; no figure-of-merit is optimised, and the solver simply tries to find + a root of the constraint-residual vector. + """ + OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") + """In this mode, the code will perform optimisation using the VMCON solver + (or a custom solver if specified) to try to find a design point that optimises + the figure of merit while satisfying the constraints. This is the default mode + of operation for PROCESS. + """ + + def __new__(cls, value: int, description: str): + """Create a new PROCESSRunMode enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this run mode. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this run mode.""" + return self._description_ + + +class FiguresOfMerit(IntEnum): + """Enumeration of the available figures of merit (FoM) that can be used as + objective functions for optimisation in PROCESS. + """ + + MAJOR_RADIUS = (1, "Plasma major radius (R₀)") + NEUTRON_WALL_LOAD = (3, "Neutron wall load") + P_TF_PLUS_P_PF = (4, "TF & PF coil power") + FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") + COST_OF_ELECTRICITY = (6, "Cost of electricity") + CAPITAL_COST = (7, "Plant capital cost") + ASPECT_RATIO = (8, "Plasma aspect ratio") + DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") + TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") + TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") + PULSE_LENGTH = (14, "Pulse length") + PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") + MIN_R0_MAX_TAU_BURN = ( + 16, + "Linear combination of major radius (minimised) and pulse length (maximised)", + ) + NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") + NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") + MAX_Q_MAX_T_PLANT_PULSE_BURN = ( + 19, + "Linear combination of big Q and pulse length (maximised)", + ) + + def __new__(cls, value: int, description: str): + """Create a new FiguresOfMerit enum member with description. + + Args: + value: The integer value of the enum member. + description: The description for this figure of merit. + + Returns + ------- + The new enum member with attached description. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Return the description for this figure of merit.""" + return self._description_ + + + +ipnvars: int = 183 """total number of variables available for iteration""" -ipeqns: int = 95 +ipeqns: int = 96 """number of constraint equations available""" ipnfoms: int = 19 @@ -768,6 +960,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 e3cee1d7d0..b37bafb0ce 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -918,6 +918,15 @@ class PhysicsData: 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 = 1.0 """tauratio /1.0/ : ratio of He and pellet particle confinement times""" diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 694e3e2b1b..924012d116 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2424,6 +2424,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 4f249138a66b1822e2e4d65f3b1634361b2af86f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 17:31:12 +0000 Subject: [PATCH 20/37] 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 6de2cd2add..a9afe720b8 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1795,10 +1795,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 @@ -1806,7 +1810,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) @@ -1818,10 +1822,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 1f20501742..92e7aca6eb 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -259,7 +259,7 @@ class IterationVariable: "molflow_plasma_fuelling_vv_injected", "physics", 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 b1ad8b846e..0c7585fcae 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -197,6 +197,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, @@ -208,9 +209,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 - ( @@ -221,9 +225,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, @@ -232,11 +238,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)) @@ -246,9 +256,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, ): @@ -256,9 +266,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): @@ -271,6 +283,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 @@ -331,11 +346,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 ), @@ -386,21 +405,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 ), @@ -410,14 +429,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", @@ -425,11 +444,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", @@ -440,7 +459,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 0bfb8449bd3e31153614a4d24059b781cdb8e948 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 21:01:12 +0000 Subject: [PATCH 21/37] Rename burn-up fraction variable to f_plasma_fuel_burnup across multiple files and update related calculations and tests --- process/core/io/plot/summary.py | 2 +- process/data_structure/physics_variables.py | 2 +- process/models/physics/physics.py | 16 ++++++++-------- process/models/stellarator/stellarator.py | 2 +- tests/unit/models/physics/test_physics.py | 4 ++-- 5 files changed, 13 insertions(+), 13 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index d0c8c8857d..5ffa1cbf60 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -2954,7 +2954,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 b37bafb0ce..597455cbc4 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -206,7 +206,7 @@ class PhysicsData: e_plasma_magnetic_stored: float = 0.0 """Plasma stored magnetic energy (J)""" - burnup: float = 0.0 + f_plasma_fuel_burnup: float = 0.0 """fractional plasma burnup""" burnup_in: float = 0.0 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 924012d116..99e0fd15a9 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -919,7 +919,7 @@ def run(self): # Calculate auxiliary physics related information sbar = 1.0e0 ( - self.data.physics.burnup, + self.data.physics.f_plasma_fuel_burnup, self.data.physics.figmer, self.data.physics.fusrat, self.data.physics.molflow_plasma_fuelling_required, @@ -1468,7 +1468,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 @@ -1500,14 +1500,14 @@ 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 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) / tauratio @@ -1519,12 +1519,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, @@ -2469,8 +2469,8 @@ def outplas(self): po.ovarrf( self.outfile, "Burn-up fraction", - "(burnup)", - self.data.physics.burnup, + "(f_plasma_fuel_burnup)", + self.data.physics.f_plasma_fuel_burnup, "OP ", ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index e9d65132f9..89a7bcb927 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2337,7 +2337,7 @@ def st_phys(self, output): sbar = 1.0e0 ( - self.data.physics.burnup, + self.data.physics.f_plasma_fuel_burnup, self.data.physics.figmer, _fusrat, self.data.physics.molflow_plasma_fuelling_required, diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index b41e11b9b8..da7f1385b3 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -2253,7 +2253,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): monkeypatch.setattr(physics.data.physics, "burnup_in", phyauxparam.burnup_in) ( - burnup, + f_plasma_fuel_burnup, figmer, fusrat, molflow_plasma_fuelling_required, @@ -2274,7 +2274,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): tauratio=phyauxparam.tauratio, ) - 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 b20c1f116471928e60546590200cbdce5fbf9f0c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sun, 22 Mar 2026 21:22:38 +0000 Subject: [PATCH 22/37] Refactor fuel burn-up rate variable to fusrat_total across multiple files and update related calculations and tests --- process/core/io/plot/summary.py | 4 ++-- process/data_structure/physics_variables.py | 6 ++---- process/models/costs/costs.py | 2 +- process/models/physics/physics.py | 18 +++++++++--------- process/models/stellarator/stellarator.py | 8 ++++---- tests/unit/models/test_costs_1990.py | 10 +++++----- 6 files changed, 23 insertions(+), 25 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 5ffa1cbf60..2b664e04ff 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -2952,8 +2952,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 597455cbc4..b5bfb4503b 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -972,11 +972,9 @@ class PhysicsData: f_nd_beam_electron: float = 0.005 """hot beam density / n_e (`iteration variable 7`)""" - f_nd_plasma_carbon_electron: float = 0.0 - """n_carbon / n_e""" - rndfuel: float = 0.0 - """fuel burnup rate (reactions/second)""" +f_nd_plasma_carbon_electron: float = None +"""n_carbon / n_e""" f_nd_plasma_iron_argon_electron: float = 0.0 """n_highZ / n_e""" diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py index 578f54d42f..f9e3599ecb 100644 --- a/process/models/costs/costs.py +++ b/process/models/costs/costs.py @@ -2338,7 +2338,7 @@ def acc2272(self): # New calculation: 2 nuclei * reactions/sec * kg/nucleus * g/kg * sec/day self.data.physics.wtgpd = ( 2.0e0 - * self.data.physics.rndfuel + * self.data.physics.fusrat_total * self.data.physics.m_fuel_amu * constants.UMASS * 1000.0e0 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 99e0fd15a9..bc494b5371 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -919,13 +919,13 @@ def run(self): # Calculate auxiliary physics related information sbar = 1.0e0 ( - self.data.physics.f_plasma_fuel_burnup, - self.data.physics.figmer, - self.data.physics.fusrat, - self.data.physics.molflow_plasma_fuelling_required, - self.data.physics.rndfuel, - self.data.physics.t_alpha_confinement, - self.data.physics.f_alpha_energy_confinement, + 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( self.data.physics.aspect, self.data.physics.nd_plasma_fuel_ions_vol_avg, @@ -2462,8 +2462,8 @@ def outplas(self): po.ovarre( self.outfile, "Fuel burn-up rate (reactions/s)", - "(rndfuel)", - self.data.physics.rndfuel, + "(fusrat_total)", + self.data.physics.fusrat_total, "OP ", ) po.ovarrf( diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 89a7bcb927..a68c005889 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2340,10 +2340,10 @@ def st_phys(self, output): self.data.physics.f_plasma_fuel_burnup, self.data.physics.figmer, _fusrat, - self.data.physics.molflow_plasma_fuelling_required, - self.data.physics.rndfuel, - self.data.physics.t_alpha_confinement, - self.data.physics.f_alpha_energy_confinement, + physics_variables.molflow_plasma_fuelling_required, + _, + physics_variables.t_alpha_confinement, + physics_variables.f_alpha_energy_confinement, ) = self.physics.phyaux( self.data.physics.aspect, self.data.physics.nd_plasma_fuel_ions_vol_avg, diff --git a/tests/unit/models/test_costs_1990.py b/tests/unit/models/test_costs_1990.py index 1975b01849..5c0c99ad0c 100644 --- a/tests/unit/models/test_costs_1990.py +++ b/tests/unit/models/test_costs_1990.py @@ -144,7 +144,7 @@ def test_acc2272(monkeypatch, costs): :param monkeypatch: Mock fixture :type monkeypatch: object """ - monkeypatch.setattr(costs.data.physics, "rndfuel", 7.158e20) + monkeypatch.setattr(costs.data.physics, "fusrat_total", 7.158e20) monkeypatch.setattr(costs.data.physics, "m_fuel_amu", 2.5) monkeypatch.setattr(costs.data.costs, "fkind", 1) monkeypatch.setattr(costs.data.costs, "c2271", 0) @@ -4384,7 +4384,7 @@ class Acc2272Param(NamedTuple): wtgpd: Any = None - rndfuel: Any = None + fusrat_total: Any = None m_fuel_amu: Any = None @@ -4410,7 +4410,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, @@ -4426,7 +4426,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, @@ -4463,7 +4463,7 @@ def test_acc2272_rut(acc2272param, monkeypatch, costs): monkeypatch.setattr(costs.data.physics, "wtgpd", acc2272param.wtgpd) - monkeypatch.setattr(costs.data.physics, "rndfuel", acc2272param.rndfuel) + monkeypatch.setattr(costs.data.physics, "fusrat_total", acc2272param.fusrat_total) monkeypatch.setattr(costs.data.physics, "m_fuel_amu", acc2272param.m_fuel_amu) From d2550bcd958685a4630a276bd7c031fa330cc35d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 11:21:44 +0000 Subject: [PATCH 23/37] Refactor plasma exhaust calculations to plasma fuelling, updating related plotting functions for tritium, deuterium, and alpha particle flow rates --- process/core/io/plot/summary.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/summary.py b/process/core/io/plot/summary.py index 2b664e04ff..bcbd310cf8 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -58,7 +58,7 @@ ElectronCyclotron, ) from process.models.physics.density_limit import DensityLimitModel -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.physics import ( diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py index 0c7585fcae..2b3a30e321 100644 --- a/process/models/physics/exhaust.py +++ b/process/models/physics/exhaust.py @@ -2,10 +2,6 @@ import logging -import matplotlib.pyplot as plt -import numpy as np - -import process.core.io.mfile as mf from process.core import constants from process.core import process_output as po from process.core.model import Model @@ -194,274 +190,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 d1c65bf31b47d3cde822ce77a69c1e08e00786c7 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 11:59:43 +0000 Subject: [PATCH 24/37] 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 65b8ba3ad89475177ec947a89e508649092523d6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 12:24:48 +0000 Subject: [PATCH 25/37] 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 | 87 +++++------------------ process/models/stellarator/stellarator.py | 25 +++---- 4 files changed, 30 insertions(+), 89 deletions(-) diff --git a/process/main.py b/process/main.py index 3d99d72545..788f2ece71 100644 --- a/process/main.py +++ b/process/main.py @@ -86,6 +86,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, ) @@ -643,6 +644,7 @@ def __init__(self, data: DataStructure): self.plasma_density_limit = PlasmaDensityLimit() self.plasma_exhaust = PlasmaExhaust() self.sauter_bootstrap_current = SauterBootstrapCurrent() + self.plasma_fuelling = PlasmaFuelling() self.plasma_bootstrap_current = PlasmaBootstrapCurrent( plasma_profile=self.plasma_profile, sauter_bootstrap=self.sauter_bootstrap_current, @@ -666,6 +668,7 @@ def __init__(self, data: DataStructure): plasma_fields=self.plasma_fields, plasma_dia_current=self.plasma_dia_current, plasma_geometry=self.plasma_geom, + 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 bc494b5371..c05a381258 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -31,7 +31,8 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust - from process.models.physics.l_h_transition import PlasmaConfinementTransition + from process.models.physics.fuelling import PlasmaFuelling +from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.plasma_current import ( PlasmaCurrent, PlasmaDiamagneticCurrent, @@ -189,6 +190,7 @@ def __init__( plasma_fields: PlasmaFields, plasma_dia_current: PlasmaDiamagneticCurrent, plasma_geometry: PlasmaGeom, + plasma_fuelling: PlasmaFuelling, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -205,6 +207,7 @@ def __init__( self.fields = plasma_fields self.dia_current = plasma_dia_current self.geometry = plasma_geometry + self.fuelling = plasma_fuelling def output(self) -> None: """Output plasma physics information.""" @@ -912,35 +915,25 @@ def run(self): self.data.physics.ind_plasma_internal_norm, ) - self.data.physics.e_plasma_magnetic_stored = ( - 0.5e0 * self.data.physics.ind_plasma * self.data.physics.plasma_current**2 + physics_variables.e_plasma_magnetic_stored = ( + 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( - self.data.physics.aspect, - self.data.physics.nd_plasma_fuel_ions_vol_avg, - self.data.physics.fusden_total, - self.data.physics.fusden_alpha_total, - self.data.physics.plasma_current, - sbar, - self.data.physics.nd_plasma_alphas_vol_avg, - self.data.physics.t_energy_confinement, - self.data.physics.vol_plasma, - self.data.physics.burnup_in, - self.data.physics.tauratio, + physics_variables.fusden_alpha_total, + physics_variables.nd_plasma_alphas_vol_avg, + 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, ) - self.data.physics.ntau, self.data.physics.nTtau = ( + physics_variables.ntau, physics_variables.nTtau = ( self.confinement.calculate_double_and_triple_product( nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, t_energy_confinement=self.data.physics.t_energy_confinement, @@ -1425,18 +1418,12 @@ def plasma_composition(self): @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, burnup_in: float, tauratio: float, - ) -> tuple[float, float, float, float, float, float, float, float]: + ) -> tuple[float, float]: """Auxiliary physics quantities Parameters @@ -1468,22 +1455,12 @@ 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 @@ -1494,41 +1471,9 @@ def phyaux( else nd_plasma_alphas_vol_avg / fusden_alpha_total ) - # 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 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) - / tauratio - ) - else: - burnup = 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 a68c005889..7b31377c55 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2335,27 +2335,18 @@ def st_phys(self, output): # Calculate auxiliary physics related information # for the rest of the code - sbar = 1.0e0 ( - self.data.physics.f_plasma_fuel_burnup, - self.data.physics.figmer, - _fusrat, - physics_variables.molflow_plasma_fuelling_required, - _, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.physics.phyaux( - self.data.physics.aspect, - self.data.physics.nd_plasma_fuel_ions_vol_avg, - self.data.physics.fusden_total, - self.data.physics.fusden_alpha_total, - self.data.physics.plasma_current, - sbar, - self.data.physics.nd_plasma_alphas_vol_avg, - self.data.physics.t_energy_confinement, - self.data.physics.vol_plasma, - self.data.physics.burnup_in, - self.data.physics.tauratio, + physics_variables.fusden_alpha_total, + physics_variables.nd_plasma_alphas_vol_avg, + physics_variables.t_energy_confinement, + ) + + 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 cc3e2de320dc0e22a2911a529087d40b4d92d3be Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 13:40:30 +0000 Subject: [PATCH 26/37] 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 ff205b40ea167b3844e63601ffc556727a35b777 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 14:56:41 +0000 Subject: [PATCH 27/37] 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 | 612 +---------------------------- 2 files changed, 102 insertions(+), 611 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 c05a381258..859c2ba21b 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2348,617 +2348,7 @@ def outplas(self): self.plasma_bootstrap_current.output() self.dia_current.output() - po.oheadr(self.outfile, "Plasma Fuelling") - po.ovarre( - self.outfile, - "Ratio of He and pellet particle confinement times", - "(tauratio)", - self.data.physics.tauratio, - ) - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - self.data.physics.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)", - self.data.physics.fusrat_total, - "OP ", - ) - po.ovarrf( - self.outfile, - "Burn-up fraction", - "(f_plasma_fuel_burnup)", - self.data.physics.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)", - self.data.reinke.impvardiv, - ) - po.ovarre( - self.outfile, - "Minimum Impurity fraction from Reinke", - "(fzmin)", - self.data.reinke.fzmin, - "OP ", - ) - po.ovarre( - self.outfile, - "Actual Impurity fraction", - "(fzactual)", - self.data.reinke.fzactual, - ) - - def output_temperature_density_profile_info(self) -> None: - """Output information about plasma temperature and density profiles.""" - po.oheadr(self.outfile, "Plasma Density and Temperature Profiles") - po.ovarrf( - self.outfile, - "Number of radial points in plasma profiles", - "(n_plasma_profile_elements)", - self.data.physics.n_plasma_profile_elements, - ) - po.ovarin( - self.outfile, - "Plasma profile model selected", - "(i_plasma_pedestal)", - self.data.physics.i_plasma_pedestal, - ) - po.ocmmnt( - self.outfile, - f"Plasma profile model selected: " - f"{PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal).description}", - ) - - if ( - PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) - == PlasmaProfileShapeType.PEDESTAL_PROFILE - ) and ( - self.data.physics.nd_plasma_electron_on_axis - < self.data.physics.nd_plasma_pedestal_electron - ): - logger.error("Central density is less than pedestal density") - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "----------------------------") - po.osubhd(self.outfile, "Temperature:") - po.ovarrf( - self.outfile, - "Temperature profile index (αₜ)", - "(alphat)", - self.data.physics.alphat, - ) - po.ovarrf( - self.outfile, - "Temperature profile index beta (βₜ)", - "(tbeta)", - self.data.physics.tbeta, - ) - po.oblnkl(self.outfile) - - if ( - PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) - == PlasmaProfileShapeType.PEDESTAL_PROFILE - ): - po.ovarrf( - self.outfile, - "Temperature pedestal r/a location (ρₜ,pedestal)", - "(radius_plasma_pedestal_temp_norm)", - self.data.physics.radius_plasma_pedestal_temp_norm, - ) - - po.ovarrf( - self.outfile, - "Electron temperature pedestal height (Tₑ,pedestal) (keV)", - "(temp_plasma_pedestal_kev)", - self.data.physics.temp_plasma_pedestal_kev, - ) - if 78 in numerics.icc: - po.ovarrf( - self.outfile, - "Electron temperature at separatrix (Tₑ,ₛₑₚ) (keV)", - "(temp_plasma_separatrix_kev)", - self.data.physics.temp_plasma_separatrix_kev, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Electron temperature at separatrix (Tₑ,ₛₑₚ) (keV)", - "(temp_plasma_separatrix_kev)", - self.data.physics.temp_plasma_separatrix_kev, - ) - po.oblnkl(self.outfile) - - po.ovarrf( - self.outfile, - "Volume averaged electron temperature (⟨Tₑ⟩) (keV)", - "(temp_plasma_electron_vol_avg_kev)", - self.data.physics.temp_plasma_electron_vol_avg_kev, - ) - po.ovarrf( - self.outfile, - "Electron temperature on axis (Tₑ₀) (keV)", - "(temp_plasma_electron_on_axis_kev)", - self.data.physics.temp_plasma_electron_on_axis_kev, - "OP ", - ) - po.ovarrf( - self.outfile, - "Volume averaged density weighted electron temperature (⟨Tₑ⟩ₙ) (keV)", - "(temp_plasma_electron_density_weighted_kev)", - self.data.physics.temp_plasma_electron_density_weighted_kev, - "OP ", - ) - po.ovarrf( - self.outfile, - "Ratio of electron density weighted to volume averaged temperature", - "(f_temp_plasma_electron_density_vol_avg)", - self.data.physics.f_temp_plasma_electron_density_vol_avg, - "OP ", - ) - po.oblnkl(self.outfile) - po.ovarrf( - self.outfile, - "Ratio of ion to electron volume-averaged temperature", - "(f_temp_plasma_ion_electron)", - self.data.physics.f_temp_plasma_ion_electron, - "IP ", - ) - po.oblnkl(self.outfile) - po.ovarrf( - self.outfile, - "Volume averaged ion temperature (⟨Tᵢ⟩) (keV)", - "(temp_plasma_ion_vol_avg_kev)", - self.data.physics.temp_plasma_ion_vol_avg_kev, - ) - po.ovarrf( - self.outfile, - "Ion temperature on axis (Tᵢ₀) (keV)", - "(temp_plasma_ion_on_axis_kev)", - self.data.physics.temp_plasma_ion_on_axis_kev, - "OP ", - ) - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "----------------------------") - - po.osubhd(self.outfile, "Density:") - po.ovarrf( - self.outfile, - "Density profile factor (αₙ)", - "(alphan)", - self.data.physics.alphan, - ) - po.oblnkl(self.outfile) - if ( - PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) - == PlasmaProfileShapeType.PEDESTAL_PROFILE - ): - po.ovarrf( - self.outfile, - "Density pedestal r/a location (ρₙ,pedestal)", - "(radius_plasma_pedestal_density_norm)", - self.data.physics.radius_plasma_pedestal_density_norm, - ) - if self.data.physics.f_nd_plasma_pedestal_greenwald >= 0e0: - po.ovarre( - self.outfile, - "Electron density pedestal height (nₑ_pedestal) (/m³)", - "(nd_plasma_pedestal_electron)", - self.data.physics.nd_plasma_pedestal_electron, - "OP ", - ) - else: - po.ovarre( - self.outfile, - "Electron density pedestal height (nₑ_pedestal) (/m³)", - "(nd_plasma_pedestal_electron)", - self.data.physics.nd_plasma_pedestal_electron, - ) - # must be assigned to their exisiting values# - fgwped_out = ( - self.data.physics.nd_plasma_pedestal_electron - / self.data.physics.nd_plasma_electron_max_array[6] - ) - fgwsep_out = ( - self.data.physics.nd_plasma_separatrix_electron - / self.data.physics.nd_plasma_electron_max_array[6] - ) - if self.data.physics.f_nd_plasma_pedestal_greenwald >= 0e0: - self.data.physics.f_nd_plasma_pedestal_greenwald = ( - self.data.physics.nd_plasma_pedestal_electron - / self.data.physics.nd_plasma_electron_max_array[6] - ) - if self.data.physics.f_nd_plasma_separatrix_greenwald >= 0e0: - self.data.physics.f_nd_plasma_separatrix_greenwald = ( - self.data.physics.nd_plasma_separatrix_electron - / self.data.physics.nd_plasma_electron_max_array[6] - ) - - po.ovarre( - self.outfile, - "Pedestal Greenwald fraction", - "(fgwped_out)", - fgwped_out, - ) - po.ovarre( - self.outfile, - "Electron density at separatrix (nₑ,ₛₑₚ) (/m³)", - "(nd_plasma_separatrix_electron)", - self.data.physics.nd_plasma_separatrix_electron, - ) - po.ovarre( - self.outfile, - "Separatrix Greenwald fraction", - "(fgwsep_out)", - fgwsep_out, - ) - po.oblnkl(self.outfile) - - po.ovarre( - self.outfile, - "Volume averaged electron number density (⟨nₑ⟩) (/m³)", - "(nd_plasma_electrons_vol_avg)", - self.data.physics.nd_plasma_electrons_vol_avg, - ) - po.ovarre( - self.outfile, - "Electron number density on axis (nₑ₀) (/m³)", - "(nd_plasma_electron_on_axis)", - self.data.physics.nd_plasma_electron_on_axis, - "OP ", - ) - po.ovarre( - self.outfile, - "Line-averaged electron number density (ñₑ) (/m³)", - "(nd_plasma_electron_line)", - self.data.physics.nd_plasma_electron_line, - "OP ", - ) - if self.data.stellarator.istell == 0: - po.ovarre( - self.outfile, - "Greenwald fraction (f_GW)", - "(dnla_gw)", - self.data.physics.nd_plasma_electron_line - / self.data.physics.nd_plasma_electron_max_array[6], - "OP ", - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Total ion volume averaged number density (⟨nᵢ⟩) (/m³)", - "(nd_plasma_ions_total_vol_avg)", - self.data.physics.nd_plasma_ions_total_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuel ion volume averaged number density (⟨n_fuel⟩) (/m³)", - "(nd_plasma_fuel_ions_vol_avg)", - self.data.physics.nd_plasma_fuel_ions_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Total impurity volume averaged number density with Z > 2 (⟨nᵢₘₚ⟩) (/m³)", - "(nd_plasma_impurities_vol_avg)", - self.data.physics.nd_plasma_impurities_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermalised alpha volume averaged number density (⟨n_alpha⟩) (/m³)", - "(nd_plasma_alphas_vol_avg)", - self.data.physics.nd_plasma_alphas_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermalised alpha to electron number density ratio (⟨n_alpha⟩/⟨nₑ⟩)", - "(f_nd_alpha_electron)", - self.data.physics.f_nd_alpha_electron, - ) - po.ovarre( - self.outfile, - "Proton volume averaged number density (⟨nₚ⟩) (/m³)", - "(nd_plasma_protons_vol_avg)", - self.data.physics.nd_plasma_protons_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Proton to electron volume averaged number density ratio (⟨nₚ⟩/⟨nₑ⟩)", - "(f_nd_protium_electrons)", - self.data.physics.f_nd_protium_electrons, - "OP ", - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Hot beam ion volume averaged number density (⟨n_beam⟩) (/m³)", - "(nd_beam_ions)", - self.data.physics.nd_beam_ions, - "OP ", - ) - po.ovarre( - self.outfile, - "Hot beam ion to electron number density ratio (⟨n_beam⟩/⟨nₑ⟩)", - "(f_nd_beam_electron)", - self.data.physics.f_nd_beam_electron, - "OP ", - ) - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "----------------------------") - - po.osubhd(self.outfile, "Pressure:") - - po.ovarrf( - self.outfile, - "Pressure profile index (αₚ)", - "(alphap)", - self.data.physics.alphap, - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Plasma thermal pressure on axis (p₀) (Pa)", - "(pres_plasma_thermal_on_axis)", - self.data.physics.pres_plasma_thermal_on_axis, - "OP ", - ) - po.ovarre( - self.outfile, - "Volume averaged plasma thermal pressure (⟨p⟩) (Pa)", - "(pres_plasma_thermal_vol_avg)", - self.data.physics.pres_plasma_thermal_vol_avg, - "OP ", - ) - # As array output is not currently supported, each element is output as a float - # instance - # Output plasma pressure profiles to mfile - for i in range(len(self.data.physics.pres_plasma_thermal_total_profile)): - po.ovarre( - self.mfile, - f"Total plasma pressure at point {i}", - f"(pres_plasma_thermal_total_profile{i})", - self.data.physics.pres_plasma_thermal_total_profile[i], - ) - for i in range(len(self.data.physics.pres_plasma_electron_profile)): - po.ovarre( - self.mfile, - f"Total plasma electron pressure at point {i}", - f"(pres_plasma_electron_profile{i})", - self.data.physics.pres_plasma_electron_profile[i], - ) - for i in range(len(self.data.physics.pres_plasma_ion_total_profile)): - po.ovarre( - self.mfile, - f"Total plasma ion pressure at point {i}", - f"(pres_plasma_ion_total_profile{i})", - self.data.physics.pres_plasma_ion_total_profile[i], - ) - for i in range(len(self.data.physics.pres_plasma_fuel_profile)): - po.ovarre( - self.mfile, - f"Total plasma fuel pressure at point {i}", - f"(pres_plasma_fuel_profile{i})", - self.data.physics.pres_plasma_fuel_profile[i], - ) - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "----------------------------") - po.oblnkl(self.outfile) - - po.ocmmnt(self.outfile, "Impurities:") - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "Plasma ion densities / electron density:") - po.oblnkl(self.outfile) - - for imp in range(N_IMPURITIES): - # MDK Update f_nd_impurity_electrons, as this will make the ITV output - # work correctly. - self.data.impurity_radiation.f_nd_impurity_electrons[imp] = ( - self.data.impurity_radiation.f_nd_impurity_electron_array[imp] - ) - str1 = ( - self.data.impurity_radiation - .impurity_arr_label[imp] - .item() - .replace("_", "") - + " concentration" - ) - str2 = f"(f_nd_impurity_electrons({imp + 1:02}))" - # MDK Add output flag for H which is calculated. - if imp == 0: - po.ovarre( - self.outfile, - str1, - str2, - self.data.impurity_radiation.f_nd_impurity_electrons[imp], - "OP ", - ) - else: - po.ovarre( - self.outfile, - str1, - str2, - self.data.impurity_radiation.f_nd_impurity_electrons[imp], - ) - if self.data.impurity_radiation.f_nd_impurity_electrons[imp] != 0.0: # noqa: RUF069 - for i in range(self.data.physics.n_plasma_profile_elements): - po.ovarre( - self.mfile, - str1 + f" at point {i}", - f"(f_nd_impurity_electrons{imp}_{i})", - self.data.impurity_radiation.f_nd_impurity_electrons[imp] - * self.plasma_profile.neprofile.profile_y[i], - "OP ", - ) - - po.oblnkl(self.outfile) - po.ovarrf( - self.outfile, - "Volume averaged plasma effective charge (⟨Zₑ⟩)", - "(n_charge_plasma_effective_vol_avg)", - self.data.physics.n_charge_plasma_effective_vol_avg, - "OP ", - ) - for i in range(len(self.data.physics.n_charge_plasma_effective_profile)): - po.ovarre( - self.mfile, - "Effective charge at point", - f"(n_charge_plasma_effective_profile{i})", - self.data.physics.n_charge_plasma_effective_profile[i], - "OP ", - ) - - for imp in range(N_IMPURITIES): - for i in range(self.data.physics.n_plasma_profile_elements): - po.ovarre( - self.mfile, - "Impurity charge at point", - f"(n_charge_plasma_profile{imp}_{i})", - self.data.impurity_radiation.n_charge_impurity_profile[imp][i], - "OP ", - ) - - po.ovarrf( - self.outfile, - "Volume averaged mass-weighted plasma effective charge (⟨Zₑ⟩ₘ)", - "(n_charge_plasma_effective_mass_weighted_vol_avg)", - self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg, - "OP ", - ) - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "----------------------------") - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "Masses:") - po.oblnkl(self.outfile) - - po.ovarre( - self.outfile, - "Average mass of all ions (amu)", - "(m_ions_total_amu)", - self.data.physics.m_ions_total_amu, - "OP ", - ) - po.ovarre( - self.outfile, - "Total mass of all ions in plasma (kg)", - "(m_plasma_ions_total)", - self.data.physics.m_plasma_ions_total, - "OP ", - ) - po.ovarre( - self.outfile, - "Average mass of all fuel ions (amu)", - "(m_fuel_amu)", - self.data.physics.m_fuel_amu, - "OP ", - ) - po.ovarre( - self.outfile, - "Total mass of all fuel ions in plasma (kg)", - "(m_plasma_fuel_ions)", - self.data.physics.m_plasma_fuel_ions, - "OP ", - ) - po.ovarre( - self.outfile, - "Average mass of all beam ions (amu)", - "(m_beam_amu)", - self.data.physics.m_beam_amu, - "OP ", - ) - po.ovarre( - self.outfile, - "Total mass of all alpha particles in plasma (kg)", - "(m_plasma_alpha)", - self.data.physics.m_plasma_alpha, - "OP ", - ) - po.ovarre( - self.outfile, - "Total mass of all electrons in plasma (kg)", - "(m_plasma_electron)", - self.data.physics.m_plasma_electron, - "OP ", - ) - po.ovarre( - self.outfile, - "Total mass of the plasma (kg)", - "(m_plasma)", - self.data.physics.m_plasma, - "OP ", - ) - - # Issue 558 - addition of constraint 76 to limit the value of - # nd_plasma_separatrix_electron, in proportion with the ballooning parameter - # and Greenwald density - if 76 in numerics.icc: - po.ovarre( - self.outfile, - "Critical ballooning parameter value", - "(alpha_crit)", - self.data.physics.alpha_crit, - ) - po.ovarre( - self.outfile, - "Critical electron density at separatrix (/m3)", - "(nd_plasma_separatrix_electron_eich_max)", - self.data.physics.nd_plasma_separatrix_electron_eich_max, - ) + self.fuelling.output_fuelling_info() @staticmethod @nb.njit(cache=True) From 6e67b2b0880c6e13f9e087b47bde7dd9e43af668 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 14:59:02 +0000 Subject: [PATCH 28/37] :fire: Remove tauratio variable from input variables and related tests for cleanup --- process/core/input.py | 1 - process/data_structure/physics_variables.py | 2 -- process/models/physics/fuelling.py | 6 ------ tests/unit/models/physics/test_physics.py | 4 ---- 4 files changed, 13 deletions(-) diff --git a/process/core/input.py b/process/core/input.py index 989f5ca256..d6653ecc7a 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -870,7 +870,6 @@ def __post_init__(self): "t_plasma_energy_confinement_max": InputVariable( "physics", float, range=(0.1, 100.0) ), - "tauratio": InputVariable("physics", float, range=(0.1, 100.0)), "n_beam_decay_lengths_core_required": InputVariable( "current_drive", float, range=(0.0, 10.0) ), diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index b5bfb4503b..cdf070b2de 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -927,8 +927,6 @@ class PhysicsData: f_molflow_plasma_fuelling_helium3: float = None """fraction of plasma fuelling that is helium-3""" - tauratio: float = 1.0 - """tauratio /1.0/ : ratio of He and pellet particle confinement times""" q95_min: float = 0.0 """lower limit for edge safety factor""" 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/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index da7f1385b3..ec9004670f 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -2148,7 +2148,6 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - tauratio: Any = None burnup_in: Any = None @@ -2195,7 +2194,6 @@ class PhyauxParam(NamedTuple): "phyauxparam", [ PhyauxParam( - tauratio=1, burnup_in=0, aspect=3, nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19, @@ -2214,7 +2212,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, @@ -2248,7 +2245,6 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr(physics.data.physics, "tauratio", phyauxparam.tauratio) monkeypatch.setattr(physics.data.physics, "burnup_in", phyauxparam.burnup_in) From 98c39b690b5a3a4e1fe0a3181f479f1f4928a238 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 15:11:36 +0000 Subject: [PATCH 29/37] Add tritium and deuterium burnup calculations and update physics variables --- process/data_structure/physics_variables.py | 9 +- process/models/physics/fuelling.py | 97 ++++++++++++++++++++- process/models/physics/physics.py | 13 +++ tests/unit/models/physics/test_physics.py | 2 +- 4 files changed, 116 insertions(+), 5 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index cdf070b2de..5988f83cfb 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -207,7 +207,14 @@ class PhysicsData: """Plasma stored magnetic energy (J)""" f_plasma_fuel_burnup: float = 0.0 - """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 = 0.0 """fractional plasma burnup user input""" 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 859c2ba21b..5d2a45f0c8 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -932,6 +932,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/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index ec9004670f..c76efb914f 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -2148,7 +2148,6 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - burnup_in: Any = None aspect: Any = None @@ -2246,6 +2245,7 @@ def test_phyaux(phyauxparam, monkeypatch, physics): """ + monkeypatch.setattr(physics.data.physics, "burnup_in", phyauxparam.burnup_in) ( From f959ac5fea0f2e12790d9e7ca7c1383e36a79d18 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 15:37:35 +0000 Subject: [PATCH 30/37] :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 | 4 +--- process/models/physics/fuelling.py | 7 ------- process/models/vacuum.py | 8 ++++---- tests/unit/models/physics/test_physics.py | 4 ++-- tests/unit/models/test_vacuum.py | 2 +- 5 files changed, 8 insertions(+), 17 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 5988f83cfb..a443de1b5d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -914,8 +914,6 @@ class PhysicsData: in which case q95 = mean edge safety factor qbar) """ - molflow_plasma_fuelling_required: float = 0.0 - """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)""" @@ -924,7 +922,7 @@ class PhysicsData: """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""" 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 f4815cfce3..de858858e6 100644 --- a/process/models/vacuum.py +++ b/process/models/vacuum.py @@ -49,7 +49,7 @@ def run(self, output: bool = False): # MDK Check this!! gasld = ( 2.0e0 - * self.data.physics.molflow_plasma_fuelling_required + * self.data.physics.molflow_plasma_fuelling_vv_injected * self.data.physics.m_fuel_amu * constants.UMASS ) @@ -117,7 +117,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 = ( - self.data.physics.molflow_plasma_fuelling_required + self.data.physics.molflow_plasma_fuelling_vv_injected / self.data.vacuum.molflow_vac_pumps ) @@ -161,8 +161,8 @@ def vacuum_simple(self, output) -> float: process_output.ovarre( self.outfile, "Plasma fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - self.data.physics.molflow_plasma_fuelling_required, + "(molflow_plasma_fuelling_vv_injected)", + self.data.physics.molflow_plasma_fuelling_vv_injected, "OP ", ) process_output.ocmmnt( diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index c76efb914f..01a782ba21 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -2252,7 +2252,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, _, @@ -2276,7 +2276,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/models/test_vacuum.py b/tests/unit/models/test_vacuum.py index 30b27cde0a..d538c4d77b 100644 --- a/tests/unit/models/test_vacuum.py +++ b/tests/unit/models/test_vacuum.py @@ -39,7 +39,7 @@ def test_simple_model(self, monkeypatch, vacuum): """ monkeypatch.setattr( vacuum.data.physics, - "molflow_plasma_fuelling_required", + "molflow_plasma_fuelling_vv_injected", 7.5745668997694112e22, ) monkeypatch.setattr(vacuum.data.physics, "a_plasma_surface", 1500.3146527709359) From df353649cb4f994177fe57f1f7c96bcdfc78a858 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 23 Mar 2026 16:01:39 +0000 Subject: [PATCH 31/37] 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 512dfa2c43e32f6a850e94e762eb6125d9507388 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 24 Mar 2026 08:58:16 +0000 Subject: [PATCH 32/37] Add plasma fuelling loss calculations and update related variables in the PlasmaFuelling class --- process/data_structure/physics_variables.py | 9 ++++ process/models/physics/fuelling.py | 59 +++++++++++++++++++++ process/models/physics/physics.py | 18 +------ 3 files changed, 69 insertions(+), 17 deletions(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index a443de1b5d..1caa85f1c4 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -924,6 +924,15 @@ class PhysicsData: 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""" 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 5d2a45f0c8..288c523c47 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -928,23 +928,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 9517a75e158bacf3ef76d464477d6ed3fcff9a05 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 24 Mar 2026 09:38:10 +0000 Subject: [PATCH 33/37] 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 51877eacf4ed456a496d307f739f08e37330ba2a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 25 Mar 2026 13:29:06 +0000 Subject: [PATCH 34/37] Refactor plotting indices in main_plot function to accommodate additional figures --- process/core/io/plot/summary.py | 35 ++++++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 9 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index bcbd310cf8..bcac51e876 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -15073,16 +15073,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 == PlasmaShapeModelType.SAUTER: 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 " f"(i_plasma_shape == {PlasmaShapeModelType.SAUTER}). " @@ -15101,10 +15118,10 @@ 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) plot_plasma_current_comparison(figs[14].add_subplot(224), m_file, scan) @@ -15325,9 +15342,9 @@ def main_plot( ax25 = figs[36].add_subplot(221) PlasmaExhaust().plot_tritium_flow_contour(ax25, m_file, scan) - ax26 = figs[31].add_subplot(222) + ax26 = figs[36].add_subplot(222) PlasmaExhaust().plot_deuterium_flow_contour(ax26, m_file, scan) - ax27 = figs[31].add_subplot(223) + ax27 = figs[36].add_subplot(223) PlasmaExhaust().plot_alpha_flow_contour(ax27, m_file, scan) From b58e4c583e9bbc4a3307f05b0c5cc718290a07f4 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 27 Mar 2026 09:07:46 +0000 Subject: [PATCH 35/37] 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 7fa71d39b0b83d666c6cdc70aa53036a49b2b2b8 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 31 Mar 2026 10:39:05 +0100 Subject: [PATCH 36/37] Refactor Physics model and unit tests to integrate PlasmaFuelling class and update parameter descriptions --- process/models/physics/physics.py | 22 +----- tests/unit/models/physics/test_physics.py | 90 +---------------------- 2 files changed, 6 insertions(+), 106 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 288c523c47..3f95e0e755 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1425,28 +1425,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). - burnup_in: float - fractional plasma burnup user input - tauratio: float - ratio of He and pellet particle confinement times + Energy confinement time (s). Returns ------- diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index 01a782ba21..5fe4ed59c3 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -6,6 +6,7 @@ import pytest from process.core import constants +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -2148,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 @@ -2193,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, ), ], @@ -2244,44 +2189,15 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - - - monkeypatch.setattr(physics.data.physics, "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, - burnup_in=phyauxparam.burnup_in, - tauratio=phyauxparam.tauratio, - ) - - 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) From 49236adce7a69062d12ade9ed8b86dc095b5af83 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 2 Jun 2026 10:12:40 +0100 Subject: [PATCH 37/37] Post rebase fixes --- process/core/io/plot/summary.py | 127 ++--- process/core/solver/constraints.py | 56 -- process/core/solver/iteration_variables.py | 8 +- process/data_structure/numerics.py | 287 ---------- process/data_structure/physics_variables.py | 85 ++- process/main.py | 1 + process/models/physics/fuelling.py | 86 ++- process/models/physics/fusion_reactions.py | 12 +- process/models/physics/physics.py | 576 +++++++++++++++++++- process/models/stellarator/stellarator.py | 16 +- tests/unit/models/physics/test_physics.py | 1 - 11 files changed, 697 insertions(+), 558 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index bcac51e876..03c0a1af0c 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -15123,42 +15123,42 @@ def main_plot( plot_ebw_ecrh_coupling_graph(figs[14].add_subplot(111), m_file, scan) - plot_bootstrap_comparison(figs[14].add_subplot(221), m_file, scan) - 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_max_normalised_beta_comparison(figs[16].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) + 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_max_normalised_beta_comparison(figs[17].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_collision_time_profile(figs[17].add_subplot(234), m_file, scan) - plot_collision_frequency_profile(figs[17].add_subplot(235), m_file, scan) - plot_mean_free_path_profile(figs[17].add_subplot(236), 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_collision_time_profile(figs[18].add_subplot(234), m_file, scan) + plot_collision_frequency_profile(figs[18].add_subplot(235), m_file, scan) + plot_mean_free_path_profile(figs[18].add_subplot(236), m_file, scan) - plot_ion_slowing_down_time_profile(figs[18].add_subplot(231), m_file, scan) + plot_ion_slowing_down_time_profile(figs[19].add_subplot(231), m_file, scan) - plot_resistivity_profile(figs[18].add_subplot(232), m_file, scan) + plot_resistivity_profile(figs[19].add_subplot(232), m_file, scan) plot_detailed_plasma_parameters( - figs[18].add_subplot(233), fig=figs[18], mfile=m_file, scan=scan + figs[19].add_subplot(233), fig=figs[19], mfile=m_file, scan=scan ) - ax_electron_freq = figs[19].add_subplot(211) + ax_electron_freq = figs[20].add_subplot(211) plot_electron_frequency_profile(ax_electron_freq, m_file, scan) - ax_ion_freq = figs[19].add_subplot(413, sharex=ax_electron_freq) + ax_ion_freq = figs[20].add_subplot(413, sharex=ax_electron_freq) plot_ion_frequency_profile(ax_ion_freq, m_file, scan) - ax_larmor = figs[19].add_subplot(414, sharex=ax_electron_freq) + ax_larmor = figs[20].add_subplot(414, sharex=ax_electron_freq) plot_larmor_radius_profile(ax_larmor, m_file, scan) - figs[19].subplots_adjust(hspace=0.5) + figs[20].subplots_adjust(hspace=0.5) # Plot poloidal cross-section poloidal_cross_section( - figs[20].add_subplot(121, aspect="equal"), + figs[21].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -15168,7 +15168,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[20].add_subplot(122, aspect="equal"), + figs[21].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -15176,27 +15176,27 @@ def main_plot( ) # Plot color key - ax17 = figs[20].add_subplot(222) + ax17 = figs[21].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[21].add_subplot(111, aspect="equal"), + figs[22].add_subplot(111, aspect="equal"), m_file, scan, radial_build, colour_scheme, ) - ax18 = figs[22].add_subplot(211) + ax18 = figs[23].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[23].add_subplot(211) + ax185 = figs[24].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[23].add_subplot(212) + ax18b = figs[24].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) @@ -15204,36 +15204,36 @@ 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) == TFConductorModel.SUPERCONDUCTING: # TF coil with WP - ax19 = figs[24].add_subplot(221, aspect="equal") + ax19 = figs[25].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[24]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[25]) if ( m_file.get("i_tf_turn_type", scan=scan) == SuperconductingTFTurnType.CROSS_CONDUCTOR ): - ax20 = figs[25].add_subplot(325, aspect="equal") + ax20 = figs[26].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_croco_turn(ax20, figs[25], m_file, scan) + plot_tf_croco_turn(ax20, figs[26], m_file, scan) elif ( m_file.get("i_tf_turn_type", scan=scan) == SuperconductingTFTurnType.CABLE_IN_CONDUIT ): # TF coil turn structure - ax20 = figs[25].add_subplot(325, aspect="equal") + ax20 = figs[26].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[25], m_file, scan) + plot_tf_cable_in_conduit_turn(ax20, figs[26], m_file, scan) if ( m_file.get("i_tf_turn_type", scan=scan) == SuperconductingTFTurnType.CROSS_CONDUCTOR ): - plot_205 = figs[25].add_subplot(223, aspect="equal") + plot_205 = figs[26].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) plot_corc_cable_geometry( plot_205, @@ -15255,8 +15255,8 @@ def main_plot( dx_hts_tape_hastelloy=m_file.get("dx_tf_hts_tape_hastelloy", scan=scan), show_legend=True, ) - plot_tf_corc_cable_summary_box(plot_205, figs[25], m_file, scan) - ax_hts_tape = figs[25].add_subplot(339) + plot_tf_corc_cable_summary_box(plot_205, figs[26], m_file, scan) + ax_hts_tape = figs[26].add_subplot(339) ax_hts_tape.set_position([0.75, 0.1, 0.2, 0.2]) plot_hts_tape_geometry( axis=ax_hts_tape, @@ -15272,9 +15272,9 @@ def main_plot( m_file.get("i_tf_turn_type", scan=scan) == SuperconductingTFTurnType.CABLE_IN_CONDUIT ): - plot_205 = figs[25].add_subplot(223, aspect="equal") + plot_205 = figs[26].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[25], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[26], m_file, scan) plot_quench_time_evolution( tau_discharge=m_file.get("t_tf_superconductor_quench", scan=scan), b_peak=m_file.get("b_tf_inboard_peak_with_ripple", scan=scan), @@ -15292,60 +15292,51 @@ def main_plot( "a_tf_turn_cable_space_no_void", scan=scan ), a_tf_turn=m_file.get("a_tf_turn", scan=scan), - axes_1=figs[26].add_subplot(211), - axes_2=figs[26].add_subplot(212), + axes_1=figs[27].add_subplot(211), + axes_2=figs[27].add_subplot(212), ) else: - ax19 = figs[24].add_subplot(211, aspect="equal") + ax19 = figs[25].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[24]) - plot_resistive_tf_info(ax19, m_file, scan, figs[24]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[25]) + plot_resistive_tf_info(ax19, m_file, scan, figs[25]) plot_tf_coil_structure( - figs[27].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[28].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[28], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[29], m_file, scan) - plot_tf_stress(figs[29].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[30].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[30].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[31].add_subplot(111), m_file, scan) - plot_cs_stress_time_profile(axis=figs[31].add_subplot(311), mfile=m_file, scan=scan) + plot_cs_stress_time_profile(axis=figs[32].add_subplot(311), mfile=m_file, scan=scan) plot_cs_coil_structure( - figs[31].add_subplot(223, aspect="equal"), figs[31], m_file, scan + figs[32].add_subplot(223, aspect="equal"), figs[32], m_file, scan ) plot_cs_turn_structure( - figs[31].add_subplot(326, aspect="equal"), figs[31], m_file, scan + figs[32].add_subplot(326, aspect="equal"), figs[32], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[32].add_subplot(221, aspect="equal"), m_file, scan + figs[33].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[32].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[32].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[33].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[33].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[33], m_file, scan) - ax_blanket = figs[33].add_subplot(122, aspect="equal") - plot_blkt_structure(ax_blanket, figs[33], m_file, scan, radial_build, colour_scheme) + plot_blkt_pipe_bends(figs[34], m_file, scan) + ax_blanket = figs[34].add_subplot(122, aspect="equal") + plot_blkt_structure(ax_blanket, figs[34], m_file, scan, radial_build, colour_scheme) plot_main_power_flow( - figs[34].add_subplot(111, aspect="equal"), m_file, scan, figs[34] + figs[35].add_subplot(111, aspect="equal"), m_file, scan, figs[35] ) - ax24 = figs[35].add_subplot(111) + ax24 = figs[36].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[35]) - - - - ax25 = figs[36].add_subplot(221) - PlasmaExhaust().plot_tritium_flow_contour(ax25, m_file, scan) - ax26 = figs[36].add_subplot(222) - PlasmaExhaust().plot_deuterium_flow_contour(ax26, m_file, scan) - ax27 = figs[36].add_subplot(223) - PlasmaExhaust().plot_alpha_flow_contour(ax27, m_file, scan) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[36]) def create_thickness_builds(m_file, scan: int): @@ -15427,7 +15418,7 @@ def plot_summary( ): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(36)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(37)] # run main_plot main_plot( diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index a9afe720b8..50038ad8a9 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1820,7 +1820,6 @@ def constraint_equation_94(constraint_registration): """ Deuterium particle balance """ - numerator = ( ( data_structure.physics_variables.eta_plasma_fuelling @@ -1848,7 +1847,6 @@ def constraint_equation_95(constraint_registration): """ Alpha particle balance """ - # Alpha particle balance numerator = ( data_structure.physics_variables.fusrat_dt_total @@ -1882,60 +1880,6 @@ def constraint_equation_96(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, - ) - - -@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, data: DataStructure): """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 92e7aca6eb..2d7c283883 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -249,12 +249,8 @@ class IterationVariable: 174: IterationVariable("triang", "physics", 0.00, 1.00), 175: IterationVariable("kappa", "physics", 0.00, 10.00), 176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30), - 177: IterationVariable( - "f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0 - ), - 178: IterationVariable( - "eta_plasma_fuelling", "physics", 0.01, 1.0 - ), + 177: IterationVariable("f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0), + 178: IterationVariable("eta_plasma_fuelling", "physics", 0.01, 1.0), 179: IterationVariable( "molflow_plasma_fuelling_vv_injected", "physics", diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index a7cfb35b26..ec0b66f65b 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -99,293 +99,6 @@ def description(self): return self._description_ -class PROCESSRunMode(IntEnum): - """Enumeration of the available PROCESS run modes, which determine the behaviour - of the code in various places. This is controlled by the `ioptimz` variable - """ - - EVALUATION = (-2, "Evaluation mode (no optimisation)") - """In this mode, the code will not perform any optimisation, and will instead - simply evaluate the constraints for the given input parameters, which is useful - for testing and for evaluating the performance of a given design point without - trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd - root-finding method from `scipy.optimize`) to seek a *consistent* solution by - varying a subset of the iteration variables until the consistency constraints - (equality constraints whose residuals must be driven to zero) are simultaneously - satisfied; no figure-of-merit is optimised, and the solver simply tries to find - a root of the constraint-residual vector. - """ - OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") - """In this mode, the code will perform optimisation using the VMCON solver - (or a custom solver if specified) to try to find a design point that optimises - the figure of merit while satisfying the constraints. This is the default mode - of operation for PROCESS. - """ - - def __new__(cls, value: int, description: str): - """Create a new PROCESSRunMode enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this run mode. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this run mode.""" - return self._description_ - - -class FiguresOfMerit(IntEnum): - """Enumeration of the available figures of merit (FoM) that can be used as - objective functions for optimisation in PROCESS. - """ - - MAJOR_RADIUS = (1, "Plasma major radius (R₀)") - NEUTRON_WALL_LOAD = (3, "Neutron wall load") - P_TF_PLUS_P_PF = (4, "TF & PF coil power") - FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") - COST_OF_ELECTRICITY = (6, "Cost of electricity") - CAPITAL_COST = (7, "Plant capital cost") - ASPECT_RATIO = (8, "Plasma aspect ratio") - DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") - TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") - TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") - PULSE_LENGTH = (14, "Pulse length") - PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") - MIN_R0_MAX_TAU_BURN = ( - 16, - "Linear combination of major radius (minimised) and pulse length (maximised)", - ) - NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") - NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") - MAX_Q_MAX_T_PLANT_PULSE_BURN = ( - 19, - "Linear combination of big Q and pulse length (maximised)", - ) - - def __new__(cls, value: int, description: str): - """Create a new FiguresOfMerit enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this figure of merit. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this figure of merit.""" - return self._description_ - - - -class PROCESSRunMode(IntEnum): - """Enumeration of the available PROCESS run modes, which determine the behaviour - of the code in various places. This is controlled by the `ioptimz` variable - """ - - EVALUATION = (-2, "Evaluation mode (no optimisation)") - """In this mode, the code will not perform any optimisation, and will instead - simply evaluate the constraints for the given input parameters, which is useful - for testing and for evaluating the performance of a given design point without - trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd - root-finding method from `scipy.optimize`) to seek a *consistent* solution by - varying a subset of the iteration variables until the consistency constraints - (equality constraints whose residuals must be driven to zero) are simultaneously - satisfied; no figure-of-merit is optimised, and the solver simply tries to find - a root of the constraint-residual vector. - """ - OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") - """In this mode, the code will perform optimisation using the VMCON solver - (or a custom solver if specified) to try to find a design point that optimises - the figure of merit while satisfying the constraints. This is the default mode - of operation for PROCESS. - """ - - def __new__(cls, value: int, description: str): - """Create a new PROCESSRunMode enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this run mode. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this run mode.""" - return self._description_ - - -class FiguresOfMerit(IntEnum): - """Enumeration of the available figures of merit (FoM) that can be used as - objective functions for optimisation in PROCESS. - """ - - MAJOR_RADIUS = (1, "Plasma major radius (R₀)") - NEUTRON_WALL_LOAD = (3, "Neutron wall load") - P_TF_PLUS_P_PF = (4, "TF & PF coil power") - FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") - COST_OF_ELECTRICITY = (6, "Cost of electricity") - CAPITAL_COST = (7, "Plant capital cost") - ASPECT_RATIO = (8, "Plasma aspect ratio") - DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") - TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") - TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") - PULSE_LENGTH = (14, "Pulse length") - PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") - MIN_R0_MAX_TAU_BURN = ( - 16, - "Linear combination of major radius (minimised) and pulse length (maximised)", - ) - NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") - NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") - MAX_Q_MAX_T_PLANT_PULSE_BURN = ( - 19, - "Linear combination of big Q and pulse length (maximised)", - ) - - def __new__(cls, value: int, description: str): - """Create a new FiguresOfMerit enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this figure of merit. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this figure of merit.""" - return self._description_ - - -class PROCESSRunMode(IntEnum): - """Enumeration of the available PROCESS run modes, which determine the behaviour - of the code in various places. This is controlled by the `ioptimz` variable - """ - - EVALUATION = (-2, "Evaluation mode (no optimisation)") - """In this mode, the code will not perform any optimisation, and will instead - simply evaluate the constraints for the given input parameters, which is useful - for testing and for evaluating the performance of a given design point without - trying to optimise it. Internally, PROCESS uses `fsolve` (a Newton-Krylov/hybrd - root-finding method from `scipy.optimize`) to seek a *consistent* solution by - varying a subset of the iteration variables until the consistency constraints - (equality constraints whose residuals must be driven to zero) are simultaneously - satisfied; no figure-of-merit is optimised, and the solver simply tries to find - a root of the constraint-residual vector. - """ - OPTIMISATION = (1, "Optimisation mode (e.g. via VMCON)") - """In this mode, the code will perform optimisation using the VMCON solver - (or a custom solver if specified) to try to find a design point that optimises - the figure of merit while satisfying the constraints. This is the default mode - of operation for PROCESS. - """ - - def __new__(cls, value: int, description: str): - """Create a new PROCESSRunMode enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this run mode. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this run mode.""" - return self._description_ - - -class FiguresOfMerit(IntEnum): - """Enumeration of the available figures of merit (FoM) that can be used as - objective functions for optimisation in PROCESS. - """ - - MAJOR_RADIUS = (1, "Plasma major radius (R₀)") - NEUTRON_WALL_LOAD = (3, "Neutron wall load") - P_TF_PLUS_P_PF = (4, "TF & PF coil power") - FUSION_GAIN_Q = (5, "Fusion gain (Qₚₗₐₛₘₐ)") - COST_OF_ELECTRICITY = (6, "Cost of electricity") - CAPITAL_COST = (7, "Plant capital cost") - ASPECT_RATIO = (8, "Plasma aspect ratio") - DIVERTOR_HEAT_LOAD = (9, "Divertor heat load") - TOROIDAL_FIELD = (10, "Plasma toroidal field on axis (B₀)") - TOTAL_INJECTED_POWER = (11, "Plasma total injected power (Pᵢₙⱼ)") - PULSE_LENGTH = (14, "Pulse length") - PLANT_AVAILABILITY_FACTOR = (15, "Plant availability factor") - MIN_R0_MAX_TAU_BURN = ( - 16, - "Linear combination of major radius (minimised) and pulse length (maximised)", - ) - NET_ELECTRICAL_OUTPUT = (17, "Plant net electrical output") - NULL_FIGURE_OF_MERIT = (18, "Null Figure of Merit") - MAX_Q_MAX_T_PLANT_PULSE_BURN = ( - 19, - "Linear combination of big Q and pulse length (maximised)", - ) - - def __new__(cls, value: int, description: str): - """Create a new FiguresOfMerit enum member with description. - - Args: - value: The integer value of the enum member. - description: The description for this figure of merit. - - Returns - ------- - The new enum member with attached description. - """ - obj = int.__new__(cls, value) - obj._value_ = value - obj._description_ = description - return obj - - @DynamicClassAttribute - def description(self): - """Return the description for this figure of merit.""" - return self._description_ - - - ipnvars: int = 183 """total number of variables available for iteration""" diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 1caa85f1c4..5fad3c4528 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -209,12 +209,11 @@ class PhysicsData: f_plasma_fuel_burnup: float = 0.0 """Total fuel burnup fraction in plasma""" + f_plasma_tritium_burnup: float = None + """Tritium 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""" + f_plasma_deuterium_burnup: float = None + """Deuterium burnup fraction in plasma""" burnup_in: float = 0.0 """fractional plasma burnup user input""" @@ -373,41 +372,27 @@ class PhysicsData: fusrat_total: float = 0.0 """fusion reaction rate, from beams and plasma (reactions/sec)""" - fusrat_plasma_dt: float = None + fusrat_plasma_dt: float = 0.0 """ D-T fusion reaction rate in plasma, (reactions/sec)""" - fusrat_dt_total: float = None + fusrat_dt_total: float = 0.0 """ Total D-T fusion reaction rate from beams and plasma, (reactions/sec)""" - fusrat_plasma_dd_helion: float = None + fusrat_plasma_dd_helion: float = 0.0 """D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" - fusrat_plasma_dd_triton: float = None + fusrat_plasma_dd_triton: float = 0.0 """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_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_dhe3: float = None -"""D-3He fusion reaction rate in plasma, (reactions/sec)""" + fusrat_plasma_dd_total: float = 0.0 + """Total D-D 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_dhe3: float = 0.0 + """D-3He fusion reaction rate in 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_neutron_production_total: float = 0.0 + """Total neutron production rate from plasma and beams (neutrons/sec)""" - fusrat_plasma_dt_profile: list[float] = field(default_factory=list) """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -914,33 +899,32 @@ class PhysicsData: in which case q95 = mean edge safety factor qbar) """ + f_plasma_particles_lcfs_recycled: float = 0.9 + """fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)""" -f_plasma_particles_lcfs_recycled: float = None -"""fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)""" + eta_plasma_fuelling: float = 0.7 + """fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)""" -eta_plasma_fuelling: float = None -"""fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)""" + molflow_plasma_fuelling_vv_injected: float = 1e21 + """plasma fuelling rate into the vacuum vessel (particles/s)""" -molflow_plasma_fuelling_vv_injected: float = None -"""plasma fuelling rate into the vacuum vessel (particles/s)""" + molflow_plasma_fuelling_vv_injected_moles: float = 0.0 + """plasma fuelling rate into the vacuum vessel (moles/s)""" -molflow_plasma_fuelling_vv_injected_moles: float = None -"""plasma fuelling rate into the vacuum vessel (moles/s)""" + molflow_plasma_fuelling_loss: float = 0.0 + """plasma fuelling rate that dosent make it to plasma (particles/s)""" -molflow_plasma_fuelling_loss: float = None -"""plasma fuelling rate that dosent make it to plasma (particles/s)""" + molflow_plasma_fuelling_loss_moles: float = 0.0 + """plasma fuelling rate that dosent make it to plasma (moles/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 = 0.5 + """fraction of plasma fuelling that is deuterium""" -f_molflow_plasma_fuelling_deuterium: float = None -"""fraction of plasma fuelling that is deuterium""" + f_molflow_plasma_fuelling_tritium: float = 0.5 + """fraction of plasma fuelling that is tritium""" -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""" + f_molflow_plasma_fuelling_helium3: float = 0.0 + """fraction of plasma fuelling that is helium-3""" q95_min: float = 0.0 """lower limit for edge safety factor""" @@ -984,9 +968,8 @@ class PhysicsData: f_nd_beam_electron: float = 0.005 """hot beam density / n_e (`iteration variable 7`)""" - -f_nd_plasma_carbon_electron: float = None -"""n_carbon / n_e""" + f_nd_plasma_carbon_electron: float = 0.0 + """n_carbon / n_e""" f_nd_plasma_iron_argon_electron: float = 0.0 """n_highZ / n_e""" diff --git a/process/main.py b/process/main.py index 788f2ece71..b364e8be5b 100644 --- a/process/main.py +++ b/process/main.py @@ -766,6 +766,7 @@ def models(self) -> tuple[Model, ...]: self.plasma_fields, self.sauter_bootstrap_current, self.plasma_transition, + self.plasma_fuelling, self.physics_detailed, self.electron_cyclotron, self.lower_hybrid, diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py index 35e8d4c6e8..3c7640baf9 100644 --- a/process/models/physics/fuelling.py +++ b/process/models/physics/fuelling.py @@ -6,54 +6,60 @@ import process.core.io.mfile as mf from process.core import constants from process.core import process_output as po +from process.core.model import Model from process.data_structure import ( numerics, - physics_variables, reinke_variables, ) logger = logging.getLogger(__name__) -class PlasmaFuelling: +class PlasmaFuelling(Model): """Class to hold plasma fuelling calculations for plasma processing.""" 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 + self.data.physics.molflow_plasma_fuelling_vv_injected_moles = ( + self.data.physics.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) + self.data.physics.molflow_plasma_fuelling_loss = ( + self.data.physics.molflow_plasma_fuelling_vv_injected + * (1 - self.data.physics.eta_plasma_fuelling) ) - physics_variables.molflow_plasma_fuelling_loss_moles = ( - physics_variables.molflow_plasma_fuelling_loss / constants.AVOGADRO_NUMBER + self.data.physics.molflow_plasma_fuelling_loss_moles = ( + self.data.physics.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, + self.data.physics.f_plasma_fuel_burnup = self.calculate_fuel_burnup_fraction( + fusrat_total=self.data.physics.fusrat_total, + molflow_plasma_fuelling_vv_injected=self.data.physics.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, + self.data.physics.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction( + fusrat_dt_total=self.data.physics.fusrat_dt_total, + molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_tritium=self.data.physics.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, + self.data.physics.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction( + fusrat_plasma_dd_total=self.data.physics.fusrat_plasma_dd_total, + molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_deuterium=self.data.physics.f_molflow_plasma_fuelling_deuterium, + fusrat_dt_total=self.data.physics.fusrat_dt_total, + fusrat_plasma_dhe3=self.data.physics.fusrat_plasma_dhe3, ) + def output(self): + """This model doesn't output to the output file, but it does generate contour + plots of plasma fuel flow rates vs recycling and fuelling efficiency. + """ + @staticmethod def calculate_fuel_burnup_fraction( fusrat_total: float, molflow_plasma_fuelling_vv_injected: float @@ -77,7 +83,6 @@ def calculate_fuel_burnup_fraction( 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 @@ -106,7 +111,6 @@ def calculate_tritium_burnup_fraction( 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 ) @@ -143,7 +147,6 @@ def calculate_deuterium_burnup_fraction( 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 + 2 * fusrat_plasma_dd_total + fusrat_plasma_dhe3) / ( molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium ) @@ -192,7 +195,6 @@ def calculate_plasma_tritium_flow_rate( Tritium flow rate in the plasma exhaust (particles/s). """ - return ( ( f_molflow_plasma_fuelling_tritium @@ -311,7 +313,6 @@ def calculate_plasma_helium3_flow_rate( Helium-3 flow rate in the plasma exhaust (particles/s). """ - return ( ( f_molflow_plasma_fuelling_helium3 @@ -357,7 +358,6 @@ def calculate_plasma_alphas_flow_rate( Alpha particle flow rate in the plasma exhaust (particles/s). """ - # Alpha particle balance return ( @@ -369,7 +369,6 @@ def calculate_plasma_alphas_flow_rate( 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))) @@ -439,7 +438,6 @@ def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): 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))) @@ -510,7 +508,6 @@ 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.""" - 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(( @@ -576,7 +573,6 @@ def plot_alpha_flow_contour(self, axis: plt.Axes, mfile: mf.MFile, scan: int): 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))) @@ -684,91 +680,91 @@ def output_fuelling_info(self): self.outfile, "Fuelling rate (nucleus-pairs/s)", "(molflow_plasma_fuelling_vv_injected)", - physics_variables.molflow_plasma_fuelling_vv_injected, + self.data.physics.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, + self.data.physics.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, + self.data.physics.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, + self.data.physics.molflow_plasma_fuelling_loss_moles, "OP ", ) po.ovarre( self.outfile, "Fraction of plasma fuelling that is deuterium", "(f_molflow_plasma_fuelling_deuterium)", - physics_variables.f_molflow_plasma_fuelling_deuterium, + self.data.physics.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, + self.data.physics.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, + self.data.physics.f_molflow_plasma_fuelling_helium3, "OP ", ) po.ovarre( self.outfile, "Fuelling efficiency", "(eta_plasma_fuelling)", - physics_variables.eta_plasma_fuelling, + self.data.physics.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, + self.data.physics.f_plasma_particles_lcfs_recycled, "OP ", ) po.ovarre( self.outfile, "Fuel burn-up rate (reactions/s)", "(fusrat_total)", - physics_variables.fusrat_total, + self.data.physics.fusrat_total, "OP ", ) po.ovarrf( self.outfile, "Total fuel burn-up fraction", "(f_plasma_fuel_burnup)", - physics_variables.f_plasma_fuel_burnup, + self.data.physics.f_plasma_fuel_burnup, "OP ", ) po.ovarrf( self.outfile, "Tritium burn-up fraction", "(f_plasma_tritium_burnup)", - physics_variables.f_plasma_tritium_burnup, + self.data.physics.f_plasma_tritium_burnup, "OP ", ) po.ovarrf( self.outfile, "Deuterium burn-up fraction", "(f_plasma_deuterium_burnup)", - physics_variables.f_plasma_deuterium_burnup, + self.data.physics.f_plasma_deuterium_burnup, "OP ", ) diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 19512707b1..5cd2399c04 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -347,8 +347,8 @@ 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 + self.data.physics.fusrat_plasma_dhe3 = ( + fusion_rate_density * self.data.physics.vol_plasma ) # Update the cumulative D-3He power density @@ -455,8 +455,8 @@ 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 + self.data.physics.fusrat_plasma_dd_helion = ( + fusion_rate_density * self.data.physics.vol_plasma ) # Update the cumulative D-D power density @@ -559,8 +559,8 @@ def dd_triton_reaction(self): # Proton production rate [particles/m³/s] proton_rate_density = fusion_rate_density - physics_variables.fusrat_plasma_dd_triton = ( - fusion_rate_density * physics_variables.vol_plasma + self.data.physics.fusrat_plasma_dd_triton = ( + fusion_rate_density * self.data.physics.vol_plasma ) # Update the cumulative D-D power density diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 3f95e0e755..fb4ea0cf9b 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -32,7 +32,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.l_h_transition import PlasmaConfinementTransition from process.models.physics.plasma_current import ( PlasmaCurrent, PlasmaDiamagneticCurrent, @@ -648,20 +648,20 @@ def run(self): self.data.physics.fusrat_total = ( self.data.physics.fusden_total * self.data.physics.vol_plasma ) - physics_variables.fusrat_plasma_dt = (physics_variables.p_plasma_dt_mw * 1e6) / ( + self.data.physics.fusrat_plasma_dt = (self.data.physics.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 + self.data.physics.fusrat_plasma_dd_total = ( + self.data.physics.fusrat_plasma_dd_helion + + self.data.physics.fusrat_plasma_dd_triton ) - physics_variables.fusrat_neutron_production_total = ( - physics_variables.fusrat_plasma_dd_helion + physics_variables.fusrat_dt_total + self.data.physics.fusrat_neutron_production_total = ( + self.data.physics.fusrat_plasma_dd_helion + self.data.physics.fusrat_dt_total ) - physics_variables.fusrat_dt_total = ( - physics_variables.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) + self.data.physics.fusrat_dt_total = ( + self.data.physics.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) ) # Create some derived values and add beam contribution to fusion power @@ -915,22 +915,22 @@ def run(self): self.data.physics.ind_plasma_internal_norm, ) - physics_variables.e_plasma_magnetic_stored = ( - 0.5e0 * physics_variables.ind_plasma * physics_variables.plasma_current**2 + self.data.physics.e_plasma_magnetic_stored = ( + 0.5e0 * self.data.physics.ind_plasma * self.data.physics.plasma_current**2 ) ( - physics_variables.t_alpha_confinement, - physics_variables.f_alpha_energy_confinement, + self.data.physics.t_alpha_confinement, + self.data.physics.f_alpha_energy_confinement, ) = self.phyaux( - physics_variables.fusden_alpha_total, - physics_variables.nd_plasma_alphas_vol_avg, - physics_variables.t_energy_confinement, + self.data.physics.fusden_alpha_total, + self.data.physics.nd_plasma_alphas_vol_avg, + self.data.physics.t_energy_confinement, ) self.fuelling.run() - physics_variables.ntau, physics_variables.nTtau = ( + self.data.physics.ntau, self.data.physics.nTtau = ( self.confinement.calculate_double_and_triple_product( nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, t_energy_confinement=self.data.physics.t_energy_confinement, @@ -1418,8 +1418,6 @@ def phyaux( fusden_alpha_total: float, nd_plasma_alphas_vol_avg: float, t_energy_confinement: float, - burnup_in: float, - tauratio: float, ) -> tuple[float, float]: """Auxiliary physics quantities @@ -1442,7 +1440,6 @@ def phyaux( parts of the code. """ - # Alpha particle confinement time (s) # Number of alphas / alpha production rate # only likely if DD is only active fusion reaction @@ -1735,14 +1732,14 @@ def outplas(self): self.outfile, "D-T Fusion rate: total (reactions/sec)", "(fusrat_dt_total)", - physics_variables.fusrat_dt_total, + self.data.physics.fusrat_dt_total, "OP ", ) po.ovarre( self.outfile, "D-T Fusion rate: plasma (reactions/sec)", "(fusrat_plasma_dt)", - physics_variables.fusrat_plasma_dt, + self.data.physics.fusrat_plasma_dt, "OP ", ) @@ -1750,56 +1747,56 @@ def outplas(self): self.outfile, "D-D -> 3He Fusion rate: plasma (reactions/sec)", "(fusrat_plasma_dd_helion)", - physics_variables.fusrat_plasma_dd_helion, + self.data.physics.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, + self.data.physics.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, + self.data.physics.fusrat_plasma_dd_total, "OP ", ) po.ovarre( self.outfile, "D-3He Fusion rate: total (reactions/sec)", "(fusrat_plasma_dhe3)", - physics_variables.fusrat_plasma_dhe3, + self.data.physics.fusrat_plasma_dhe3, "OP ", ) po.ovarre( self.outfile, "Neutron production rate: total (particles/sec)", "(fusrat_neutron_production_total)", - physics_variables.fusrat_neutron_production_total, + self.data.physics.fusrat_neutron_production_total, "OP ", ) po.ovarre( self.outfile, "D-3He Fusion rate: total (reactions/sec)", "(fusrat_plasma_dhe3)", - physics_variables.fusrat_plasma_dhe3, + self.data.physics.fusrat_plasma_dhe3, "OP ", ) po.ovarre( self.outfile, "Neutron production rate: total (particles/sec)", "(fusrat_neutron_production_total)", - physics_variables.fusrat_neutron_production_total, + self.data.physics.fusrat_neutron_production_total, "OP ", ) po.ovarre( self.outfile, "D-D Fusion rate: total (reactions/sec)", "(fusrat_plasma_dd_total)", - physics_variables.fusrat_plasma_dd_total, + self.data.physics.fusrat_plasma_dd_total, "OP ", ) po.ovarre( @@ -2331,6 +2328,525 @@ def outplas(self): self.fuelling.output_fuelling_info() + def output_temperature_density_profile_info(self) -> None: + """Output information about plasma temperature and density profiles.""" + po.oheadr(self.outfile, "Plasma Density and Temperature Profiles") + po.ovarrf( + self.outfile, + "Number of radial points in plasma profiles", + "(n_plasma_profile_elements)", + self.data.physics.n_plasma_profile_elements, + ) + po.ovarin( + self.outfile, + "Plasma profile model selected", + "(i_plasma_pedestal)", + self.data.physics.i_plasma_pedestal, + ) + po.ocmmnt( + self.outfile, + f"Plasma profile model selected: " + f"{PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal).description}", + ) + + if ( + PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) + == PlasmaProfileShapeType.PEDESTAL_PROFILE + ) and ( + self.data.physics.nd_plasma_electron_on_axis + < self.data.physics.nd_plasma_pedestal_electron + ): + logger.error("Central density is less than pedestal density") + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.osubhd(self.outfile, "Temperature:") + po.ovarrf( + self.outfile, + "Temperature profile index (αₜ)", + "(alphat)", + self.data.physics.alphat, + ) + po.ovarrf( + self.outfile, + "Temperature profile index beta (βₜ)", + "(tbeta)", + self.data.physics.tbeta, + ) + po.oblnkl(self.outfile) + + if ( + PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) + == PlasmaProfileShapeType.PEDESTAL_PROFILE + ): + po.ovarrf( + self.outfile, + "Temperature pedestal r/a location (ρₜ,pedestal)", + "(radius_plasma_pedestal_temp_norm)", + self.data.physics.radius_plasma_pedestal_temp_norm, + ) + + po.ovarrf( + self.outfile, + "Electron temperature pedestal height (Tₑ,pedestal) (keV)", + "(temp_plasma_pedestal_kev)", + self.data.physics.temp_plasma_pedestal_kev, + ) + if 78 in numerics.icc: + po.ovarrf( + self.outfile, + "Electron temperature at separatrix (Tₑ,ₛₑₚ) (keV)", + "(temp_plasma_separatrix_kev)", + self.data.physics.temp_plasma_separatrix_kev, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Electron temperature at separatrix (Tₑ,ₛₑₚ) (keV)", + "(temp_plasma_separatrix_kev)", + self.data.physics.temp_plasma_separatrix_kev, + ) + po.oblnkl(self.outfile) + + po.ovarrf( + self.outfile, + "Volume averaged electron temperature (⟨Tₑ⟩) (keV)", + "(temp_plasma_electron_vol_avg_kev)", + self.data.physics.temp_plasma_electron_vol_avg_kev, + ) + po.ovarrf( + self.outfile, + "Electron temperature on axis (Tₑ₀) (keV)", + "(temp_plasma_electron_on_axis_kev)", + self.data.physics.temp_plasma_electron_on_axis_kev, + "OP ", + ) + po.ovarrf( + self.outfile, + "Volume averaged density weighted electron temperature (⟨Tₑ⟩ₙ) (keV)", + "(temp_plasma_electron_density_weighted_kev)", + self.data.physics.temp_plasma_electron_density_weighted_kev, + "OP ", + ) + po.ovarrf( + self.outfile, + "Ratio of electron density weighted to volume averaged temperature", + "(f_temp_plasma_electron_density_vol_avg)", + self.data.physics.f_temp_plasma_electron_density_vol_avg, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarrf( + self.outfile, + "Ratio of ion to electron volume-averaged temperature", + "(f_temp_plasma_ion_electron)", + self.data.physics.f_temp_plasma_ion_electron, + "IP ", + ) + po.oblnkl(self.outfile) + po.ovarrf( + self.outfile, + "Volume averaged ion temperature (⟨Tᵢ⟩) (keV)", + "(temp_plasma_ion_vol_avg_kev)", + self.data.physics.temp_plasma_ion_vol_avg_kev, + ) + po.ovarrf( + self.outfile, + "Ion temperature on axis (Tᵢ₀) (keV)", + "(temp_plasma_ion_on_axis_kev)", + self.data.physics.temp_plasma_ion_on_axis_kev, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + + po.osubhd(self.outfile, "Density:") + po.ovarrf( + self.outfile, + "Density profile factor (αₙ)", + "(alphan)", + self.data.physics.alphan, + ) + po.oblnkl(self.outfile) + if ( + PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) + == PlasmaProfileShapeType.PEDESTAL_PROFILE + ): + po.ovarrf( + self.outfile, + "Density pedestal r/a location (ρₙ,pedestal)", + "(radius_plasma_pedestal_density_norm)", + self.data.physics.radius_plasma_pedestal_density_norm, + ) + if self.data.physics.f_nd_plasma_pedestal_greenwald >= 0e0: + po.ovarre( + self.outfile, + "Electron density pedestal height (nₑ_pedestal) (/m³)", + "(nd_plasma_pedestal_electron)", + self.data.physics.nd_plasma_pedestal_electron, + "OP ", + ) + else: + po.ovarre( + self.outfile, + "Electron density pedestal height (nₑ_pedestal) (/m³)", + "(nd_plasma_pedestal_electron)", + self.data.physics.nd_plasma_pedestal_electron, + ) + # must be assigned to their exisiting values# + fgwped_out = ( + self.data.physics.nd_plasma_pedestal_electron + / self.data.physics.nd_plasma_electron_max_array[6] + ) + fgwsep_out = ( + self.data.physics.nd_plasma_separatrix_electron + / self.data.physics.nd_plasma_electron_max_array[6] + ) + if self.data.physics.f_nd_plasma_pedestal_greenwald >= 0e0: + self.data.physics.f_nd_plasma_pedestal_greenwald = ( + self.data.physics.nd_plasma_pedestal_electron + / self.data.physics.nd_plasma_electron_max_array[6] + ) + if self.data.physics.f_nd_plasma_separatrix_greenwald >= 0e0: + self.data.physics.f_nd_plasma_separatrix_greenwald = ( + self.data.physics.nd_plasma_separatrix_electron + / self.data.physics.nd_plasma_electron_max_array[6] + ) + + po.ovarre( + self.outfile, + "Pedestal Greenwald fraction", + "(fgwped_out)", + fgwped_out, + ) + po.ovarre( + self.outfile, + "Electron density at separatrix (nₑ,ₛₑₚ) (/m³)", + "(nd_plasma_separatrix_electron)", + self.data.physics.nd_plasma_separatrix_electron, + ) + po.ovarre( + self.outfile, + "Separatrix Greenwald fraction", + "(fgwsep_out)", + fgwsep_out, + ) + po.oblnkl(self.outfile) + + po.ovarre( + self.outfile, + "Volume averaged electron number density (⟨nₑ⟩) (/m³)", + "(nd_plasma_electrons_vol_avg)", + self.data.physics.nd_plasma_electrons_vol_avg, + ) + po.ovarre( + self.outfile, + "Electron number density on axis (nₑ₀) (/m³)", + "(nd_plasma_electron_on_axis)", + self.data.physics.nd_plasma_electron_on_axis, + "OP ", + ) + po.ovarre( + self.outfile, + "Line-averaged electron number density (ñₑ) (/m³)", + "(nd_plasma_electron_line)", + self.data.physics.nd_plasma_electron_line, + "OP ", + ) + if self.data.stellarator.istell == 0: + po.ovarre( + self.outfile, + "Greenwald fraction (f_GW)", + "(dnla_gw)", + self.data.physics.nd_plasma_electron_line + / self.data.physics.nd_plasma_electron_max_array[6], + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total ion volume averaged number density (⟨nᵢ⟩) (/m³)", + "(nd_plasma_ions_total_vol_avg)", + self.data.physics.nd_plasma_ions_total_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuel ion volume averaged number density (⟨n_fuel⟩) (/m³)", + "(nd_plasma_fuel_ions_vol_avg)", + self.data.physics.nd_plasma_fuel_ions_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Total impurity volume averaged number density with Z > 2 (⟨nᵢₘₚ⟩) (/m³)", + "(nd_plasma_impurities_vol_avg)", + self.data.physics.nd_plasma_impurities_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermalised alpha volume averaged number density (⟨n_alpha⟩) (/m³)", + "(nd_plasma_alphas_vol_avg)", + self.data.physics.nd_plasma_alphas_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermalised alpha to electron number density ratio (⟨n_alpha⟩/⟨nₑ⟩)", + "(f_nd_alpha_electron)", + self.data.physics.f_nd_alpha_electron, + ) + po.ovarre( + self.outfile, + "Proton volume averaged number density (⟨nₚ⟩) (/m³)", + "(nd_plasma_protons_vol_avg)", + self.data.physics.nd_plasma_protons_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Proton to electron volume averaged number density ratio (⟨nₚ⟩/⟨nₑ⟩)", + "(f_nd_protium_electrons)", + self.data.physics.f_nd_protium_electrons, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Hot beam ion volume averaged number density (⟨n_beam⟩) (/m³)", + "(nd_beam_ions)", + self.data.physics.nd_beam_ions, + "OP ", + ) + po.ovarre( + self.outfile, + "Hot beam ion to electron number density ratio (⟨n_beam⟩/⟨nₑ⟩)", + "(f_nd_beam_electron)", + self.data.physics.f_nd_beam_electron, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + + po.osubhd(self.outfile, "Pressure:") + + po.ovarrf( + self.outfile, + "Pressure profile index (αₚ)", + "(alphap)", + self.data.physics.alphap, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Plasma thermal pressure on axis (p₀) (Pa)", + "(pres_plasma_thermal_on_axis)", + self.data.physics.pres_plasma_thermal_on_axis, + "OP ", + ) + po.ovarre( + self.outfile, + "Volume averaged plasma thermal pressure (⟨p⟩) (Pa)", + "(pres_plasma_thermal_vol_avg)", + self.data.physics.pres_plasma_thermal_vol_avg, + "OP ", + ) + # As array output is not currently supported, each element is output as a float + # instance + # Output plasma pressure profiles to mfile + for i in range(len(self.data.physics.pres_plasma_thermal_total_profile)): + po.ovarre( + self.mfile, + f"Total plasma pressure at point {i}", + f"(pres_plasma_thermal_total_profile{i})", + self.data.physics.pres_plasma_thermal_total_profile[i], + ) + for i in range(len(self.data.physics.pres_plasma_electron_profile)): + po.ovarre( + self.mfile, + f"Total plasma electron pressure at point {i}", + f"(pres_plasma_electron_profile{i})", + self.data.physics.pres_plasma_electron_profile[i], + ) + for i in range(len(self.data.physics.pres_plasma_ion_total_profile)): + po.ovarre( + self.mfile, + f"Total plasma ion pressure at point {i}", + f"(pres_plasma_ion_total_profile{i})", + self.data.physics.pres_plasma_ion_total_profile[i], + ) + for i in range(len(self.data.physics.pres_plasma_fuel_profile)): + po.ovarre( + self.mfile, + f"Total plasma fuel pressure at point {i}", + f"(pres_plasma_fuel_profile{i})", + self.data.physics.pres_plasma_fuel_profile[i], + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.oblnkl(self.outfile) + + po.ocmmnt(self.outfile, "Impurities:") + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Plasma ion densities / electron density:") + po.oblnkl(self.outfile) + + for imp in range(N_IMPURITIES): + # MDK Update f_nd_impurity_electrons, as this will make the ITV output + # work correctly. + self.data.impurity_radiation.f_nd_impurity_electrons[imp] = ( + self.data.impurity_radiation.f_nd_impurity_electron_array[imp] + ) + str1 = ( + self.data.impurity_radiation + .impurity_arr_label[imp] + .item() + .replace("_", "") + + " concentration" + ) + str2 = f"(f_nd_impurity_electrons({imp + 1:02}))" + # MDK Add output flag for H which is calculated. + if imp == 0: + po.ovarre( + self.outfile, + str1, + str2, + self.data.impurity_radiation.f_nd_impurity_electrons[imp], + "OP ", + ) + else: + po.ovarre( + self.outfile, + str1, + str2, + self.data.impurity_radiation.f_nd_impurity_electrons[imp], + ) + if self.data.impurity_radiation.f_nd_impurity_electrons[imp] != 0.0: # noqa: RUF069 + for i in range(self.data.physics.n_plasma_profile_elements): + po.ovarre( + self.mfile, + str1 + f" at point {i}", + f"(f_nd_impurity_electrons{imp}_{i})", + self.data.impurity_radiation.f_nd_impurity_electrons[imp] + * self.plasma_profile.neprofile.profile_y[i], + "OP ", + ) + + po.oblnkl(self.outfile) + po.ovarrf( + self.outfile, + "Volume averaged plasma effective charge (⟨Zₑ⟩)", + "(n_charge_plasma_effective_vol_avg)", + self.data.physics.n_charge_plasma_effective_vol_avg, + "OP ", + ) + for i in range(len(self.data.physics.n_charge_plasma_effective_profile)): + po.ovarre( + self.mfile, + "Effective charge at point", + f"(n_charge_plasma_effective_profile{i})", + self.data.physics.n_charge_plasma_effective_profile[i], + "OP ", + ) + + for imp in range(N_IMPURITIES): + for i in range(self.data.physics.n_plasma_profile_elements): + po.ovarre( + self.mfile, + "Impurity charge at point", + f"(n_charge_plasma_profile{imp}_{i})", + self.data.impurity_radiation.n_charge_impurity_profile[imp][i], + "OP ", + ) + + po.ovarrf( + self.outfile, + "Volume averaged mass-weighted plasma effective charge (⟨Zₑ⟩ₘ)", + "(n_charge_plasma_effective_mass_weighted_vol_avg)", + self.data.physics.n_charge_plasma_effective_mass_weighted_vol_avg, + "OP ", + ) + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "----------------------------") + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Masses:") + po.oblnkl(self.outfile) + + po.ovarre( + self.outfile, + "Average mass of all ions (amu)", + "(m_ions_total_amu)", + self.data.physics.m_ions_total_amu, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of all ions in plasma (kg)", + "(m_plasma_ions_total)", + self.data.physics.m_plasma_ions_total, + "OP ", + ) + po.ovarre( + self.outfile, + "Average mass of all fuel ions (amu)", + "(m_fuel_amu)", + self.data.physics.m_fuel_amu, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of all fuel ions in plasma (kg)", + "(m_plasma_fuel_ions)", + self.data.physics.m_plasma_fuel_ions, + "OP ", + ) + po.ovarre( + self.outfile, + "Average mass of all beam ions (amu)", + "(m_beam_amu)", + self.data.physics.m_beam_amu, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of all alpha particles in plasma (kg)", + "(m_plasma_alpha)", + self.data.physics.m_plasma_alpha, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of all electrons in plasma (kg)", + "(m_plasma_electron)", + self.data.physics.m_plasma_electron, + "OP ", + ) + po.ovarre( + self.outfile, + "Total mass of the plasma (kg)", + "(m_plasma)", + self.data.physics.m_plasma, + "OP ", + ) + + # Issue 558 - addition of constraint 76 to limit the value of + # nd_plasma_separatrix_electron, in proportion with the ballooning parameter + # and Greenwald density + if 76 in numerics.icc: + po.ovarre( + self.outfile, + "Critical ballooning parameter value", + "(alpha_crit)", + self.data.physics.alpha_crit, + ) + po.ovarre( + self.outfile, + "Critical electron density at separatrix (/m3)", + "(nd_plasma_separatrix_electron_eich_max)", + self.data.physics.nd_plasma_separatrix_electron_eich_max, + ) + @staticmethod @nb.njit(cache=True) def calculate_plasma_masses( diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 7b31377c55..840e98d6ef 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2336,17 +2336,17 @@ def st_phys(self, output): # for the rest of the code ( - physics_variables.t_alpha_confinement, - physics_variables.f_alpha_energy_confinement, + self.data.physics.t_alpha_confinement, + self.data.physics.f_alpha_energy_confinement, ) = self.physics.phyaux( - physics_variables.fusden_alpha_total, - physics_variables.nd_plasma_alphas_vol_avg, - physics_variables.t_energy_confinement, + self.data.physics.fusden_alpha_total, + self.data.physics.nd_plasma_alphas_vol_avg, + self.data.physics.t_energy_confinement, ) - 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, + self.data.physics.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction( + fusrat_total=self.data.physics.fusrat_total, + molflow_plasma_fuelling_vv_injected=self.data.physics.molflow_plasma_fuelling_vv_injected, ) # Calculate the neoclassical sanity check with PROCESS parameters diff --git a/tests/unit/models/physics/test_physics.py b/tests/unit/models/physics/test_physics.py index 5fe4ed59c3..bf412ff2ad 100644 --- a/tests/unit/models/physics/test_physics.py +++ b/tests/unit/models/physics/test_physics.py @@ -6,7 +6,6 @@ import pytest from process.core import constants -from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics,