diff --git a/documentation/proc-pages/physics-models/plasma.md b/documentation/proc-pages/physics-models/plasma.md index fff2c5ee3d..9ae5863f15 100644 --- a/documentation/proc-pages/physics-models/plasma.md +++ b/documentation/proc-pages/physics-models/plasma.md @@ -267,6 +267,11 @@ simulations. If `ipedestal` = 1 or 2 then the pedestal density `neped` is set as a fraction `fgwped` of the Greenwald density (providing `fgwped` >= 0). The default value of `fgwped` is 0.8[^7]. +!!! warning " Un-realistic profiles" + + If `ipedestal >= 1` it is highly recommended to use constraint equation 81 (icc=81). This enforces solutions in which $n_0$ has to be greater than $n_{ped}$. + Negative $n_0$ values can also arise during iteration, so it is important to be weary on how low the lower bound for $n_e (\mathtt{dene})$ is set. + ## Beta Limit The plasma beta limit[^7] is given by diff --git a/process/profiles.py b/process/profiles.py index 339c30f763..5e8947048b 100644 --- a/process/profiles.py +++ b/process/profiles.py @@ -3,10 +3,7 @@ from scipy import integrate from abc import ABC, abstractmethod -from process.fortran import ( - maths_library, - physics_variables, -) +from process.fortran import maths_library, physics_variables, error_handling logger = logging.getLogger(__name__) # Logging handler for console output @@ -127,7 +124,9 @@ def calculate_profile_y(self, rho, rhopedn, n0, nped, nsep, alphan): ) @staticmethod - def ncore(rhopedn, nped, nsep, nav, alphan): + def ncore( + rhopedn: float, nped: float, nsep: float, nav: float, alphan: float + ) -> float: """This routine calculates the core denesity of a pedestalised profile. :param rhopedn: normalised minor radius pedestal position @@ -138,11 +137,10 @@ def ncore(rhopedn, nped, nsep, nav, alphan): :type nsep: float :param nav: electron density (/m3) :type nav: float - :param alphan: _description_ :param alphan: density peaking parameter :type alphan: float :return: Core density - :rtype: numpy.array + :type: float """ ncore = ( @@ -155,15 +153,11 @@ def ncore(rhopedn, nped, nsep, nav, alphan): ) ) - if ncore < 0: - # Prevent ncore from going negative (and terminating the optimisation) by - # kludging to small positive value. Allows solver to continue and - # hopefully be constrained away from this point (e.g. constraint 81, ne0 > neped) - logger.exception("ncore is negative. Kludging to 1e-6.") + if ncore < 0.0: + # Allows solver to continue and + # warns the user to raise the lower bound on dene if the run did not converge + error_handling.report_error(282) ncore = 1.0e-6 - # raise ValueError( - # f"Core density is negative: {ncore}. {nped = }, {nsep = }, {nav = }" - # ) return ncore def set_physics_variables(self): diff --git a/process/utilities/errorlist.json b/process/utilities/errorlist.json index 7aa7d672be..d39f74ed14 100644 --- a/process/utilities/errorlist.json +++ b/process/utilities/errorlist.json @@ -8,7 +8,7 @@ "comment2": [ "Increment n_errortypes if an error is added to this list" ], - "n_errortypes": 281, + "n_errortypes": 282, "errors": [ { "no": 1, @@ -1414,6 +1414,11 @@ "no": 281, "level": 3, "message": "CHECK: Cannot have i_tf_bucking >= 2 when tf_in_cs = 1" + }, + { + "no": 282, + "level": 2, + "message": "CHECK: ncore is going negative when solving. Please raise the value of dene and or its lower limit." } ] } diff --git a/source/fortran/iteration_variables.f90 b/source/fortran/iteration_variables.f90 index 72df52766f..08070f13b8 100755 --- a/source/fortran/iteration_variables.f90 +++ b/source/fortran/iteration_variables.f90 @@ -148,7 +148,7 @@ subroutine init_itv_6 use numerics, only: lablxc, boundl, boundu implicit none lablxc(6) = 'dene ' - boundl(6) = 1.00D19 + boundl(6) = 2.00D19 boundu(6) = 1.00D21 end subroutine init_itv_6 @@ -3244,8 +3244,8 @@ subroutine init_itv_145 use numerics, only: lablxc, boundl, boundu implicit none lablxc(145) = 'fgwped ' - boundl(145) = 0.500D0 - boundu(145) = 1.000D0 + boundl(145) = 0.100D0 + boundu(145) = 0.9D0 end subroutine init_itv_145 real(kind(1.d0)) function itv_145() @@ -3377,7 +3377,7 @@ subroutine init_itv_152 implicit none lablxc(152) = 'fgwsep ' boundl(152) = 0.001D0 - boundu(152) = 1.000D0 + boundu(152) = 0.5D0 end subroutine init_itv_152 real(kind(1.d0)) function itv_152()