Warnings for function ncore

[jonmaddock]

In GitLab by @stuartmuldrew on Jul 9, 2018, 08:32

The function ncore in plasma_profiles.f90 calculates the central density of a pedestal profile, however there is currently no warning if the central density is less than the pedestal density. Additionally, under certain situations it returns a negative central density without stopping. I will add a warning to the error list.

[jonmaddock]

In GitLab by @stuartmuldrew on Jul 10, 2018, 13:12

created branch 735-ncore

[jonmaddock]

In GitLab by @stuartmuldrew on Jul 10, 2018, 14:09

mentioned in commit 90f7ddd5ddac024534591074ee09ef7a2b9d4003

[jonmaddock]

In GitLab by @stuartmuldrew on Jul 10, 2018, 14:44

mentioned in commit 18af683e2aa401094b5940ef73071087d3c9cdc1

[jonmaddock]

In GitLab by @stuartmuldrew on Jul 10, 2018, 17:01

mentioned in commit e9a3810ba87f8859447297072f733a9322d0008a

[jonmaddock]

In GitLab by @mkovari on Jul 17, 2018, 11:51

I don't think it should be possible for the central density to be less than the pedestal density unless the density profile exponent has a silly value. Is there a bug in ncore?

[jonmaddock]

In GitLab by @stuartmuldrew on Aug 24, 2018, 09:07

The equation for the core density is given by:

ncore = 1.0D0 / (3.0D0*rhopedn**2) * ( &
         3.0D0*nav*(1.0D0 + alphan) &
         + nsep*(1.0D0 + alphan)*(-2.0D0 + rhopedn + rhopedn**2) &
         - nped*( (1.0D0 + alphan)*(1.0D0 + rhopedn) + &
         (alphan - 2.0D0)*rhopedn**2 ) )

I guess if the pedestal density is set to be very high in the input compared to the volume averaged, then the core density can end up low. I've added the warning so maybe this will inform the user not to do that.

The reference for that equation is the Johner (2011) HELIOS paper, however it doesn't appear in that form. In the paper the central density is written in a very convoluted way.

[jonmaddock]

In GitLab by @stuartmuldrew on Oct 9, 2018, 16:03

Attached is Hanni's Mathematica file ped_averages.nb

[jonmaddock]

In GitLab by @mkovari on Oct 10, 2018, 08:50

Can you attach this file as an image?

[jonmaddock]

In GitLab by @stuartmuldrew on Oct 24, 2018, 16:59

Unfortunately that is the only file Hanni could find and I don't have access to Mathematica to run it. Is there a machine with Mathematica installed?

[jonmaddock]

In GitLab by @mkovari on Oct 24, 2018, 17:05

Not that I know of

[jonmaddock]

In GitLab by @skahn on Jun 10, 2019, 11:46

Hi all

Investigation results

It is really possible for Hanni's formula to be right and the issue to be caused by the fact that the nav (volumic average electrons temperature also dene) value of the iteration is too low.

It appears that even though raising the flag and the appropriate action is commented in the terminal :
The 2 CHECK: boundl(6) has been raised to ensure dene > neped
nothing is really done in the value, and therefore dene is at some point low enough to get a negative central density value.

Proposed solutions

It seems that all of these problems comes from the fact that the pedestal and the separatrix densities are defined in a different than the average density value :

• The average electronic density is a optimization variable with only an upper limit using fGW on it.
• The separatrix and the pedestal electronic densities are defined with an exact fraction of the greenwald limit

In this situation, the pedestal densities can be set at a relatively high value while the average densities are free to vary to lower values, causing from time to time negative central densities.

1. Find a way to really enforce the action printed in the warning messages
2. Rethink the way the density profile is parametrized to make sure that the density profile keeps the desired shape during the optimization process.
3. Add an inequality constraint equation make sure ne(0) > ne(ped) (already done in local version)

Solution 1

