From 3fec599ef3a369901da71a66d24ffb5e7302bdd7 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:01:19 +0100 Subject: [PATCH 01/10] Rename NeProfile to ElectronDensityProfile and update references throughout the codebase --- .../physics-models/profiles/plasma_density_profile.md | 2 +- .../source/physics-models/profiles/plasma_profiles.md | 8 ++++---- process/main.py | 4 ++-- process/models/physics/plasma_profiles.py | 5 +++-- process/models/physics/profiles.py | 2 +- 5 files changed, 11 insertions(+), 10 deletions(-) diff --git a/documentation/source/physics-models/profiles/plasma_density_profile.md b/documentation/source/physics-models/profiles/plasma_density_profile.md index 97a118b04d..44516e6ca9 100644 --- a/documentation/source/physics-models/profiles/plasma_density_profile.md +++ b/documentation/source/physics-models/profiles/plasma_density_profile.md @@ -1,4 +1,4 @@ -# Density Profile | `NeProfile(Profile)` +# Density Profile | `ElectronDensityProfile(Profile)` The density profile class is organised around a central runner function that is called each time the plasma is parameterised by the parent [`PlasmaProfile()`](./plasma_profiles.md) class. It is called by [`pedestal_parameterisation()`](plasma_profiles.md#pedestal_parameterisation) and [`parabolic parameterisation()`](./plasma_profiles.md#parabolic_parameterisation). The sequence of the runner function can be seen below along with explanation of the following calculations. diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index 6032e13737..85bf35c275 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -4,7 +4,7 @@ In `PROCESS` the density, temperature and current profiles of the plasma for electrons and ions can take two forms depending on the switch value for `i_plasma_pedestal`. Either without a [pedestal](http://fusionwiki.ciemat.es/wiki/Pedestal), `i_plasma_pedestal == 0` or with a pedestal `i_plasma_pedestal == 1`. `i_plasma_pedestal == 0` is better suited for modelling L-mode plasmas, while `i_plasma_pedestal == 1` is better suited for modelling [H-mode](https://en.wikipedia.org/wiki/High-confinement_mode) plasmas. -The files responsible for calculating and storing the profiles are `plasma_profiles.py` and `profiles.py`. A central plasma profile object is created from the [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) class that contains attributes for the plasma density and temperature. The density and temperature profiles are in themselves objects of the [`Profile`](./plasma_profiles_abstract_class.md) abstract base class. [`Profile`](./plasma_profiles_abstract_class.md), [`NeProfile`](plasma_density_profile.md) and [`TeProfile`](./plasma_temperature_profile.md) are all defined in `profiles.py`. [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) is exclusively in `plasma_profiles.py` +The files responsible for calculating and storing the profiles are `plasma_profiles.py` and `profiles.py`. A central plasma profile object is created from the [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) class that contains attributes for the plasma density and temperature. The density and temperature profiles are in themselves objects of the [`Profile`](./plasma_profiles_abstract_class.md) abstract base class. [`Profile`](./plasma_profiles_abstract_class.md), [`ElectronDensityProfile`](plasma_density_profile.md) and [`TeProfile`](./plasma_temperature_profile.md) are all defined in `profiles.py`. [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) is exclusively in `plasma_profiles.py`
![UML of profiles](./uml_classes_PlasmaProfile.png){height="1000px"} @@ -238,7 +238,7 @@ The graph below is for a standard pedestal profile. You can vary its attributes ## Plasma Profile Class | `PlasmaProfile` ### Initialization | `__init__()` -The parent plasma profile class is `PlasmaProfile`. Initialization sets the profile class size and `neprofile` and `teprofile` to [`NeProfile`](plasma_density_profile.md) & [`TeProfile`](plasma_temperature_profile.md) objects from `profiles.py` +The parent plasma profile class is `PlasmaProfile`. Initialization sets the profile class size and `neprofile` and `teprofile` to [`ElectronDensityProfile`](plasma_density_profile.md) & [`TeProfile`](plasma_temperature_profile.md) objects from `profiles.py` ???+ Note @@ -266,7 +266,7 @@ Depending on the value of `i_plasma_pedestal` different functions will be ran, t If pedestal profile values are set they are reset to have values that agree with the original form of the parabolic profiles. Such that $\rho_{\text{ped}} = 1$ and that pedestal and separatrix densities and temperatures are zero. This will then warn the user in the terminal. -The density and temperature profile runner function [`TeProfile/NeProfile.run()`](plasma_density_profile.md#runner-function-run) is then called to re-calculate the profile and core values. +The density and temperature profile runner function [`TeProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is then called to re-calculate the profile and core values. Ratio of density-weighted to volume-averaged temperature factor is calculated: @@ -593,7 +593,7 @@ $$ ##### `pedestal_parameterisation()` -The density and temperature profile runner function [`TeProfile/NeProfile.run()`](plasma_density_profile.md#runner-function-run) is firstly called to re-calculate the profile and core values. +The density and temperature profile runner function [`TeProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is firstly called to re-calculate the profile and core values. Perform integrations to calculate ratio of density-weighted to volume-averaged temperature, etc. Density-weighted temperature = $\frac{\int{nT \ dV}}{\int{n \ dV}}$, which is approximately equal to the ratio $\frac{\int{\rho \ n(\rho) T(\rho) \ d\rho}}{\int{\rho \ n(\rho) \ d\rho}}$ diff --git a/process/main.py b/process/main.py index 4f2c7ccfa8..ed33770088 100644 --- a/process/main.py +++ b/process/main.py @@ -97,7 +97,7 @@ from process.models.physics.plasma_fields import PlasmaFields from process.models.physics.plasma_geometry import PlasmaGeom from process.models.physics.plasma_profiles import PlasmaProfile -from process.models.physics.profiles import NeProfile, TeProfile +from process.models.physics.profiles import ElectronDensityProfile, TeProfile from process.models.physics.scrape_off_layer import ScrapeOffLayer from process.models.power import Power from process.models.pulse import Pulse @@ -671,7 +671,7 @@ def __init__(self, data: DataStructure): self.pulse = Pulse() self.shield = Shield() self.ife = IFE(availability=self.availability, costs=self.costs) - self.ne_profile = NeProfile() + self.ne_profile = ElectronDensityProfile() self.te_profile = TeProfile() self.plasma_profile = PlasmaProfile(self.ne_profile, self.te_profile) self.fw = FirstWall() diff --git a/process/models/physics/plasma_profiles.py b/process/models/physics/plasma_profiles.py index d553319984..f18e7056f9 100644 --- a/process/models/physics/plasma_profiles.py +++ b/process/models/physics/plasma_profiles.py @@ -29,7 +29,8 @@ def __init__(self, ne_profile, te_profile): Args: profile_size (int): The size of the plasma profile. outfile (str): The output file path. - neprofile (NeProfile): An instance of the NeProfile class. + neprofile (ElectronDensityProfile): An instance of the ElectronDensityProfile + class. teprofile (TeProfile): An instance of the TeProfile class. """ # Default profile_size = 201, but it's possible to experiment with this value. @@ -195,7 +196,7 @@ def pedestal_parameterisation(self): `temp_plasma_electron_density_weighted_kev` and `temp_plasma_ion_density_weighted_kev`. """ - # Run TeProfile and NeProfile class methods: + # Run TeProfile and ElectronDensityProfile class methods: # Re-calculate core and profile values self.teprofile.run() diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 53342d9abc..69732e575f 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -139,7 +139,7 @@ def description(self): return self._description_ -class NeProfile(Profile): +class ElectronDensityProfile(Profile): """Electron density profile class. Contains a function to calculate the electron density profile and store the data. """ From 358dc02072a217b8fbb5dede4a4e8a47b5b4fc49 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:18:06 +0100 Subject: [PATCH 02/10] Refine ElectronDensityProfile class documentation and parameter names for clarity --- process/models/physics/profiles.py | 63 +++++++++++++++--------------- 1 file changed, 31 insertions(+), 32 deletions(-) diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 69732e575f..d54f4f4aa5 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -140,7 +140,7 @@ def description(self): class ElectronDensityProfile(Profile): - """Electron density profile class. Contains a function to calculate the electron + """Electron density (nₑ) profile class. Contains a function to calculate the electron density profile and store the data. """ @@ -151,12 +151,12 @@ def run(self): self.calculate_profile_dx() self.set_physics_variables() self.calculate_profile_y( - self.profile_x, - self.data.physics.radius_plasma_pedestal_density_norm, - self.data.physics.nd_plasma_electron_on_axis, - self.data.physics.nd_plasma_pedestal_electron, - self.data.physics.nd_plasma_separatrix_electron, - self.data.physics.alphan, + rho=self.profile_x, + radius_plasma_pedestal_density_norm=self.data.physics.radius_plasma_pedestal_density_norm, + nd_on_axis=self.data.physics.nd_plasma_electron_on_axis, + nd_pedestal=self.data.physics.nd_plasma_pedestal_electron, + nd_separatrix=self.data.physics.nd_plasma_separatrix_electron, + alphan=self.data.physics.alphan, ) self.integrate_profile_y() @@ -164,53 +164,52 @@ def calculate_profile_y( self, rho: np.array, radius_plasma_pedestal_density_norm: float, - n0: float, - nped: float, - nsep: float, + nd_on_axis: float, + nd_pedestal: float, + nd_separatrix: float, alphan: float, - ): - """Calculates the density at each normalised minor radius position - rho for a HELIOS-type density pedestal profile (neprofile). + ) -> None: + """Calculates the number density at each normalised minor radius (ρ) position. Parameters ---------- rho : - Normalised minor radius vector. + Normalised minor radius (ρ) vector. radius_plasma_pedestal_density_norm : - Normalised minor radius pedestal position. - n0 : - Central density (/m3). - nped : - Pedestal density (/m3). - nsep : - Separatrix density (/m3) + Normalised minor radius pedestal position (ρₙ,pedestal). + nd_on_axis : + Central number density (n₀) [/m³]. + nd_pedestal : + Pedestal density (n_pedestal) [/m³]. + nd_separatrix : + Separatrix density (n_sep) [/m³]. alphan : - Density peaking parameter. - """ + Density peaking parameter (αₙ). + """ # noqa: RUF002 if ( PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PARABOLIC_PROFILE ): - self.profile_y = n0 * (1 - rho**2) ** alphan + self.profile_y = nd_on_axis * (1 - rho**2) ** alphan # Input checks - if n0 < nped: + if nd_on_axis < nd_pedestal: logger.info( "NPROFILE: density pedestal is higher than core density. %s, %s", - nped, - n0, + nd_pedestal, + nd_on_axis, ) rho_index = rho <= radius_plasma_pedestal_density_norm self.profile_y[rho_index] = ( - nped - + (n0 - nped) + nd_pedestal + + (nd_on_axis - nd_pedestal) * (1 - (rho[rho_index] / radius_plasma_pedestal_density_norm) ** 2) ** alphan ) # Invert the rho_index - self.profile_y[~rho_index] = nsep + (nped - nsep) * (1 - rho[~rho_index]) / ( - 1 - radius_plasma_pedestal_density_norm - ) + self.profile_y[~rho_index] = nd_separatrix + (nd_pedestal - nd_separatrix) * ( + 1 - rho[~rho_index] + ) / (1 - radius_plasma_pedestal_density_norm) @staticmethod def ncore( From 39856bd4bcd542584a74f0bf1c833bed8699c3df Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:20:45 +0100 Subject: [PATCH 03/10] Rename ncore() to calculate_core_on_axis_density() for clarity and update references --- .../profiles/plasma_density_profile.md | 2 +- process/models/physics/profiles.py | 16 +++++++++------- .../unit/models/physics/test_plasma_profiles.py | 2 +- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/documentation/source/physics-models/profiles/plasma_density_profile.md b/documentation/source/physics-models/profiles/plasma_density_profile.md index 44516e6ca9..6f72c7c0b7 100644 --- a/documentation/source/physics-models/profiles/plasma_density_profile.md +++ b/documentation/source/physics-models/profiles/plasma_density_profile.md @@ -19,7 +19,7 @@ $$ n_{\text{i0}} = \left(\frac{n_\text{i}}{n_\text{e}}\right)n_{\text{e0}} $$ -#### Electron core density of a pedestalised profile | `ncore()` +#### Electron core density of a pedestalised profile | `calculate_core_on_axis_density()` This function calculates the core electron density for a pedestalsied profile (`i_plasma_pedestal == 1`). It takes in values of diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index d54f4f4aa5..feb73e2152 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -212,7 +212,7 @@ def calculate_profile_y( ) / (1 - radius_plasma_pedestal_density_norm) @staticmethod - def ncore( + def calculate_core_on_axis_density( radius_plasma_pedestal_density_norm: float, nped: float, nsep: float, @@ -345,12 +345,14 @@ def set_physics_variables(self): PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PEDESTAL_PROFILE ): - self.data.physics.nd_plasma_electron_on_axis = self.ncore( - self.data.physics.radius_plasma_pedestal_density_norm, - self.data.physics.nd_plasma_pedestal_electron, - self.data.physics.nd_plasma_separatrix_electron, - self.data.physics.nd_plasma_electrons_vol_avg, - self.data.physics.alphan, + self.data.physics.nd_plasma_electron_on_axis = ( + self.calculate_core_on_axis_density( + self.data.physics.radius_plasma_pedestal_density_norm, + self.data.physics.nd_plasma_pedestal_electron, + self.data.physics.nd_plasma_separatrix_electron, + self.data.physics.nd_plasma_electrons_vol_avg, + self.data.physics.alphan, + ) ) self.data.physics.nd_plasma_ions_on_axis = ( self.data.physics.nd_plasma_ions_total_vol_avg diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py index 2bcc5b8d59..70180da28e 100644 --- a/tests/unit/models/physics/test_plasma_profiles.py +++ b/tests/unit/models/physics/test_plasma_profiles.py @@ -74,7 +74,7 @@ def test_ncore(monkeypatch, plasmaprofile): nsep = 3.4294618459618943e19 nav = 7.4321e19 alphan = 1.0 - assert neprofile.ncore( + assert neprofile.calculate_core_on_axis_density( radius_plasma_pedestal_density_norm, nped, nsep, nav, alphan ) == pytest.approx(9.7756974320342041e19) From ce2b4005b4bf5e136a60b31c8a73fec461e8d802 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:33:27 +0100 Subject: [PATCH 04/10] Refactor ElectronDensityProfile: update parameter names and documentation for clarity --- process/models/physics/profiles.py | 80 ++++++++----------- .../models/physics/test_plasma_profiles.py | 6 +- 2 files changed, 37 insertions(+), 49 deletions(-) diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index feb73e2152..38a82e9322 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -178,11 +178,11 @@ def calculate_profile_y( radius_plasma_pedestal_density_norm : Normalised minor radius pedestal position (ρₙ,pedestal). nd_on_axis : - Central number density (n₀) [/m³]. + Central number density (n₀) [m⁻³]. nd_pedestal : - Pedestal density (n_pedestal) [/m³]. + Pedestal density (n_pedestal) [m⁻³]. nd_separatrix : - Separatrix density (n_sep) [/m³]. + Separatrix density (n_sep) [m⁻³]. alphan : Density peaking parameter (αₙ). """ # noqa: RUF002 @@ -214,58 +214,44 @@ def calculate_profile_y( @staticmethod def calculate_core_on_axis_density( radius_plasma_pedestal_density_norm: float, - nped: float, - nsep: float, - nav: float, + nd_pedestal: float, + nd_separatrix: float, + nd_vol_average: float, alphan: float, ) -> float: - """Calculates the core density of a pedestalised profile. - The solution comes from integrating and summing the two separate density profiles - for the core and pedestal region within their bounds. This has to be multiplied - by the torus volume element before integration which leads to an added rho term - in each part of the profile. When dividing by the volume of integration to get - the average density the simplification leads to a factor of 2 having to be - multiplied on to each of the integration results. This function for the average - density can then be re-arranged to calculate the central plasma density - n_0 / ncore. + """Calculates the core density (n₀) of a pedestalised profile. Parameters ---------- radius_plasma_pedestal_density_norm : - The normalised minor radius pedestal position. - nped : - The pedestal density (/m3). - nsep : - The separatrix density (/m3). - nav : - The electron density (/m3). - alphan : - The density peaking parameter + Normalised minor radius pedestal position (ρₙ,pedestal). + nd_pedestal: float, + The pedestal density (n_pedestal) [m⁻³]. + nd_separatrix: float, + The separatrix density (n_sep) [m⁻³]. + nd_vol_average: float, + The volume averaged density (⟨n⟩) [m⁻³]. + alphan: float, + The density peaking parameter (αₙ). Returns ------- : - The core density. - - References - ---------- - Jean, J. (2011). HELIOS: A Zero-Dimensional Tool for Next Step and Reactor - Studies. Fusion Science and Technology, 59(2), 308-349. - https://doi.org/10.13182/FST11-A11650 + The core on-axis density (n₀) [/m³]. """ - ncore = ( + nd_on_axis = ( 1 / (3 * radius_plasma_pedestal_density_norm**2) * ( - 3 * nav * (1 + alphan) - + nsep + 3 * nd_vol_average * (1 + alphan) + + nd_separatrix * (1 + alphan) * ( -2 + radius_plasma_pedestal_density_norm + radius_plasma_pedestal_density_norm**2 ) - - nped + - nd_pedestal * ( (1 + alphan) * (1 + radius_plasma_pedestal_density_norm) + (alphan - 2) * radius_plasma_pedestal_density_norm**2 @@ -273,16 +259,16 @@ def calculate_core_on_axis_density( ) ) - if ncore < 0.0: + if nd_on_axis < 0.0: # Allows solver to continue and # warns the user to raise the lower bound on nd_plasma_electrons_vol_avg # if the run did not converge logger.error( - "ncore is going negative when solving. Please raise the value of " - "nd_plasma_electrons_vol_avg and or its lower limit." + "nd_on_axis is going negative when solving. Please raise the value of " + "nd_plasma_electrons_vol_avg (⟨nₑ⟩) and or its lower limit." ) - ncore = 1.0e-6 - return ncore + nd_on_axis = 1.0e-6 + return nd_on_axis def set_pedestal_and_separatrix_values(self): """Sets the pedestal and separatrix density values based on the user input @@ -345,14 +331,12 @@ def set_physics_variables(self): PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PEDESTAL_PROFILE ): - self.data.physics.nd_plasma_electron_on_axis = ( - self.calculate_core_on_axis_density( - self.data.physics.radius_plasma_pedestal_density_norm, - self.data.physics.nd_plasma_pedestal_electron, - self.data.physics.nd_plasma_separatrix_electron, - self.data.physics.nd_plasma_electrons_vol_avg, - self.data.physics.alphan, - ) + self.data.physics.nd_plasma_electron_on_axis = self.calculate_core_on_axis_density( + radius_plasma_pedestal_density_norm=self.data.physics.radius_plasma_pedestal_density_norm, + nd_pedestal=self.data.physics.nd_plasma_pedestal_electron, + nd_separatrix=self.data.physics.nd_plasma_separatrix_electron, + nd_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, + alphan=self.data.physics.alphan, ) self.data.physics.nd_plasma_ions_on_axis = ( self.data.physics.nd_plasma_ions_total_vol_avg diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py index 70180da28e..98eae0eb08 100644 --- a/tests/unit/models/physics/test_plasma_profiles.py +++ b/tests/unit/models/physics/test_plasma_profiles.py @@ -75,7 +75,11 @@ def test_ncore(monkeypatch, plasmaprofile): nav = 7.4321e19 alphan = 1.0 assert neprofile.calculate_core_on_axis_density( - radius_plasma_pedestal_density_norm, nped, nsep, nav, alphan + radius_plasma_pedestal_density_norm=radius_plasma_pedestal_density_norm, + nd_pedestal=nped, + nd_separatrix=nsep, + nd_vol_avg=nav, + alphan=alphan, ) == pytest.approx(9.7756974320342041e19) From fbe62f072dd0102606c16d4e1af4622d6d0b2e0b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:36:03 +0100 Subject: [PATCH 05/10] Rename calculate_core_on_axis_density to calculate_pedestal_profile_on_axis_density and update references --- .../profiles/plasma_density_profile.md | 2 +- process/models/physics/profiles.py | 4 ++-- tests/unit/models/physics/test_plasma_profiles.py | 12 +++++++++--- 3 files changed, 12 insertions(+), 6 deletions(-) diff --git a/documentation/source/physics-models/profiles/plasma_density_profile.md b/documentation/source/physics-models/profiles/plasma_density_profile.md index 6f72c7c0b7..2c3875975f 100644 --- a/documentation/source/physics-models/profiles/plasma_density_profile.md +++ b/documentation/source/physics-models/profiles/plasma_density_profile.md @@ -19,7 +19,7 @@ $$ n_{\text{i0}} = \left(\frac{n_\text{i}}{n_\text{e}}\right)n_{\text{e0}} $$ -#### Electron core density of a pedestalised profile | `calculate_core_on_axis_density()` +#### Electron core density of a pedestalised profile | `calculate_pedestal_profile_on_axis_density()` This function calculates the core electron density for a pedestalsied profile (`i_plasma_pedestal == 1`). It takes in values of diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 38a82e9322..0c246719c4 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -212,7 +212,7 @@ def calculate_profile_y( ) / (1 - radius_plasma_pedestal_density_norm) @staticmethod - def calculate_core_on_axis_density( + def calculate_pedestal_profile_on_axis_density( radius_plasma_pedestal_density_norm: float, nd_pedestal: float, nd_separatrix: float, @@ -331,7 +331,7 @@ def set_physics_variables(self): PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PEDESTAL_PROFILE ): - self.data.physics.nd_plasma_electron_on_axis = self.calculate_core_on_axis_density( + self.data.physics.nd_plasma_electron_on_axis = self.calculate_pedestal_profile_on_axis_density( radius_plasma_pedestal_density_norm=self.data.physics.radius_plasma_pedestal_density_norm, nd_pedestal=self.data.physics.nd_plasma_pedestal_electron, nd_separatrix=self.data.physics.nd_plasma_separatrix_electron, diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py index 98eae0eb08..2bcc2faeb9 100644 --- a/tests/unit/models/physics/test_plasma_profiles.py +++ b/tests/unit/models/physics/test_plasma_profiles.py @@ -3,6 +3,8 @@ import numpy as np import pytest +from process.models.physics.profiles import PlasmaProfile + @pytest.fixture def plasmaprofile(process_models): @@ -51,7 +53,9 @@ class NeProfileParam(NamedTuple): ], ids=["baseline_2018"], ) -def test_neprofile(neprofileparam: ProfileParam, monkeypatch, plasmaprofile): +def test_neprofile( + neprofileparam: ProfileParam, monkeypatch, plasmaprofile: PlasmaProfile +): monkeypatch.setattr( plasmaprofile.data.physics, "n_plasma_profile_elements", @@ -62,7 +66,9 @@ def test_neprofile(neprofileparam: ProfileParam, monkeypatch, plasmaprofile): assert neprofile.profile_y == pytest.approx(neprofileparam.expected_neprofile) -def test_ncore(monkeypatch, plasmaprofile): +def test_calculate_pedestal_profile_on_axis_density( + monkeypatch, plasmaprofile: PlasmaProfile +): monkeypatch.setattr( plasmaprofile.data.physics, "n_plasma_profile_elements", @@ -74,7 +80,7 @@ def test_ncore(monkeypatch, plasmaprofile): nsep = 3.4294618459618943e19 nav = 7.4321e19 alphan = 1.0 - assert neprofile.calculate_core_on_axis_density( + assert neprofile.calculate_pedestal_profile_on_axis_density( radius_plasma_pedestal_density_norm=radius_plasma_pedestal_density_norm, nd_pedestal=nped, nd_separatrix=nsep, From c5ddc2480c8ed9c8ae7ec63801312051e11efb4a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:40:58 +0100 Subject: [PATCH 06/10] Rename TeProfile to ElectronTemperatureProfile and update references throughout the codebase --- .../source/physics-models/profiles/plasma_profiles.md | 8 ++++---- .../physics-models/profiles/plasma_temperature_profile.md | 2 +- process/main.py | 7 +++++-- process/models/physics/plasma_profiles.py | 4 ++-- process/models/physics/profiles.py | 2 +- 5 files changed, 13 insertions(+), 10 deletions(-) diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index 85bf35c275..c551e97507 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -4,7 +4,7 @@ In `PROCESS` the density, temperature and current profiles of the plasma for electrons and ions can take two forms depending on the switch value for `i_plasma_pedestal`. Either without a [pedestal](http://fusionwiki.ciemat.es/wiki/Pedestal), `i_plasma_pedestal == 0` or with a pedestal `i_plasma_pedestal == 1`. `i_plasma_pedestal == 0` is better suited for modelling L-mode plasmas, while `i_plasma_pedestal == 1` is better suited for modelling [H-mode](https://en.wikipedia.org/wiki/High-confinement_mode) plasmas. -The files responsible for calculating and storing the profiles are `plasma_profiles.py` and `profiles.py`. A central plasma profile object is created from the [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) class that contains attributes for the plasma density and temperature. The density and temperature profiles are in themselves objects of the [`Profile`](./plasma_profiles_abstract_class.md) abstract base class. [`Profile`](./plasma_profiles_abstract_class.md), [`ElectronDensityProfile`](plasma_density_profile.md) and [`TeProfile`](./plasma_temperature_profile.md) are all defined in `profiles.py`. [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) is exclusively in `plasma_profiles.py` +The files responsible for calculating and storing the profiles are `plasma_profiles.py` and `profiles.py`. A central plasma profile object is created from the [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) class that contains attributes for the plasma density and temperature. The density and temperature profiles are in themselves objects of the [`Profile`](./plasma_profiles_abstract_class.md) abstract base class. [`Profile`](./plasma_profiles_abstract_class.md), [`ElectronDensityProfile`](plasma_density_profile.md) and [`ElectronTemperatureProfile`](./plasma_temperature_profile.md) are all defined in `profiles.py`. [`PlasmaProfile`](plasma_profiles.md#plasma-profile-class-plasmaprofile) is exclusively in `plasma_profiles.py`
![UML of profiles](./uml_classes_PlasmaProfile.png){height="1000px"} @@ -238,7 +238,7 @@ The graph below is for a standard pedestal profile. You can vary its attributes ## Plasma Profile Class | `PlasmaProfile` ### Initialization | `__init__()` -The parent plasma profile class is `PlasmaProfile`. Initialization sets the profile class size and `neprofile` and `teprofile` to [`ElectronDensityProfile`](plasma_density_profile.md) & [`TeProfile`](plasma_temperature_profile.md) objects from `profiles.py` +The parent plasma profile class is `PlasmaProfile`. Initialization sets the profile class size and `neprofile` and `teprofile` to [`ElectronDensityProfile`](plasma_density_profile.md) & [`ElectronTemperatureProfile`](plasma_temperature_profile.md) objects from `profiles.py` ???+ Note @@ -266,7 +266,7 @@ Depending on the value of `i_plasma_pedestal` different functions will be ran, t If pedestal profile values are set they are reset to have values that agree with the original form of the parabolic profiles. Such that $\rho_{\text{ped}} = 1$ and that pedestal and separatrix densities and temperatures are zero. This will then warn the user in the terminal. -The density and temperature profile runner function [`TeProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is then called to re-calculate the profile and core values. +The density and temperature profile runner function [`ElectronTemperatureProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is then called to re-calculate the profile and core values. Ratio of density-weighted to volume-averaged temperature factor is calculated: @@ -593,7 +593,7 @@ $$ ##### `pedestal_parameterisation()` -The density and temperature profile runner function [`TeProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is firstly called to re-calculate the profile and core values. +The density and temperature profile runner function [`ElectronTemperatureProfile/ElectronDensityProfile.run()`](plasma_density_profile.md#runner-function-run) is firstly called to re-calculate the profile and core values. Perform integrations to calculate ratio of density-weighted to volume-averaged temperature, etc. Density-weighted temperature = $\frac{\int{nT \ dV}}{\int{n \ dV}}$, which is approximately equal to the ratio $\frac{\int{\rho \ n(\rho) T(\rho) \ d\rho}}{\int{\rho \ n(\rho) \ d\rho}}$ diff --git a/documentation/source/physics-models/profiles/plasma_temperature_profile.md b/documentation/source/physics-models/profiles/plasma_temperature_profile.md index dbea93848e..32b1b34a8d 100644 --- a/documentation/source/physics-models/profiles/plasma_temperature_profile.md +++ b/documentation/source/physics-models/profiles/plasma_temperature_profile.md @@ -1,4 +1,4 @@ -# Temperature Profile | `TeProfile(Profile)` +# Temperature Profile | `ElectronTemperatureProfile(Profile)` The temperature profile class is organised around a central runner function that is called each time the plasma is parameterised by the parent [`PlasmaProfile()`](plasma_profiles.md#plasma-profile-class-plasmaprofile) class. It is called by [`pedestal_parameterisation()`](plasma_profiles.md#pedestal_parameterisation) and [`parabolic parameterisation()`](plasma_profiles.md#parabolic_parameterisation). The sequence of the runner function can be seen below along with explanation of the following calculations. diff --git a/process/main.py b/process/main.py index ed33770088..412c092907 100644 --- a/process/main.py +++ b/process/main.py @@ -97,7 +97,10 @@ from process.models.physics.plasma_fields import PlasmaFields from process.models.physics.plasma_geometry import PlasmaGeom from process.models.physics.plasma_profiles import PlasmaProfile -from process.models.physics.profiles import ElectronDensityProfile, TeProfile +from process.models.physics.profiles import ( + ElectronDensityProfile, + ElectronTemperatureProfile, +) from process.models.physics.scrape_off_layer import ScrapeOffLayer from process.models.power import Power from process.models.pulse import Pulse @@ -672,7 +675,7 @@ def __init__(self, data: DataStructure): self.shield = Shield() self.ife = IFE(availability=self.availability, costs=self.costs) self.ne_profile = ElectronDensityProfile() - self.te_profile = TeProfile() + self.te_profile = ElectronTemperatureProfile() self.plasma_profile = PlasmaProfile(self.ne_profile, self.te_profile) self.fw = FirstWall() self.blanket_library = BlanketLibrary(fw=self.fw) diff --git a/process/models/physics/plasma_profiles.py b/process/models/physics/plasma_profiles.py index f18e7056f9..ab638a8961 100644 --- a/process/models/physics/plasma_profiles.py +++ b/process/models/physics/plasma_profiles.py @@ -31,7 +31,7 @@ def __init__(self, ne_profile, te_profile): outfile (str): The output file path. neprofile (ElectronDensityProfile): An instance of the ElectronDensityProfile class. - teprofile (TeProfile): An instance of the TeProfile class. + teprofile (ElectronTemperatureProfile): An instance of the ElectronTemperatureProfile class. """ # Default profile_size = 201, but it's possible to experiment with this value. # See `n_plasma_profile_elements` @@ -196,7 +196,7 @@ def pedestal_parameterisation(self): `temp_plasma_electron_density_weighted_kev` and `temp_plasma_ion_density_weighted_kev`. """ - # Run TeProfile and ElectronDensityProfile class methods: + # Run ElectronTemperatureProfile and ElectronDensityProfile class methods: # Re-calculate core and profile values self.teprofile.run() diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 0c246719c4..1a1937d116 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -345,7 +345,7 @@ def set_physics_variables(self): ) -class TeProfile(Profile): +class ElectronTemperatureProfile(Profile): """Electron temperature profile class. Contains a function to calculate the temperature profile and store the data. """ From 1d68fa74214c0af86c5a5d25c1d3612ea9f8026c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 15:52:49 +0100 Subject: [PATCH 07/10] Refactor ElectronTemperatureProfile: enhance parameter names and documentation for clarity --- process/models/physics/profiles.py | 81 ++++++++++++++---------------- 1 file changed, 37 insertions(+), 44 deletions(-) diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 1a1937d116..42053c04d5 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -196,7 +196,7 @@ def calculate_profile_y( if nd_on_axis < nd_pedestal: logger.info( - "NPROFILE: density pedestal is higher than core density. %s, %s", + "NPROFILE: Pedestal density is higher than core density. %s, %s", nd_pedestal, nd_on_axis, ) @@ -346,7 +346,7 @@ def set_physics_variables(self): class ElectronTemperatureProfile(Profile): - """Electron temperature profile class. Contains a function to calculate the + """Electron temperature (Tₑ) profile class. Contains a function to calculate the temperature profile and store the data. """ @@ -357,13 +357,13 @@ def run(self): self.calculate_profile_dx() self.set_physics_variables() self.calculate_profile_y( - self.profile_x, - self.data.physics.radius_plasma_pedestal_temp_norm, - self.data.physics.temp_plasma_electron_on_axis_kev, - self.data.physics.temp_plasma_pedestal_kev, - self.data.physics.temp_plasma_separatrix_kev, - self.data.physics.alphat, - self.data.physics.tbeta, + rho=self.profile_x, + radius_plasma_pedestal_temp_norm=self.data.physics.radius_plasma_pedestal_temp_norm, + temp_on_axis_kev=self.data.physics.temp_plasma_electron_on_axis_kev, + temp_pedestal_kev=self.data.physics.temp_plasma_pedestal_kev, + temp_separatrix_kev=self.data.physics.temp_plasma_separatrix_kev, + alphat=self.data.physics.alphat, + tbeta=self.data.physics.tbeta, ) self.integrate_profile_y() @@ -371,31 +371,29 @@ def calculate_profile_y( self, rho: np.array, radius_plasma_pedestal_temp_norm: float, - t0: float, - temp_plasma_pedestal_kev: float, - temp_plasma_separatrix_kev: float, + temp_on_axis_kev: float, + temp_pedestal_kev: float, + temp_separatrix_kev: float, alphat: float, tbeta: float, - ): - """Calculates the temperature at a normalised minor radius position rho for a - pedestalised profile (teprofile). - If i_plasma_pedestal = 0 the original parabolic profile form is used instead. + ) -> None: + """Calculates the temperature at each normalised minor radius (ρ) position. Parameters ---------- - rho : np.array - Normalised minor radius. - radius_plasma_pedestal_temp_norm : float - Normalised minor radius pedestal position. - t0 : float - Central temperature (keV). - temp_plasma_pedestal_kev : float - Pedestal temperature (keV). - temp_plasma_separatrix_kev : float - Separatrix temperature (keV). - alphat : float - Temperature peaking parameter. - tbeta : float + rho : + Normalised minor radius (ρ) vector + radius_plasma_pedestal_temp_norm : + Normalised minor radius pedestal position (ρₜ,pedestal). + temp_on_axis_kev : + Central on-axis temperature (T₀) [keV]. + temp_pedestal_kev : + Pedestal temperature (T_pedestal) [keV]. + temp_separatrix_kev : + Separatrix temperature (T_separatrix) [keV]. + alphat : + Temperature peaking parameter (αₜ). + tbeta : Second temperature exponent. Raises @@ -403,12 +401,7 @@ def calculate_profile_y( ProcessValueError If negative temperature in plasma profile - References - ---------- - Jean, J. (2011). HELIOS: A Zero-Dimensional Tool for Next Step and Reactor - Studies. Fusion Science and Technology, 59(2), 308-349. - https://doi.org/10.13182/FST11-A11650 - """ + """ # noqa: RUF002 if ( PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PARABOLIC_PROFILE @@ -416,25 +409,25 @@ def calculate_profile_y( # profile values of 0 cause divide by 0 errors so ensure the profile value # is at least 1e-8 # which is small enough that it won't make a difference to any calculations - self.profile_y = np.maximum(t0 * (1 - rho**2) ** alphat, 1e-8) + self.profile_y = np.maximum(temp_on_axis_kev * (1 - rho**2) ** alphat, 1e-8) return - if t0 < temp_plasma_pedestal_kev: + if temp_on_axis_kev < temp_pedestal_kev: logger.info( - "TPROFILE: temperature pedestal is higher than core temperature. %s, %s", - temp_plasma_pedestal_kev, - t0, + "TPROFILE: Pedestal temperature is higher than core temperature. %s, %s", + temp_pedestal_kev, + temp_on_axis_kev, ) rho_index = rho <= radius_plasma_pedestal_temp_norm self.profile_y[rho_index] = ( - temp_plasma_pedestal_kev - + (t0 - temp_plasma_pedestal_kev) + temp_pedestal_kev + + (temp_on_axis_kev - temp_pedestal_kev) * (1 - (rho[rho_index] / radius_plasma_pedestal_temp_norm) ** tbeta) ** alphat ) - self.profile_y[~rho_index] = temp_plasma_separatrix_kev + ( - temp_plasma_pedestal_kev - temp_plasma_separatrix_kev + self.profile_y[~rho_index] = temp_separatrix_kev + ( + temp_pedestal_kev - temp_separatrix_kev ) * (1 - rho[~rho_index]) / (1 - radius_plasma_pedestal_temp_norm) # Check for any negative temperature in profile: always fatal in From e8bdee6552b379e3fb9fa5e59587475540d3315e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 14 Aug 2026 16:00:31 +0100 Subject: [PATCH 08/10] Rename tcore to calculate_pedestal_profile_on_axis_temperature and update references in the ElectronTemperatureProfile class and tests --- .../profiles/plasma_temperature_profile.md | 2 +- process/models/physics/profiles.py | 77 ++++++++----------- .../models/physics/test_plasma_profiles.py | 14 +++- 3 files changed, 45 insertions(+), 48 deletions(-) diff --git a/documentation/source/physics-models/profiles/plasma_temperature_profile.md b/documentation/source/physics-models/profiles/plasma_temperature_profile.md index 32b1b34a8d..00b129e0b1 100644 --- a/documentation/source/physics-models/profiles/plasma_temperature_profile.md +++ b/documentation/source/physics-models/profiles/plasma_temperature_profile.md @@ -14,7 +14,7 @@ The temperature profile class is organised around a central runner function that The core electron temperature is calculated using the [`tcore`](plasma_temperature_profile.md#electron-core-density-of-a-pedestalised-profile--tcore) method. -#### Electron core density of a pedestalised profile | `tcore()` +#### Electron core density of a pedestalised profile | `calculate_pedestal_profile_on_axis_temperature()` This function calculates the core electron density for a pedestalsied profile in $\text{keV}$. The inclusion of a new $\beta_T$ exponent term allows a more accurate description of temperature profiles with a triangular shape or a strong gradient near the pedestal (characteristic of regimes with an [internal transport barrier](https://wiki.fusion.ciemat.es/wiki/Internal_Transport_Barrier)). diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 42053c04d5..76989ae4ab 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -394,7 +394,7 @@ def calculate_profile_y( alphat : Temperature peaking parameter (αₜ). tbeta : - Second temperature exponent. + Second temperature exponent (βₜ). Raises ------ @@ -436,64 +436,51 @@ def calculate_profile_y( raise ProcessValueError("Negative temperature in plasma profile") @staticmethod - def tcore( + def calculate_pedestal_profile_on_axis_temperature( radius_plasma_pedestal_temp_norm: float, - temp_plasma_pedestal_kev: float, - temp_plasma_separatrix_kev: float, - tav: float, + temp_pedestal_kev: float, + temp_separatrix_kev: float, + temp_vol_avg_kev: float, alphat: float, tbeta: float, ) -> float: - """Calculates the core temperature (keV) - of a pedestalised profile. The solution comes from integrating and summing the - two separate temperature profiles for the core and pedestal region within their - bounds. This has to be multiplied by the torus volume element before integration - which leads to an added rho term in each part of the profile. When dividing by - the volume of integration to get the average temperature the simplification - leads to a factor of 2 having to be multiplied on to each of the integration - results. This function for the average temperature can then be re-arranged to - calculate the central plasma temperature T_0 / tcore. + """Calculates the core density (T₀) of a pedestalised profile. Parameters ---------- - radius_plasma_pedestal_temp_norm : float - Normalised minor radius pedestal position. - temp_plasma_pedestal_kev : float - Pedestal temperature (keV). - temp_plasma_separatrix_kev : float - Separatrix temperature (keV). - tav : float - Volume average temperature (keV). - alphat : float - Temperature peaking parameter. - tbeta : float - Second temperature exponent. + radius_plasma_pedestal_temp_norm : + Normalised minor radius pedestal position (ρₜ,pedestal). + temp_pedestal_kev : + Pedestal temperature (T_pedestal) [keV]. + temp_separatrix_kev : + Separatrix temperature (T_separatrix) [keV]. + temp_vol_avg_kev : + Volume average temperature (⟨T⟩) [keV]. + alphat : + Temperature peaking parameter (αₜ). + tbeta : + Second temperature exponent (βₜ). Returns ------- - float - Core temperature. + : + The core on-axis temperature (T₀) [keV] - References - ---------- - Jean, J. (2011). HELIOS: A Zero-Dimensional Tool for Next Step and Reactor - Studies. Fusion Science and Technology, 59(2), 308-349. - https://doi.org/10.13182/FST11-A11650 """ # Calculate core temperature - return temp_plasma_pedestal_kev + ( + return temp_pedestal_kev + ( ( tbeta * ( - 3 * tav - + temp_plasma_separatrix_kev + 3 * temp_vol_avg_kev + + temp_separatrix_kev * ( -2.0 + radius_plasma_pedestal_temp_norm + radius_plasma_pedestal_temp_norm**2 ) - - temp_plasma_pedestal_kev + - temp_pedestal_kev * ( 1 + radius_plasma_pedestal_temp_norm @@ -522,13 +509,15 @@ def set_physics_variables(self): PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PEDESTAL_PROFILE ): - self.data.physics.temp_plasma_electron_on_axis_kev = self.tcore( - self.data.physics.radius_plasma_pedestal_temp_norm, - self.data.physics.temp_plasma_pedestal_kev, - self.data.physics.temp_plasma_separatrix_kev, - self.data.physics.temp_plasma_electron_vol_avg_kev, - self.data.physics.alphat, - self.data.physics.tbeta, + self.data.physics.temp_plasma_electron_on_axis_kev = ( + self.calculate_pedestal_profile_on_axis_temperature( + self.data.physics.radius_plasma_pedestal_temp_norm, + self.data.physics.temp_plasma_pedestal_kev, + self.data.physics.temp_plasma_separatrix_kev, + self.data.physics.temp_plasma_electron_vol_avg_kev, + self.data.physics.alphat, + self.data.physics.tbeta, + ) ) self.data.physics.temp_plasma_ion_on_axis_kev = ( diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py index 2bcc2faeb9..d2531621ab 100644 --- a/tests/unit/models/physics/test_plasma_profiles.py +++ b/tests/unit/models/physics/test_plasma_profiles.py @@ -139,7 +139,10 @@ def test_teprofile(teprofileparam: ProfileParam, monkeypatch, plasmaprofile): assert teprofile.profile_y == pytest.approx(teprofileparam.expected_teprofile) -def test_tcore(plasmaprofile): +def test_calculate_pedestal_profile_on_axis_temperature(plasmaprofile: PlasmaProfile): + """ + Test the calculation of the pedestal profile on-axis temperature. + """ teprofile = plasmaprofile.teprofile radius_plasma_pedestal_temp_norm = 0.94 tped = 3.7775374842470044 @@ -148,8 +151,13 @@ def test_tcore(plasmaprofile): alphat = 1.45 tbeta = 2.0 - assert teprofile.tcore( - radius_plasma_pedestal_temp_norm, tped, tsep, tav, alphat, tbeta + assert teprofile.calculate_pedestal_profile_on_axis_temperature( + radius_plasma_pedestal_temp_norm=radius_plasma_pedestal_temp_norm, + temp_pedestal_kev=tped, + temp_separatrix_kev=tsep, + temp_vol_avg_kev=tav, + alphat=alphat, + tbeta=tbeta, ) == pytest.approx(28.09093632260765) From 475ef4c4b39a1a6f706224d3caa0c5090dd96dc9 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 24 Aug 2026 15:47:15 +0100 Subject: [PATCH 09/10] Add function to calculate volume average of a radially normalized profile --- process/models/physics/profiles.py | 38 ++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 76989ae4ab..b0c9c1f9b3 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -525,3 +525,41 @@ def set_physics_variables(self): / self.data.physics.temp_plasma_electron_vol_avg_kev * self.data.physics.temp_plasma_electron_on_axis_kev ) + + +def calculate_vol_avg_of_profile(profile_x: np.ndarray, profile_y: np.ndarray) -> float: + """Calculate the volume averaged value (⟨profile_y⟩) of a radially normalised profile. + + Parameters + ---------- + profile_x : np.ndarray + The x-values of the profile. + profile_y : np.ndarray + The y-values of the profile. + + Notes + ----- + - The 2 factor in the calculation arises using both sides of the profile and the + radial normalisation of the profile. + + Returns + ------- + float: The volume-averaged value (⟨profile_y⟩) of the profile. + """ + if profile_x.ndim != 1: + raise ValueError("profile_x must be a 1D array.") + + if profile_x.size < 2: + raise ValueError("profile_x must contain at least 2 points.") + + if not np.isclose(profile_x[0], 0.0) or not np.isclose(profile_x[-1], 1.0): + raise ValueError("profile_x must span from 0 to 1.") + + if np.any(np.diff(profile_x) <= 0): + raise ValueError("profile_x must be strictly increasing.") + + return 2.0 * sp.integrate.simpson( + profile_y * profile_x, + x=profile_x, + dx=profile_x[1] - profile_x[0], + ) From 73dbd766029795d4e4ceb4cd505198d7c6e90019 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 24 Aug 2026 16:09:21 +0100 Subject: [PATCH 10/10] Add documentation for volume average calculation and update variable name for clarity --- .../profiles/plasma_profiles.md | 41 +++++++++++++++++++ process/models/physics/profiles.py | 2 +- 2 files changed, 42 insertions(+), 1 deletion(-) diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index c551e97507..39498e0717 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -636,6 +636,47 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found ----- + + +## Calculate profile volume average | `calculate_vol_avg_of_profile()` + +General method to calculate the volume averaged value of any normalised radial profile. + +### Derivation + +The integral $I$ of any quantity $G(\rho)$ over the volume of the plasma where $\rho$ is the normalised minor radius is: + +$$ +I=\int_{}^{} G(\rho) \ dV +$$ + +The volume element is the surface area of the toroidal surface times the thickness $dV=2\pi R_0 2\pi\kappa r dr $ + +$$ +I=\int_{0}^{1} G(\rho) 2\pi R_0 2\pi\kappa r \ dr +$$ + +Substituting the normalised minor radius coordinate $\rho=r/a$, and bringing the constants out in front, + +$$ +I=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1} G(\rho) \rho \ d\rho +$$ + +However, if we calculate the volume $V$, the exact same integral is used where $G=1$, + +$$ +V=4 \pi^2 R_0 \kappa a^2 \int_{0}^{1} \rho \ d\rho = \frac{1}{2} 4 \pi^2 R_0 \kappa a^2 +$$ + +Therefore the general volume integral of $G(\rho)$, where $V$ is the volume, is + +$$ +I=2 V \int_{0}^{1} G(\rho) \rho \ d\rho +$$ + +If $G$ is a density then $2 \int_{0}^{1} G(\rho) \rho \ d\rho$ is the volume averaged density. + + ## Key Constraints -------- diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index b0c9c1f9b3..e320e30caf 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -335,7 +335,7 @@ def set_physics_variables(self): radius_plasma_pedestal_density_norm=self.data.physics.radius_plasma_pedestal_density_norm, nd_pedestal=self.data.physics.nd_plasma_pedestal_electron, nd_separatrix=self.data.physics.nd_plasma_separatrix_electron, - nd_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, + nd_vol_average=self.data.physics.nd_plasma_electrons_vol_avg, alphan=self.data.physics.alphan, ) self.data.physics.nd_plasma_ions_on_axis = (