🐛 Fix error in volume averaged thermal plasma pressure to be prof…

documentation/physics-models/profiles/plasma_profiles.md

@@ -271,7 +271,8 @@ The density and temperature profile runner function [`TeProfile/NeProfile.run()`
 Ratio of density-weighted to volume-averaged temperature factor is calculated:
 
 $$
-\texttt{pcoef} = \frac{(1+\alpha_n)(1+\alpha_T)}{(1+\alpha_T + \alpha_n)}
+\texttt{f_temp_plasma_electron_density_vol_avg} =\frac{\langle T_{\text{e}} \rangle_{\text{n}}}{\underbrace{\langle T_{\text{e}} \rangle_{\text{V}}}_{\texttt{temp_plasma_electrons_vol_avg}}}\\
+ = \frac{(1+\alpha_n)(1+\alpha_T)}{(1+\alpha_T + \alpha_n)}
 $$
 
 --------
@@ -310,8 +311,10 @@ $\blacksquare$
 The density weighted temperatures are set:
 
 $$\begin{aligned}
-\texttt{temp_plasma_electron_density_weighted_kev} = \texttt{pcoef} \times T_\text{e} \\
-\texttt{temp_plasma_ion_density_weighted_kev} = \texttt{pcoef}\times T_\text{i}
+\overbrace{\langle T_{\text{e}} \rangle_{\text{n}}}^{\texttt{temp_plasma_electron_density_weighted_kev}} = \\
+\texttt{f_temp_plasma_electron_density_vol_avg} \times \langle T_\text{e} \rangle  \\
+\overbrace{\langle T_{\text{i}} \rangle_{\text{n}}}^{\texttt{temp_plasma_ion_density_weighted_kev}} = \\
+\texttt{f_temp_plasma_electron_density_vol_avg}\times \langle T_\text{i} \rangle 
 \end{aligned}$$
 
 ------
@@ -392,34 +395,43 @@ $$\begin{aligned}
 
 ##### `calculate_profile_factors()`
 
-The central plasma pressure is calculated from the ideal gas law.
+The central thermal plasma pressure is calculated from the ideal gas law.
 
 $$
-p_0 = (\texttt{nd_plasma_electron_on_axis} \times \texttt{temp_plasma_electron_on_axis_kev}\\
-+\texttt{nd_plasma_ions_on_axis}\times \texttt{temp_plasma_ion_on_axis_kev})\times (1000 \times \text{e})
+p_{\text{thermal,0}} = (\overbrace{n_{\text{e,0}}}^{\texttt{nd_plasma_electron_on_axis}} \times \overbrace{T_{\text{e,0}}}^{\texttt{temp_plasma_electron_on_axis_kev}}\\
++\underbrace{n_{n\text{i,0}}}_{\texttt{nd_plasma_ions_on_axis}}\times \underbrace{T_{\text{i,0}}}_{\texttt{temp_plasma_ion_on_axis_kev})}\times (\times \text{keV})
 $$
 
 With the coefficients used to turn the temperature from $\text{keV}$ back to Joules.
 
-Since we multiply the density and temperature profiles together to get the pressure we can find the pressure profile factor by adding the profile factors for temperature and density.
+
+We calculate the volume averaged thermal plasma pressure using the volume averaged density and the density weighted volume averaged temperature for the electrons and ions. The density weighted temperature is used as $\langle nT \rangle \neq \langle n \rangle \langle T \rangle$. This can be seen also in the [equality constraint for thermal $\beta$](../plasma_beta/plasma_beta.md#beta-consistency)
 
 $$
-\alpha_p = \alpha_n + \alpha_T
+\langle p_{\text{thermal}} \rangle = \overbrace{\langle n_{\text{e}} \rangle_{\text{V}}}^{\texttt{nd_plasma_electrons_vol_avg}} \times \overbrace{\langle T_{\text{e}} \rangle_{\text{n}}}^{\texttt{temp_plasma_electron_density_weighted_kev}}\\
+ + \underbrace{\langle n_{\text{i}} \rangle_{\text{V}}}_{\texttt{nd_plasma_ions_total_vol_avg}} \times \underbrace{\langle T_{\text{i}} \rangle_{\text{n}}}_{\texttt{temp_plasma_ion_density_weighted_kev}}
 $$
 
-The volume averaged pressure can then be set if we assume the pressure also has a parabolic profile. Using the standard relation used for both density and temeprature we can set the volume averaged pressure as:
+???+ note "Parabolic pressure profile factor"
 
-$$
-\langle p \rangle =   \frac{p_0}{\alpha_p+1}
-$$
+    For a prabolic profile, since we multiply the density and temperature profiles together to get the pressure we can find the pressure profile factor by adding the profile factors for temperature and density.
 
-???+ note "Pressure profile factor"
+    $$
+    \alpha_p = \alpha_n + \alpha_T
+    $$
+
+    The volume averaged pressure can then be set if we assume the pressure also has a parabolic profile. Using the standard relation used for both density and temeprature we can set the volume averaged pressure as:
+
+    $$
+    \langle p \rangle =   \frac{p_0}{\alpha_p+1}
+    $$
 
     The calculation of $\alpha_p$ is only valid assuming a parabolic profile case. The calculatio of $p_0$ is still true as the core values are calculated independantly for each profile type.
 
     $p_0$ is NOT equal to $\langle p \rangle \times (1 + \alpha_p)$, but $p(\rho) = n(\rho)T(\rho)$ and $\langle p \rangle = \langle n \rangle$ $T_n$, where $T_n$ is the
     density-weighted temperature.
 
+
 The on-axis plasma current density ($J_0$) is then calculated by assuming a parabolic profile and using the calculated plasma current ($I_{\text{plasma}}$) and the plasmas poloidal cross-sectional area ($A_{\text{pol}}$).
 
 $$
@@ -588,24 +600,22 @@ The density and temperature profile runner function [`TeProfile/NeProfile.run()`
  Density weighted temperatures are thus set as such
 
 $$
-\texttt{temp_plasma_electron_density_weighted_kev} = \\
-\frac{\int_0^1{\rho \ n(\rho) T(\rho) \  d\rho}}{\int_0^1{\rho \ n(\rho) \ d\rho}} \\
-\texttt{temp_plasma_ion_density_weighted_kev} = \\
-\texttt{temp_plasma_electron_density_weighted_kev}\times \frac{T_\text{i}}{T_\text{e}}
+ \overbrace{\langle T_{\text{e}} \rangle_{\text{n}}}^{\texttt{temp_plasma_electron_density_weighted_kev}} = \frac{\int_0^1{\rho \ n_{\text{e}}(\rho) T_{\text{e}}(\rho) \  d\rho}}{\int_0^1{\rho \ n_{\text{e}}(\rho) \ d\rho}} \\
+\overbrace{\langle T_{\text{i}} \rangle_{\text{n}}}^{\texttt{temp_plasma_ion_density_weighted_kev}} = \langle T_{\text{e}} \rangle_{\text{n}}\times \frac{\langle T_{\text{i}} \rangle_{\text{V}}}{\langle T_{\text{e}} \rangle_{\text{V}}}
 $$
 
 The above is done numerically with [scipy.integrate.simpson](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.simpson.html) integration.
 
 Set profile factor which is the ratio of density-weighted to volume-averaged temperature
 
 $$
-\texttt{pcoef} = \frac{\texttt{temp_plasma_electron_density_weighted_kev}}{T_\text{e}}
+\texttt{f_temp_plasma_electron_density_vol_avg} =\frac{\langle T_{\text{e}} \rangle_{\text{n}}}{\underbrace{\langle T_{\text{e}} \rangle_{\text{V}}}_{\texttt{temp_plasma_electrons_vol_avg}}}
 $$
 
 Calculate the line averaged electron density by integrating the normalised profile using the class [`integrate_profile_y()`](./plasma_profiles_abstract_class.md#calculate-the-profile-integral-value-integrate_profile_y) function
 
 $$
-\texttt{nd_plasma_electron_line} = \int_0^1{n(\rho) \ d\rho}
+ \overbrace{\bar{n_{\text{e}}}}^{\texttt{nd_plasma_electron_line}} = \int_0^1{n(\rho) \ d\rho}
 $$
 
 A divertor variable `prn1` is set to be equal to the separatrix density over the mean density:

process/io/plot_proc.py

@@ -10567,6 +10567,18 @@ def plot_plasma_pressure_profiles(axis, mfile_data, scan):
         color="green",
         label="Total",
     )
+
+    # Plot horizontal line for volume-average thermal pressure (converted to kPa)
+    p_vol_kpa = mfile_data.data["pres_plasma_thermal_vol_avg"].get_scan(scan) / 1000.0
+    axis.axhline(
+        p_vol_kpa,
+        color="black",
+        linestyle="--",
+        linewidth=1.2,
+        label="Volume avg",
+        zorder=5,
+    )
+
     axis.set_xlabel("$\\rho$ [r/a]")
     axis.set_ylabel("Thermal Pressure [kPa]")
     axis.minorticks_on()
@@ -10577,6 +10589,22 @@ def plot_plasma_pressure_profiles(axis, mfile_data, scan):
     axis.set_ylim(bottom=0)
     axis.legend()
 
+    textstr_pressure = "\n".join((
+        rf"$p_0$: {mfile_data.data['pres_plasma_thermal_on_axis'].get_scan(scan) / 1000:,.3f} kPa",
+        rf"$\langle p_{{\text{{total}}}} \rangle_\text{{V}}$: {mfile_data.data['pres_plasma_thermal_vol_avg'].get_scan(scan) / 1000:,.3f} kPa",
+    ))
+
+    axis.text(
+        0.5,
+        1.2,
+        textstr_pressure,
+        transform=axis.transAxes,
+        fontsize=9,
+        verticalalignment="top",
+        horizontalalignment="center",
+        bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5},
+    )
+
 
 def plot_plasma_pressure_gradient_profiles(axis, mfile_data, scan):
     # Get the plasma pressure profiles

process/plasma_profiles.py

@@ -273,15 +273,11 @@ def calculate_profile_factors(self) -> None:
         #  Central pressure (Pa), from ideal gas law : p = nkT
 
         physics_variables.pres_plasma_thermal_on_axis = (
-            (
-                physics_variables.nd_plasma_electron_on_axis
-                * physics_variables.temp_plasma_electron_on_axis_kev
-                + physics_variables.nd_plasma_ions_on_axis
-                * physics_variables.temp_plasma_ion_on_axis_kev
-            )
-            * 1.0e3
-            * constants.ELECTRON_CHARGE
-        )
+            physics_variables.nd_plasma_electron_on_axis
+            * physics_variables.temp_plasma_electron_on_axis_kev
+            + physics_variables.nd_plasma_ions_on_axis
+            * physics_variables.temp_plasma_ion_on_axis_kev
+        ) * constants.KILOELECTRON_VOLT
 
         # Electron pressure profile (Pa)
         physics_variables.pres_plasma_electron_profile = self.neprofile.profile_y * (
@@ -314,19 +310,21 @@ def calculate_profile_factors(self) -> None:
             * physics_variables.f_temp_plasma_ion_electron
         )
 
-        #  Pressure profile index (N.B. no pedestal effects included here)
+        #  Pressure profile index (only true for a parabolic profile)
         #  N.B. pres_plasma_thermal_on_axis is NOT equal to <p> * (1 + alphap), but p(rho) = n(rho)*T(rho)
         #  and <p> = <n>.T_n where <...> denotes volume-averages and T_n is the
         #  density-weighted temperature
 
         physics_variables.alphap = physics_variables.alphan + physics_variables.alphat
 
-        # Shall assume that the pressure profile is parabolic. Can find volume average from
-        # profile index and core value the same as for density and temperature
+        # Calculate the volume averaged plasma thermal pressure from the density-weighted temperatures
+        # Density-weighted temperatures are used as <nT> != <n>*<T>
         physics_variables.pres_plasma_thermal_vol_avg = (
-            physics_variables.pres_plasma_thermal_on_axis
-            / (physics_variables.alphap + 1)
-        )
+            physics_variables.nd_plasma_electrons_vol_avg
+            * physics_variables.temp_plasma_electron_density_weighted_kev
+            + physics_variables.nd_plasma_ions_total_vol_avg
+            * physics_variables.temp_plasma_ion_density_weighted_kev
+        ) * constants.KILOELECTRON_VOLT
 
         # Central plasma current density (A/m^2)
         # Assumes a parabolic profile for the current density