The constraint or the warning check is performed with the initial value of neped. If the neped is set with the greenwald fraction, the defauld PROCESS value is used for the test : 0.4E20.
If there is no lower boundary, there will be effectively a lower bound of 0.404E20.
In the STEP run, nesep is set with 0.65*nGW ~ 1.8E20. Therefore this security is not enough.

The solution would be to have a dynamic boundary setting for the average density, by the structure of the code is currently not adapted to this. The solution 3 therefore may be better for now.

Solution 2 proposal

If solution 2 is proposed, it can be parametrized as

• n0/nped set to a number
• nped/nsep set to a number
• average densities set with the Greenwald fractions.

All the three would be set as exact values or iteration variables, with properly defined limit ranges

Solution 3

This is done in a local version, and although this is not perfectly avoiding the negative ne(0) values, it clearly improves the situation

Cheers

Seb

[jonmaddock]

In GitLab by @jmorris-uk on Jul 17, 2019, 18:03

removed milestone

[jonmaddock]

In GitLab by @skahn on Jul 29, 2019, 09:24

This has lead to a contraint, we can close this issue.

[jonmaddock]

In GitLab by @skahn on Jul 29, 2019, 09:24

closed

[jonmaddock]

In GitLab by @mkovari on Jul 29, 2019, 09:26

Hi Seb,
Are you saying that no action is required so the issue can be closed? If so, why?

[jonmaddock]

In GitLab by @skahn on Jul 29, 2019, 09:28

I am saying that action has been taken (solution 1), and fully implemented in PROCESS develop.

[jonmaddock]

In GitLab by @jmorris-uk on Jul 30, 2019, 14:06

Constraint implemented as constraint number 81. See attached image. Constraint enforces ne0 to be above neped.

[jonmaddock]

In GitLab by @mkovari on Jul 30, 2019, 14:08

I assume that you have since changed the pasted comments in the code.

[jonmaddock]

In GitLab by @jmorris-uk on Jul 30, 2019, 14:30

... as of 5 seconds ago.

[mkovari]

This problem raised its head again in December 2022:

The graph shows ncore as a function of nsep, for different values of nped.
n.0 does indeed go negative, but only when nped > n.av

n.0 is greater than n.av except when n.ped=n.av

The whole problem arises because we have the option of setting n.ped and n.sep in terms of the Greenwald density, but we don’t do the same with n.av.

There is a minus sign here.

You may notice that my paper doesn’t give any references for this formula, but I don’t know why. We may need to ask Hanni for the reference.

I though the reference might be the attached paper. However, the formula there looks very different:

where Theta1 and Gamma1 are given by separate formulae.

Johner_2010_HELIOS.pdf

[mkovari]

This issue was originally raised by Stuart but it doesn't seem that Seb fixed it fully since Jon discovered the same problem in December 2022. I have reopened it, leaving it as an issue (rather than a discussion) as it is a bug rather than a feature request. @jonmaddock @jmorris-uk

[chris-ashe]

I have made a quick plotting tool that allows all of the profile and pedestal values to be played around with so that it is easier to see how the negative central density problem comes about.
https://www.desmos.com/calculator/l9d9vtuzqe

A quick look from a file that presents such warnings shows that raising the lower limit of dene quickly pushes the solver away from a low enough value that ncore will never become negative.

A possible quick fix solution would be just to warn the user to raise the lower limit of dene

[mkovari]

Can you suggest a value for the lower bound of dene? For example, if you expect that the converged solution will have volume-averaged electron density dene~1e20 m-3 then set the lower bound to 8e19?

Also recommend that people use constraint number 81.

[chris-ashe]

@mkovari
Can easily be done. Though will need to have some margin as $n_{\text{ped}}$ and $n_{\text{sep}}$ can be iterated as well during a run and have influence on ncore. It would be best to present it as a warning and not enforce it internally to prevent how the profiles are solved becoming a black box.

Will also have a look at the ncore value that is enforced at the moment as it is re-set to 1E-6 if it becomes negative. Will see if pushing it back up to normal central density values pushes it back into the right direction for solving

[mkovari]

So can you suggest the lower bound as a fraction of the expected value of dene? I just gave 80% as an example.

