diff --git a/documentation/source/physics-models/profiles/plasma_density_profile.md b/documentation/source/physics-models/profiles/plasma_density_profile.md index b65a30674c..87742457d9 100644 --- a/documentation/source/physics-models/profiles/plasma_density_profile.md +++ b/documentation/source/physics-models/profiles/plasma_density_profile.md @@ -148,4 +148,31 @@ $$\begin{aligned} 5. Profile is then integrated with `integrate_profile_y()` using Simpsons integration from the profile abstract base class +----------------------- + +### Setting pedestal and separatrix values | `set_pedestal_and_separatrix_values()` + +The switch `i_nd_plasma_pedestal_separatrix` controls how the values of the density pedestal and separatrix are set. + +#### User input + +If `i_nd_plasma_pedestal_separatrix == 0` then the values of `nd_plasma_pedestal_electron` and `nd_plasma_separatrix_electron` are taken directly from the input file + +#### Fraction of Greenwald Limit + +If `i_nd_plasma_pedestal_separatrix == 1`, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as: + +$$ +n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14} +$$ + +$$ +n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14} +$$ + + +`f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald` can be set as iteration variables respectively by using `ixc = 145` +and `ixc = 152` respectively + + [^1]: Jean, J. (2011). *HELIOS: A Zero-Dimensional Tool for Next Step and Reactor Studies*. Fusion Science and Technology, 59(2), 308–349. diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index c6514b5d51..985d975768 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -277,10 +277,10 @@ $$ -------- -The line averaged density is then calculated for the profile paramaters +The line averaged density/temperature is then calculated for the profile paramaters -###### Line averaged density derivation -Line averaged electron density is calculated by integrating the profile across the normalised width of the profile and then dividing by the width of the integration bounds +###### Line averaged density/temperature derivation +Line averaged electron density/temperature is calculated by integrating the profile across the normalised width of the profile and then dividing by the width of the integration bounds. Using the density as an example: $$ \overbrace{\bar{n_{\text{e}}}}^{\texttt{nd_plasma_electron_line}} = \frac{\int^1_0 n_0(1-\rho^2)^{\alpha_n} \ d\rho}{\rho} @@ -612,12 +612,16 @@ $$ \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 +Calculate the line averaged electron density and temperature 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 $$ \overbrace{\bar{n_{\text{e}}}}^{\texttt{nd_plasma_electron_line}} = \int_0^1{n(\rho) \ d\rho} $$ +$$ + \overbrace{\bar{T_{\text{e}}}}^{\texttt{temp_plasma_electron_line_avg_kev}} = \int_0^1{T(\rho) \ d\rho} +$$ + A divertor variable `prn1` is set to be equal to the separatrix density over the mean density: $$ @@ -636,24 +640,6 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found -------- -### Setting pedestal values as fractions of the Greenwald limit - -By default, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as: - -$$ -n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14} -$$ - -$$ -n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14} -$$ - -To set the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ directly, the user can input the value of $\texttt{f_nd_plasma_pedestal_greenwald}$ or $\texttt{f_nd_plasma_separatrix_greenwald}$ to be less than 0.0 (i.e negative) to prevent the Greenwald fraction value being set. - -$\texttt{f_nd_plasma_pedestal_greenwald}$ and $\texttt{f_nd_plasma_separatrix_greenwald}$ can be set as iteration variables respectively by using `ixc = 45` -and `ixc = 152` respectively - ------- ### Pedestal Density Upper limit diff --git a/process/core/init.py b/process/core/init.py index d298b30f90..c567aa3c78 100644 --- a/process/core/init.py +++ b/process/core/init.py @@ -21,6 +21,7 @@ from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode from process.data_structure.physics_variables import DivertorNumberModels from process.data_structure.scan_variables import init_scan_variables +from process.models.physics.profiles import DensityProfilePedestalType from process.models.stellarator.initialization import st_init from process.models.superconductors import ( SuperconductorMaterial, @@ -421,41 +422,49 @@ def check_process(inputs, data): # noqa: ARG001 ) # Density checks - # Case where pedestal density is set manually + # Issue #589: Pedestal density is lower than separatrix density + pedestal_type = DensityProfilePedestalType( + data.physics.i_nd_plasma_pedestal_separatrix + ) if ( - data.physics.f_nd_plasma_pedestal_greenwald < 0 - or not ( - data_structure.numerics.ixc[: data_structure.numerics.nvar] == 145 - ).any() + pedestal_type == DensityProfilePedestalType.USER_INPUT + and data.physics.nd_plasma_pedestal_electron + < data.physics.nd_plasma_separatrix_electron + ) or ( + pedestal_type == DensityProfilePedestalType.GREENWALD_FRACTION + and data.physics.f_nd_plasma_pedestal_greenwald + < data.physics.f_nd_plasma_separatrix_greenwald ): - # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron but it is less than nd_plasma_separatrix_electron. - if ( - data.physics.nd_plasma_pedestal_electron - < data.physics.nd_plasma_separatrix_electron - ): - raise ProcessValidationError( - "Density pedestal is lower than separatrix density", - nd_plasma_pedestal_electron=data.physics.nd_plasma_pedestal_electron, - nd_plasma_separatrix_electron=data.physics.nd_plasma_separatrix_electron, - ) + raise ProcessValidationError( + "Density pedestal is lower than separatrix density", + **( + { + "nd_plasma_pedestal_electron": data.physics.nd_plasma_pedestal_electron, + "nd_plasma_separatrix_electron": data.physics.nd_plasma_separatrix_electron, + } + if pedestal_type == DensityProfilePedestalType.USER_INPUT + else { + "f_nd_plasma_pedestal_greenwald": data.physics.f_nd_plasma_pedestal_greenwald, + "f_nd_plasma_separatrix_greenwald": data.physics.f_nd_plasma_separatrix_greenwald, + } + ), + ) - # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron, - # but pedestal width = 0. - if ( - abs(data.physics.radius_plasma_pedestal_density_norm - 1.0) <= 1e-7 - and ( - data.physics.nd_plasma_pedestal_electron - - data.physics.nd_plasma_separatrix_electron - ) - >= 1e-7 - ): - warn( - "Density pedestal is at plasma edge " - f"({data.physics.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron " - f"({data.physics.nd_plasma_pedestal_electron}) differs from " - f"nd_plasma_separatrix_electron ({data.physics.nd_plasma_separatrix_electron})", - stacklevel=2, - ) + if ( + abs(data.physics.radius_plasma_pedestal_density_norm - 1.0) <= 1e-7 + and ( + data.physics.nd_plasma_pedestal_electron + - data.physics.nd_plasma_separatrix_electron + ) + >= 1e-7 + ): + warn( + "Density pedestal is at plasma edge " + f"({data.physics.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron " + f"({data.physics.nd_plasma_pedestal_electron}) differs from " + f"nd_plasma_separatrix_electron ({data.physics.nd_plasma_separatrix_electron})", + stacklevel=2, + ) # Issue #862 : Variable nd_plasma_electron_on_axis/nd_plasma_pedestal_electron ratio without constraint eq 81 (nd_plasma_electron_on_axis>nd_plasma_pedestal_electron) # -> Potential hollowed density profile diff --git a/process/core/input.py b/process/core/input.py index bfc45f88c3..1bb6a01323 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -149,9 +149,9 @@ def __post_init__(self): "f_p_div_lower": InputVariable("physics", float, range=(0.0, 1.0)), "f_plasma_fuel_deuterium": InputVariable("physics", float, range=(0.0, 1.0)), "ffwal": InputVariable("physics", float, range=(0.0, 10.0)), - "f_nd_plasma_pedestal_greenwald": InputVariable("physics", float, range=(-1.0, 5.0)), + "f_nd_plasma_pedestal_greenwald": InputVariable("physics", float, range=(0.1, 1.5)), "f_nd_plasma_separatrix_greenwald": InputVariable( - "physics", float, range=(-1.0, 5.0) + "physics", float, range=(0.001, 0.9) ), "f_plasma_fuel_helium3": InputVariable("physics", float, range=(-1.0, 5.0)), # TODO: does f_nd_impurity_electrons require an additional range? @@ -578,6 +578,9 @@ def __post_init__(self): ), "nd_plasma_pedestal_electron": InputVariable("physics", float, range=(0.0, 1e21)), "nd_plasma_separatrix_electron": InputVariable("physics", float, range=(0.0, 1e21)), + "i_nd_plasma_pedestal_separatrix": InputVariable( + data_structure.physics_variables, int, choices=[0, 1] + ), "nflutfmax": InputVariable("constraints", float, range=(0.0, 1e24)), "j_tf_coil_full_area": InputVariable("tfcoil", float, range=(10000.0, 1000000000.0)), "f_a_cs_turn_steel": InputVariable("pf_coil", float, range=(0.001, 0.999)), diff --git a/process/core/io/mfile/comparison.py b/process/core/io/mfile/comparison.py index c55e95d80f..9ed37ecae4 100644 --- a/process/core/io/mfile/comparison.py +++ b/process/core/io/mfile/comparison.py @@ -117,7 +117,7 @@ "temp_plasma_electron_on_axis_kev", "nd_plasma_electrons_vol_avg", "nd_plasma_electron_on_axis", - "dnla_gw", + "f_nd_plasma_greenwald", "temp_plasma_separatrix_kev", "nd_plasma_separatrix_electron", "temp_plasma_pedestal_kev", diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 410d3a3620..e3b4e684d6 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -3172,7 +3172,7 @@ def plot_main_plasma_information( f"$\\mathbf{{Density \\ limit:}}$\n" f"({DensityLimitModel(int(mfile.get('i_density_limit', scan=scan))).full_name})\n" f"$n_{{\\text{{e,limit}}}}: {mfile.get('nd_plasma_electrons_max', scan=scan):.3e} \\ m^{{-3}}$\n" - f"$f_{{\\text{{GW}}}}$: {mfile.get('dnla_gw', scan=scan):.4f}" + f"$f_{{\\text{{GW}}}}$: {mfile.get('f_nd_plasma_greenwald', scan=scan):.4f}" ) axis.text( @@ -3866,8 +3866,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): nd_ions_total = mfile.get("nd_plasma_ions_total_vol_avg", scan=scan) nd_fuel_ions = mfile.get("nd_plasma_fuel_ions_vol_avg", scan=scan) alphan = mfile.get("alphan", scan=scan) - fgwped_out = mfile.get("fgwped_out", scan=scan) - fgwsep_out = mfile.get("fgwsep_out", scan=scan) + f_nd_plasma_pedestal_greenwald = mfile.get( + "f_nd_plasma_pedestal_greenwald", scan=scan + ) + f_nd_plasma_separatrix_greenwald = mfile.get( + "f_nd_plasma_separatrix_greenwald", scan=scan + ) nd_plasma_electrons_vol_avg = mfile.get("nd_plasma_electrons_vol_avg", scan=scan) # find impurity densities imp_frac = np.array([ @@ -4057,7 +4061,10 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): # Add text box with density profile parameters textstr_density = "\n".join(( - rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$", + ( + rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$" + rf"$\hspace{{4}} \overline{{n_{{e}}}}$: {mfile.get('nd_plasma_electron_line', scan=scan):.3e} m$^{{-3}}$" + ), ( rf"$n_{{\text{{e,0}}}}$: {ne0:.3e} m$^{{-3}}$" rf"$\hspace{{4}} \alpha_{{\text{{n}}}}$: {alphan:.3f}" @@ -4068,7 +4075,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): f"{nd_fuel_ions / nd_plasma_electrons_vol_avg:.3f}" ), ( - rf"$f_{{\text{{GW e,ped}}}}$: {fgwped_out:.3f}" + rf"$f_{{\text{{GW e,ped}}}}$: {f_nd_plasma_pedestal_greenwald:.3f}" r"$ \hspace{7} \frac{n_{e,0}}{\langle n_e \rangle}$: " f"{ne0 / nd_plasma_electrons_vol_avg:.3f}" ), @@ -4078,12 +4085,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): f"{mfile.get('nd_plasma_electron_line', scan=scan) / mfile.get('nd_plasma_electron_max_array(7)', scan=scan):.3f}" ), rf"$n_{{\text{{e,sep}}}}$: {nd_plasma_separatrix_electron:.3e} m$^{{-3}}$", - rf"$f_{{\text{{GW e,sep}}}}$: {fgwsep_out:.3f}", + rf"$f_{{\text{{GW e,sep}}}}$: {f_nd_plasma_separatrix_greenwald:.3f}", )) props_density = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5} ax_main.text( - 0.0, + -0.05, -0.175, textstr_density, transform=ax_impurity.transAxes, @@ -4302,33 +4309,34 @@ def plot_t_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int): # Add text box with temperature profile parameters textstr_temperature = "\n".join(( ( - rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$: {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV" - rf"$\hspace{{3}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV" + rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$: {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV" + rf"$\hspace{{2}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV" + rf"$\hspace{{2}} \overline{{T_{{e}}}}$: {mfile.get('temp_plasma_electron_line_avg_kev', scan=scan):.3f} keV" ), ( - rf"$T_{{\text{{e,0}}}}$: {te0:.3f} keV" - rf"$\hspace{{4}} \alpha_{{\text{{T}}}}$: {alphat:.3f}" + rf"$T_{{\text{{e,0}}}}$: {te0:.3f} keV" + rf"$\hspace{{2}} \alpha_{{\text{{T}}}}$: {alphat:.3f}" ), ( rf"$T_{{\text{{e,ped}}}}$: {temp_plasma_pedestal_kev:.3f} keV" - r"$ \hspace{4} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: " + r"$ \hspace{3} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: " f"{f_temp_plasma_ion_electron:.3f}" ), ( rf"$\rho_{{\text{{ped,T}}}}$: {radius_plasma_pedestal_temp_norm:.3f}" - r"$ \hspace{6} \frac{T_{e,0}}{\langle T_e \rangle}$: " + r"$ \hspace{5} \frac{T_{e,0}}{\langle T_e \rangle}$: " f"{te0 / te:.3f}" ), ( rf"$T_{{\text{{e,sep}}}}$: {temp_plasma_separatrix_kev:.3f} keV" - r"$ \hspace{4} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: " + r"$ \hspace{3} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: " f"{mfile.get('f_temp_plasma_electron_density_vol_avg', scan=scan):.3f}" ), )) props_temperature = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5} prof.text( - 0.0, + -0.1, -0.125, textstr_temperature, transform=prof.transAxes, diff --git a/process/core/io/variable_metadata.py b/process/core/io/variable_metadata.py index b0bf93117d..aed8636100 100644 --- a/process/core/io/variable_metadata.py +++ b/process/core/io/variable_metadata.py @@ -187,7 +187,7 @@ class VariableMetadata: description="Average electron density", units="m^{-3}", ), - "dnla_gw": VariableMetadata( + "f_nd_plasma_greenwald": VariableMetadata( latex=r"$f_{\mathrm{GW}}$", description="Greenwald fraction", units="" ), "normalised_toroidal_beta": VariableMetadata( diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index a3add624a0..4b858c7a33 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -210,8 +210,8 @@ class IterationVariable: 1.0e17, 1.0e20, ), - 145: IterationVariable("f_nd_plasma_pedestal_greenwald", "physics", 0.1, 0.9), - 152: IterationVariable("f_nd_plasma_separatrix_greenwald", "physics", 0.001, 0.5), + 145: IterationVariable("f_nd_plasma_pedestal_greenwald", "physics", 0.1, 1.5), + 152: IterationVariable("f_nd_plasma_separatrix_greenwald", "physics", 0.001, 0.9), 155: IterationVariable("pfusife", "ife", 5.0e2, 3.0e3), 156: IterationVariable("rrin", "ife", 1.0, 1.0e1), 158: IterationVariable( diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 6ded110158..43af54e6a7 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -336,16 +336,16 @@ class PhysicsData: load calculation (`i_pflux_fw_neutron=1`) """ + f_nd_plasma_greenwald: float = None + """Greenwald fraction of the line averaged electron density. The classic Greenwald + limit value""" + f_nd_plasma_pedestal_greenwald: float = 0.85 - """fraction of Greenwald density to set as pedestal-top density. If `<0`, pedestal-top - density set manually using nd_plasma_pedestal_electron (`i_plasma_pedestal==1`). - (`iteration variable 145`) + """Greenwald fraction of the pedestal density """ - f_nd_plasma_separatrix_greenwald: float = 0.50 - """fraction of Greenwald density to set as separatrix density. If `<0`, separatrix - density set manually using nd_plasma_separatrix_electron (`i_plasma_pedestal==1`). - (`iteration variable 152`) + f_nd_plasma_separatrix_greenwald: float = 0.5 + """Greenwald fraction of the separatrix density """ f_plasma_fuel_helium3: float = 0.0 @@ -488,16 +488,22 @@ class PhysicsData: """ nd_plasma_pedestal_electron: float = 4.0e19 - """electron density of pedestal [m-3] (`i_plasma_pedestal==1)""" + """electron density of pedestal [m⁻³] (`i_plasma_pedestal==1)""" nd_plasma_separatrix_electron: float = 3.0e19 - """electron density at separatrix [m-3] (`i_plasma_pedestal==1)""" + """electron density at separatrix [m⁻³] (`i_plasma_pedestal==1)""" + + i_nd_plasma_pedestal_separatrix: int = 1 + """switch for pedestal and separatrix density calculation: + - =0 User input pedestal and separatrix density + - =1 Calculate pedestal and separatrix density as fraction of Greenwald limit (see `f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald`) + """ alpha_crit: float = 0.0 """critical ballooning parameter value""" nd_plasma_separatrix_electron_eich_max: float = 0.0 - """Eich critical electron density at separatrix [m-3]""" + """Eich critical electron density at separatrix [m⁻³]""" plasma_res_factor: float = 1.0 """plasma resistivity pre-factor""" @@ -983,6 +989,9 @@ class PhysicsData: temp_plasma_electron_density_weighted_kev: float = 0.0 """density weighted average electron temperature (keV)""" + temp_plasma_electron_line_avg_kev: float = None + """Line averaged electron temperature (keV)""" + temp_plasma_ion_vol_avg_kev: float = 12.9 """volume averaged ion temperature (keV). N.B. calculated from temp_plasma_electron_vol_avg_kev if `f_temp_plasma_ion_electron > 0.0`""" diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py index 2a5d486197..e917382ebf 100644 --- a/process/models/physics/density_limit.py +++ b/process/models/physics/density_limit.py @@ -87,7 +87,9 @@ def run(self): zeff=self.data.physics.n_charge_plasma_effective_vol_avg, ) - # Calculate beta_norm_max based on i_beta_norm_max + # Convert the chosen density limit model to the actual density limit value and + # store in physics variables. This is the value that will be used for all + # comparisons to the plasma density in the rest of the code. try: model = DensityLimitModel(int(self.data.physics.i_density_limit)) self.data.physics.nd_plasma_electrons_max = self.get_density_limit_value( @@ -100,6 +102,12 @@ def run(self): i_density_limit=self.data.physics.i_density_limit, ) from None + # Assign the Greenwald fraction for the rest of the code + self.data.physics.f_nd_plasma_greenwald = ( + self.data.physics.nd_plasma_electron_line + / self.data.physics.nd_plasma_electron_max_array[6] + ) + @staticmethod def get_density_limit_value( model: DensityLimitModel, nd_plasma_electron_max_array: np.ndarray diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 573ab73acf..118d749801 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -21,7 +21,10 @@ from process.data_structure.impurity_radiation_variables import N_IMPURITIES from process.models.physics import impurity_radiation from process.models.physics.plasma_geometry import PlasmaGeom -from process.models.physics.profiles import PlasmaProfileShapeType +from process.models.physics.profiles import ( + DensityProfilePedestalType, + PlasmaProfileShapeType, +) if TYPE_CHECKING: from process.data_structure.physics_variables import PhysicsData @@ -356,24 +359,8 @@ def run(self): if ( PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) == PlasmaProfileShapeType.PEDESTAL_PROFILE - ) and (self.data.physics.f_nd_plasma_pedestal_greenwald >= 0e0): - self.data.physics.nd_plasma_pedestal_electron = ( - self.data.physics.f_nd_plasma_pedestal_greenwald - * 1.0e14 - * self.data.physics.plasma_current - / (np.pi * self.data.physics.rminor * self.data.physics.rminor) - ) - - if ( - PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) - == PlasmaProfileShapeType.PEDESTAL_PROFILE - ) and (self.data.physics.f_nd_plasma_separatrix_greenwald >= 0e0): - self.data.physics.nd_plasma_separatrix_electron = ( - self.data.physics.f_nd_plasma_separatrix_greenwald - * 1.0e14 - * self.data.physics.plasma_current - / (np.pi * self.data.physics.rminor * self.data.physics.rminor) - ) + ): + self.plasma_profile.neprofile.set_pedestal_and_separatrix_values() self.plasma_profile.run() @@ -2461,6 +2448,12 @@ def output_temperature_density_profile_info(self) -> None: self.data.physics.temp_plasma_electron_on_axis_kev, "OP ", ) + po.ovarrf( + self.outfile, + "Line averaged electron temperature (keV)", + "(temp_plasma_electron_line_avg_kev)", + self.data.physics.temp_plasma_electron_line_avg_kev, + ) po.ovarrf( self.outfile, "Volume averaged density weighted electron temperature (⟨Tₑ⟩ₙ) (keV)", @@ -2509,68 +2502,90 @@ def output_temperature_density_profile_info(self) -> None: ) po.oblnkl(self.outfile) if ( - PlasmaProfileShapeType(self.data.physics.i_plasma_pedestal) + self.data.physics.i_plasma_pedestal == PlasmaProfileShapeType.PEDESTAL_PROFILE ): + po.ovarin( + self.outfile, + "Pedestal and separatrix density model selected", + "(i_nd_plasma_pedestal_separatrix)", + self.data.physics.i_nd_plasma_pedestal_separatrix, + ) + po.ocmmnt( + self.outfile, + f"Pedestal and separatrix values set by: " + f"{DensityProfilePedestalType(self.data.physics.i_nd_plasma_pedestal_separatrix).description}", + ) + po.oblnkl(self.outfile) + 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( + if ( + self.data.physics.i_nd_plasma_pedestal_separatrix + == DensityProfilePedestalType.USER_INPUT + ): + po.ovarin( self.outfile, "Electron density pedestal height (nₑ_pedestal) (/m³)", "(nd_plasma_pedestal_electron)", self.data.physics.nd_plasma_pedestal_electron, + ) + po.ovarin( + self.outfile, + "Electron separatrix density (nₑ,ₛₑₚ) (/m³)", + "(nd_plasma_separatrix_electron)", + self.data.physics.nd_plasma_separatrix_electron, + ) + po.ovarre( + self.outfile, + "Pedestal Greenwald fraction", + "(f_nd_plasma_pedestal_greenwald)", + self.data.physics.f_nd_plasma_pedestal_greenwald, "OP ", ) - else: + po.ovarre( + self.outfile, + "Separatrix Greenwald fraction", + "(f_nd_plasma_separatrix_greenwald)", + self.data.physics.f_nd_plasma_separatrix_greenwald, + "OP ", + ) + + elif ( + self.data.physics.i_nd_plasma_pedestal_separatrix + == DensityProfilePedestalType.GREENWALD_FRACTION + ): po.ovarre( self.outfile, "Electron density pedestal height (nₑ_pedestal) (/m³)", "(nd_plasma_pedestal_electron)", self.data.physics.nd_plasma_pedestal_electron, + "OP ", ) - # 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] + po.ovarre( + self.outfile, + "Electron separatrix density (nₑ,ₛₑₚ) (/m³)", + "(nd_plasma_separatrix_electron)", + self.data.physics.nd_plasma_separatrix_electron, + "OP ", ) - 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.ovarin( + self.outfile, + "Pedestal Greenwald fraction", + "(f_nd_plasma_pedestal_greenwald)", + self.data.physics.f_nd_plasma_pedestal_greenwald, + ) + po.ovarin( + self.outfile, + "Separatrix Greenwald fraction", + "(f_nd_plasma_separatrix_greenwald)", + self.data.physics.f_nd_plasma_separatrix_greenwald, ) - 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( @@ -2597,9 +2612,8 @@ def output_temperature_density_profile_info(self) -> None: 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], + "(f_nd_plasma_greenwald)", + self.data.physics.f_nd_plasma_greenwald, "OP ", ) po.oblnkl(self.outfile) diff --git a/process/models/physics/plasma_profiles.py b/process/models/physics/plasma_profiles.py index 79e2027596..b9ca7cfd42 100644 --- a/process/models/physics/plasma_profiles.py +++ b/process/models/physics/plasma_profiles.py @@ -140,6 +140,14 @@ def parabolic_paramterisation(self): / sp.special.gamma(self.data.physics.alphan + 1.5) ) + self.data.physics.temp_plasma_electron_line_avg_kev = ( + self.data.physics.temp_plasma_electron_vol_avg_kev + * (1.0 + self.data.physics.alphat) + * (sp.special.gamma(0.5) / 2.0) + * sp.special.gamma(self.data.physics.alphat + 1.0) + / sp.special.gamma(self.data.physics.alphat + 1.5) + ) + # Density-weighted temperatures self.data.physics.temp_plasma_electron_density_weighted_kev = ( @@ -151,7 +159,7 @@ def parabolic_paramterisation(self): * self.data.physics.f_temp_plasma_electron_density_vol_avg ) - # Central values for temperature (keV) and density (m**-3) + # Central values for temperature (keV) and density (m^-3) self.data.physics.temp_plasma_electron_on_axis_kev = ( self.data.physics.temp_plasma_electron_vol_avg_kev @@ -189,9 +197,9 @@ def pedestal_parameterisation(self): # Perform integrations to calculate ratio of density-weighted # to volume-averaged temperature, etc. - # Density-weighted temperature = integral(n.T dV) / integral(n dV) + # Density-weighted temperature = ∫(n.T dV) / ∫(n dV) # which is approximately equal to the ratio - # integral(rho.n(rho).T(rho) drho) / integral(rho.n(rho) drho) + # ∫(ρ.n(ρ).T(ρ) dρ) / ∫(ρ.n(ρ) dρ) # noqa: RUF003 drho = self.neprofile.profile_dx dens = self.neprofile.profile_y @@ -219,11 +227,15 @@ def pedestal_parameterisation(self): / self.data.physics.temp_plasma_electron_vol_avg_kev ) - # Line-averaged electron density - # = integral(n(rho).drho) + # Line-averaged electron density and temperature + # = ∫(n(ρ).dρ) # noqa: RUF003 self.data.physics.nd_plasma_electron_line = self.neprofile.profile_integ + self.data.physics.temp_plasma_electron_line_avg_kev = ( + self.teprofile.profile_integ + ) + # Scrape-off density / volume averaged density # (Input value is used if i_plasma_pedestal = 0) @@ -283,7 +295,7 @@ def calculate_profile_factors(self): # Pressure profile index (only true for a parabolic profile) # N.B. pres_plasma_thermal_on_axis is NOT equal to

