gs2 initialization failure with gridgen error using geqdsk input

[Jason Francis Parisi]

When using efit_eq = .true. (reading in a geqdsk), initialization sometimes fails.

I have seen several error types, the most common of which states “Bad behaviour spotted in theta_grid_gridgen”.

I have replicated this error on two machines: Princeton’s Portal and Stellar clusters, so I don’t think it’s device-dependent.

Sometimes, increasing ntheta or changing equal_arc = true/false fixes the error, but often not. My current workaround is to run a script that changes ntheta / equal_arc until gs2 initializes successfully.

I have also tried changing input parameters epsknob, cv_fraction, and deltaw, but to no avail. Compiling and running with debug = .true., we don’t get any additional useful error messages.

I have included an input file, geqdsk file, a slurm submission file, and a slurm .out file, for an erroneous case.

[Jason Francis Parisi]

  <div class="preview-container wiki-content"><!-- loaded via ajax --></div>
  <div class="mask"></div>
</div>

</div> </form>

[Jason Francis Parisi]

I have attached a python script (and its results) that calculates which ntheta values gridgen fails for, for this particular input file and equilibrium. There is clustering of successful and unsuccessful ntheta values at lower ntheta, although at higher (ntheta>~100) ntheta values, gridgen almost always fails.

[David Dickinson]

Yes gridgen is a bit temperamental at high ntheta. I started work on a replacement a while ago but didn’t get around to finishing it off. In the meantime there are a few options to try (most of which I think you’re probably aware of and have tested):

1. set alknob = 1.0 (or much higher) – this adds more significance to the pitch angle resolution metrics, removing the dominance of the theta resolution metrics.
2. Try playing with epsknob and bpknob.
3. Set skip_gridgen = T (subject to some assumptions being true).
4. Use the file theta grid type to specify the geometry (e.g. adapt something like https://github.com/rahulgaur104/VMEC2GK for geqdsk).

‌

I’ll try and run your specific case to see if I can identify the source of the problem.

[David Dickinson]

You’re looking quite close to the edge and the geqdsk file isn’t the highest resolution. Using matplotlib to extract the psiN=0.96 flux surface gives Bmag vs theta:

The orange dot represents the minimum (not at theta=0, I’ll have to check if this is an assumption). Note also the wiggles on the outboard side. For me this comes from Bpol:

Clearly my quick python isn’t doing the same thing as GS2, but I think it probably highlights some of the challenges that the geometry calculation and gridgen may be facing.

[Jason Francis Parisi]

Thank you David. These B, Bpol (theta) plots are indeed concerning. The equilibria are generated for NSTX-like equilibria, where |B| is not necessary a minimum at theta = 0.

I noticed when running skip_gridgen = T, that phi → NaN; could you comment when skip_gridgen = T is suitable?

I will explore 1),2),4) more.

[David Dickinson]

Ah @Jason Parisi it seems skipping gridgen requires an assumption of symmetry :

When active we are effectively forced to assumed that the input bmag is symmetric and monotonic. If this is not true then we should not trust the results.

I think we probably don’t assume theta=0 corresponds to min(B), rather than max(B) is at the ends of the grid. I’ll have a dig when I get time.

[David Dickinson]

Just to confirm, but with skip_gridgen=T (so that I can dump the grids) I find bmag looks like the below as a function of theta grid index.

I think the reason for NaNs in the case we skip gridgen is that the resulting pitch angle grid is not monotonic.

Here’s cvdrift as a function of theta. It’s hard to imagine that gridgen would be able to improve this structure significantly. I suspect we may be hitting up against the limits of the equilibrium resolution here.

‌

[David Dickinson]

Here’s what we get for bmag with ntheta=512, whilst it technically worked I’m not sure you’d be happy with this grid!

‌

[Jason Francis Parisi]

@David Dickinson Thanks for investigating this in so much detail.

I think that the solution is for me to generate geqdsk files with higher resolution. I have been using a relatively new tool – EFIT-AI – which may require more refined user input on my end.

I will also compare this with CGYRO, to see how it handles comparably noisy geometric quantities.

Strange that cvdrift reasonable for |theta| ~< pi, and then goes crazy beyond that.

[David Dickinson]

I’d certainly be interested to hear how well CGYRO copes with the same equilibrium file!

‌

[David Dickinson]

Just to note, I’ve revisited this slightly following some of the recent changes in PR #944. Whilst the same issue is still occurring I have a better handle on the logic leading to this now, I think.

1. If we end up with the nthetaout is not odd error then this indicates that we have removed a set of theta points (corresponding to a unique value of B) which contains an odd number of points. Assuming the input B is periodic, a set of points with an odd number of members must include (an odd number of ) local extrema. We bias towards keeping such sets by giving them an artificially high resolution metric. We should only decide to delete such sets when we have already removed all the other sets – i.e. all those sets without a local extrema (or with an even number of local extrema).
2. For each set containing one or more local extrema we can count how many theta grid points we have here. We could (but don’t yet) warn if the sum of these point counts exceeds the requested number of theta grid points as this is likely to lead to the odd nthetaout error.
3. Simply increasing ntheta sometimes doesn’t help because it can increase the number of local extrema found in the B(theta) spline.

Some possible solutions include:

1. Smooth the input B(theta). We don’t currently offer a smoothing spline so this might require some additional code.
2. Increase the spline tension – increasing this discourages the spline from containing new local extrema – large tension makes the spline effectively a linear interpolation.

A secondary problem is where we remove too many points such that our resulting theta grid has fewer points than requested. This happens because we choose which set of points to remove purely by looking at the resolution metric. We could also count the number of points in the set and use this as an additional condition:

In other words replace idel = minimum_index(resolution_metric, bmagset /= 0.0, precision_bound) with idel = minimum_index(resolution_metric, bmagset /= 0.0 .and. npts_set <= ntheta_left-ntheta, precision_bound). Currently I’ve just added an onscreen warning when we are about to remove a set with more than ntheta_left - ntheta points.