I definitely don't like the 1E-6 kludge.

[chris-ashe]

Two possible solutions at the moment, one quick fix and one long term

Quick Fix
If ncore ever goes to zero or negative raise a level 2/3 error asking the user the please raise their lower limit of dene to prevent ncore going negative. This still allows full control and prevents implementing a black box solution

Long term
Invert the solving of ncore and have it as an iteration variable instead of dene . With dene instead being calculated for. Thus the lower bound of ncore could be specified to never go near zero. This will need to be stated clearly as ncore and dene cannot be iteration variables at the same time

@mkovari

[mkovari]

Qick fix - yes, so what lower bound do you suggest?

Other options

1. Allow ncore as an iteration variable instead of dene.

You pointed out that we often want to match a published scenario with specified volume-averaged density, so this option may not be ideal. However, it is certainly an option, and if we really want to we could add a constraint for a specified value of dene.

2. A better option would be to allow the ratio of ncore to neped as an iteration variable, with default lower bound = 1.

3. Or, allow the ratio of dene to neped as an iteration variable, with default lower bound = 1.

4. Or, specify all the densities as a parameter times the Greenwald density, and do not allow any to be iteration variables. This might interfere with convergence.
   $nesep = A \times N_{GW}$
   $neped = B \times N_{GW}$
   $dene = C \times N_{GW}$

5. Or, the previous option, but with an additional iteration variable f:
   $nesep = fA \times N_{GW}$
   $neped = fB \times N_{GW}$
   $dene = fC \times N_{GW}$

[mkovari]

We do already have a check at initialisation:

PROCESS/source/fortran/initial.f90

Lines 431 to 436 in a422025

	! Issue #862 : Variable ne0/neped ratio without constraint eq 81 (ne0>neped)
	! -> Potential hollowed density profile
	if ( (ioptimz >= 0) .and. (.not.any(icc==81)) ) then
	if ( any(ixc == 145 )) call report_error(154)
	if ( any(ixc == 6 )) call report_error(155)
	end if

This prints the following messages:

      "no": 154,
      "level": 2,
      "message": "CHECK: neped set with fgwped without constraint eq 81 (neped<ne0)"

      "no": 155,
      "level": 2,
      "message": "CHECK: dene used as iteration variable without constraint 81 (neped<ne0)"

[mkovari]

Compare this density issue with a similar problem that can arise with temperature, for which we do have a check at initialisation that actually changes the value of te and its lower bound:

PROCESS/source/fortran/initial.f90

Lines 384 to 410 in a422025

	! Temperature checks
	if (teped < tesep) then
	fdiags(1) = teped ; fdiags(2) = tesep
	call report_error(146)
	end if
	if ((abs(rhopedt-1.0D0) <= 1.0D-7).and.((teped-tesep) >= 1.0D-7)) then
	fdiags(1) = rhopedt ; fdiags(2) = teped ; fdiags(3) = tesep
	call report_error(147)
	end if
	! Core temperature should always be calculated (later) as being
	! higher than the pedestal temperature, if and only if the
	! volume-averaged temperature never drops below the pedestal
	! temperature. Prevent this by adjusting te, and its lower bound
	! (which will only have an effect if this is an optimisation run)
	if (te <= teped) then
	fdiags(1) = te ; fdiags(2) = teped
	te = teped*1.001D0
	call report_error(149)
	end if
	if ((ioptimz >= 0).and.(any(ixc == 4)).and.(boundl(4) < teped*1.001D0)) then
	call report_error(150)
	boundl(4) = teped*1.001D0
	boundu(4) = max(boundu(4), boundl(4))
	end if

[mkovari]

It is true that there are papers claiming that stable H mode can be obtained as long as ne,sep  ≤  nGW/2. However, this is the outer midplane density. I don't think this is representative of the average density over the separatrix. This is not something we have ever discussed, really. But if you need to then obviously go ahead and make sure the default bounds are high enough to permit the regression test to run without modification.

As you recall there are several other options available.