* (1 + alphap), - # but p(rho) = n(rho)*T(rho) + # but p(ρ) = n(ρ)*T(ρ) # noqa: RUF003 # and

= .T_n where <...> denotes volume-averages and T_n is the # density-weighted temperature @@ -299,7 +311,7 @@ def calculate_profile_factors(self): * self.data.physics.temp_plasma_ion_density_weighted_kev ) * constants.KILOELECTRON_VOLT - # Central plasma current density (A/m^2) + # Central plasma current density (A/m²) # Assumes a parabolic profile for the current density self.data.physics.j_plasma_on_axis = ( (self.data.physics.plasma_current) diff --git a/process/models/physics/profiles.py b/process/models/physics/profiles.py index 3fc53ccaea..5b46410ee2 100644 --- a/process/models/physics/profiles.py +++ b/process/models/physics/profiles.py @@ -13,6 +13,7 @@ import scipy as sp from process.core.model import Model +from process.models.physics.density_limit import PlasmaDensityLimit logger = logging.getLogger(__name__) @@ -109,6 +110,31 @@ def integrate_profile_y(self): ) +class DensityProfilePedestalType(IntEnum): + """Enum for i_nd_plasma_pedestal_separatrix types""" + + USER_INPUT = (0, "User input direct values") + GREENWALD_FRACTION = (1, "Fractions of the Greenwald limit") + + def __new__(cls, value: int, description: str): + """Create a new DensityProfilePedestalType instance. + + Parameters + ---------- + value: Integer value for the enum member. + description: Human-readable description for the enum member. + """ + obj = int.__new__(cls, value) + obj._value_ = value + obj._description_ = description + return obj + + @DynamicClassAttribute + def description(self): + """Get the description of the plasma profile shape.""" + return self._description_ + + class NeProfile(Profile): """Electron density profile class. Contains a function to calculate the electron density profile and store the data. @@ -255,6 +281,51 @@ def ncore( ncore = 1.0e-6 return ncore + def set_pedestal_and_separatrix_values(self): + """Sets the pedestal and separatrix density values based on the user input or greenwald fraction method.""" + i_nd_plasma_pedestal_separatrix = DensityProfilePedestalType( + self.data.physics.i_nd_plasma_pedestal_separatrix + ) + + if i_nd_plasma_pedestal_separatrix == DensityProfilePedestalType.USER_INPUT: + self.data.physics.f_nd_plasma_pedestal_greenwald = ( + self.data.physics.nd_plasma_pedestal_electron + / ( + PlasmaDensityLimit.calculate_greenwald_density_limit( + c_plasma=self.data.physics.plasma_current, + rminor=self.data.physics.rminor, + ) + ) + ) + + self.data.physics.f_nd_plasma_separatrix_greenwald = ( + self.data.physics.nd_plasma_separatrix_electron + / ( + PlasmaDensityLimit.calculate_greenwald_density_limit( + c_plasma=self.data.physics.plasma_current, + rminor=self.data.physics.rminor, + ) + ) + ) + elif ( + i_nd_plasma_pedestal_separatrix + == DensityProfilePedestalType.GREENWALD_FRACTION + ): + self.data.physics.nd_plasma_pedestal_electron = ( + self.data.physics.f_nd_plasma_pedestal_greenwald + * PlasmaDensityLimit.calculate_greenwald_density_limit( + c_plasma=self.data.physics.plasma_current, + rminor=self.data.physics.rminor, + ) + ) + self.data.physics.nd_plasma_separatrix_electron = ( + self.data.physics.f_nd_plasma_separatrix_greenwald + * PlasmaDensityLimit.calculate_greenwald_density_limit( + c_plasma=self.data.physics.plasma_current, + rminor=self.data.physics.rminor, + ) + ) + def set_physics_variables(self): """Calculates and sets physics variables required for the profile.""" if ( diff --git a/tests/unit/models/physics/test_plasma_profiles.py b/tests/unit/models/physics/test_plasma_profiles.py index e836abb508..57383c32c0 100644 --- a/tests/unit/models/physics/test_plasma_profiles.py +++ b/tests/unit/models/physics/test_plasma_profiles.py @@ -228,6 +228,8 @@ class PlasmaProfilesParam(NamedTuple): expected_nd_electron_line: float = 0.0 + expected_temp_plasma_electron_line_avg_kev: float = 0.0 + expected_ti: float = 0.0 @@ -277,6 +279,7 @@ class PlasmaProfilesParam(NamedTuple): expected_ne0=1.0585658890823703e20, expected_ti0=27.370104119511087, expected_nd_electron_line=8.8687354645836431e19, + expected_temp_plasma_electron_line_avg_kev=17.582796426026842, expected_ti=13.07, ), PlasmaProfilesParam( @@ -322,6 +325,7 @@ class PlasmaProfilesParam(NamedTuple): expected_ne0=1.0585658890823703e20, expected_ti0=27.370104119511087, expected_nd_electron_line=8.8687354645836431e19, + expected_temp_plasma_electron_line_avg_kev=17.582796426026842, expected_ti=13.07, ), ], @@ -548,6 +552,14 @@ def test_plasma_profiles(plasmaprofilesparam, monkeypatch, plasmaprofile): plasmaprofilesparam.expected_nd_electron_line ) + assert plasmaprofile.data.physics.temp_plasma_electron_line_avg_kev == pytest.approx( + plasmaprofilesparam.expected_temp_plasma_electron_line_avg_kev + ) + + assert plasmaprofile.data.physics.temp_plasma_electron_line_avg_kev == pytest.approx( + plasmaprofilesparam.expected_temp_plasma_electron_line_avg_kev + ) + assert plasmaprofile.data.physics.temp_plasma_ion_vol_avg_kev == pytest.approx( plasmaprofilesparam.expected_ti )