Fix EQDSK geometry coupling and add a new ITER hybrid test case.

[copybara-service]

Fix EQDSK geometry coupling and add a new ITER hybrid test case.

New ITERhybrid EQDSK file generated from CHEASE data. Good comparison with the
simulation using CHEASE geometry.

Summary of fixes:

1. vpr definition bug leading to major impact on plasma volume
2. Avoiding LCFS flux-surface-integral issues for diverted cases by taking the last contour just inside the LCFS.
3. Rmaj taken now as Rgeo (LCFS Rmaj), and local Rmaj and Rmin used for delta (triangularity) calculations
4. StandardGeometryIntermediate values at the magnetic axis are prescribed, since a contour cannot be defined there.

The number of surfaces to calculate, and the location of the LCFS contour with respect to the boundary psi, are user-inputs, with reasonable defaults set but need to be checked against more varied sets of EQDSK files.

The EQDSK converter has only been tested against CHEASE-generated EQDSK which is at a specific COCOS. Therefore issues may arise if the input EQDSK is in a different COCOS. In near-future work we should demand that the EQDSK COCOS is an input, and we use the eqdsk convertor to convert to the TORAX COCOS, which itself also needs to be made fully self-consistent in a future PR.

This important caveat is mentioned in the documentation, as well as a user-warning when running TORAX with EQDSK geometry.

[jcitrin]

@theo-brown , @fusionby2030

First look has exposed a number of issues:

1. For the eqdsk library to even ingest the file, need the suffix .eqdsk and to remove the footer of CHEASE-added info (useful info, but the eqdsk library is unhappy with it). That's done here

2. In geometry.py

    # TODO(b/375696414) explore replacing with a more robust method
    epsilon = constants.CONSTANTS.eps  # 1e-7

    psi_interpolant = np.linspace(
        eqfile['psimag'] + epsilon, eqfile['psibdry'] - epsilon, 100
    )

    surfaces = []
    cg_psi = contourpy.contour_generator(X, Z, masked_psi_eqdsk_2dgrid)

    for _, _psi in enumerate(psi_interpolant):
      vertices = cg_psi.create_contour(_psi)
      x_surface, z_surface = vertices[0].T[0], vertices[0].T[1]
      surfaces.append((x_surface, z_surface))

Indeed this TODO is important. For this EQDSK, epsilon=1e-7 is too small and no vertices are found (vertices is empty) for the first iteration of for loop, and vertices[0] crashes. I think it may be better to special-case the first "surface" on the magnetic axis with the expected geometry Intermediates for the magnetic axis, and then start the iteration at a position clearly not on-axis.

I can also imagine there being issues with vertices if the number of surfaces requested is much higher than what the eqdsk file resolution can support. Therefore, perhaps the number of surfaces can be a geometry config input option, and can even imagine something fancy like trying to find the vertices of the first non-axis point and if it values, to reduce iteratively the number of requested surfaces until success (and output a user warning).

By hacking epsilon, got this to work. Unfortunately geometry outputs looked weird, i.e. volumes and areas much less than expected for ITER. Probably not a big deal but haven't looked deeper due to lack of time. Will pick this up again later. If someone else would like a crack in the meantime please go ahead, and can push to this PR.

[theo-brown]

For the eqdsk library to even ingest the file, need the suffix .eqdsk and to remove the footer of CHEASE-added info (useful info, but the eqdsk library is unhappy with it). That's done here

IMO this is ok, as it's on the user to provide a valid EQDSK file. Were you flagging this as a problem or as a general thing to note?

[theo-brown]

For this EQDSK, epsilon=1e-7 is too small and no vertices are found (vertices is empty) for the first iteration of for loop, and vertices[0] crashes

Unfortunately geometry outputs looked weird, i.e. volumes and areas much less than expected for ITER.

Yes, got the same error (regarding vertices[0]) and slightly odd geometry outputs on a STEP EQDSK file.

[jcitrin]

IMO this is ok, as it's on the user to provide a valid EQDSK file. Were you flagging this as a problem or as a general thing to note?

I agree. It was just to note.

[jcitrin]

Another small thing I noticed that should be fixed is that Rmaj for the geometry Intermediates is taken as the axis R, as opposed to the geometric R.

[fusionby2030]

I have a pile of TSVV and WPTE tasks to finish before end of year, so it is unclear how quickly I can help out here.

Although, I am interested in this, and will look into generating equilibria from CHEASE and HELENA (both output EQDSK as well) for a size & current scan to investigate the systematic differences. Targeting JET/ASDEX/TCV.

But other commitments mean I can't promise anything quickly on my side.

[jcitrin]

@theo-brown did you take a further look at this? If not I'll take a look tomorrow (Fri)

[jcitrin]

I have a pile of TSVV and WPTE tasks to finish before end of year, so it is unclear how quickly I can help out here.
Although, I am interested in this, and will look into generating equilibria from CHEASE and HELENA
(both output EQDSK as well) for a size & current scan to investigate the systematic differences. Targeting JET/ASDEX/TCV.

Thanks and no problem @fusionby2030 !

[theo-brown]

I haven't had a chance to, sorry, and realistically I won't be able to for another few days. Feel free to have a crack at it in the meantime!

[jcitrin]

No problem. I took a look now.

Major issues

1. Most egregious issue was minor bug in the vpr definition that impacted the plasma volume by several orders of magnitude :)
2. Initial psi is off by quite a bit. Maybe a COCOS thing? Will look deeper

Minor issues

1. Rmaj should be defined based on the LCFS flux surface, not the magnetic axis (due to Shafranov shift)
2. Something seems off with the triangularity comparison with CHEASE, probably a definition issue, will look closer
3. The values for the magnetic axis should be set based on expected values rather than trying to draw a valid contour near the axis as it is now (contour finding may fail causing a crash).
4. The number of contours to find should be a user input, as opposed to hard-coded 100 as now. Should then put a guard where if the first inner contour finding fails, which could happen if the EQDSK file resolution is insufficient compared to the number of surfaces requested, that an informative error message is provided, or alternatively the requested number of surfaces automatically reduced until a contour is found (with a user-warning).
5. Issues arise with the LCFS contour and integrated values for some of the variables, which is expected for diverted geometry. Need to take some special care there.

Will do this and update.

[jcitrin]

Comparison after fixing the volume (vpr) and Rmaj

This is for the test_iterhybrid_predictor_corrector
Initial condition: EQDSK solid, CHEASE dashed

It's quite informative actually.

1. Note the very different psi initial condition, the wrong EQDSK one has a lot more current in the plasma due to this, initially. However, the psi boundary conditions are set the same, so the EQDSK case has this massive initial negative Vloop trying to reduce the plasma current
2. The LCFS q and magnetic shear in the EQDSK case is completely off due to the awry flux surface averaging there, but that doesn't actually impact the internal evolution so just "looks bad" more than anything
   

In spite of this, the two cases end up reasonably close later in the simulation after the current has diffused. Nice to see this robustness. Will carry on fixing now.

[jcitrin]

Result post fixes. Full summary now in PR description