Non-urgent, lowest priority: Using only the SOLPS geometry, extrapolate flux surface boundaries into the core

[eldond]

This development is only useful for the case where we do not have a valid equilibrium file that covers the flux map into the core. Since we now have that for the primary sample we're working with, this became a lot less important. But it still might be nice someday and we talked about how to do it, so I'm writing down the design ideas so they don't get lost.

1. Start with all the core cells. Completing Populate grid subsets for core, sol, and divertor cells SOLPS2imas.jl#5 first will be very helpful.
2. Use convex hull to identify the cells that are on the outer border of the group of core cells.
   • In IMAS, it's not clear that there's a hard rule for how to enter cells in ggd. In the native SOLPS files like b2fgmtry, it's possible to take a subset of the arrays using information about the cuts. So, there's not an obvious way to select one layer of the mesh at a time and guarantee that it will always work the same way for everything that gets represented in ggd.
3. Record this outer layer of cells and use them to identify the boundary of the outer flux surface. Then remove them from the set being considered, then use convex hull to identify the next layer. Repeat until all the flux surfaces are known.
4. For each flux surface, calculate geometric quantities:
   • Rmax, Rmin, Zmax, Zmin: the extrema of the boundary path
   • Rupper is the R coordinate associated with Zmax, and Rlower is the coordinate of Zmin
   • a = (Rmax - Rmin) / 2: minor radius
   • Rgeo = (Rmax + Rmin) / 2: geometric center
   • elongation is (Zmax - Zmin)/2a
   • Upper and lower triangularity are (Rgeo - Rupper)/a and (Rgeo - Rlower)/a
   • Shafranov shift is Rgeo - Rgeo[LCFS]
5. Look at a real equilibrium file and plot all of these quantities from the magnetic axis to the LCFS. Use the trends to make up reasonable extrapolation rules.
6. Take the geometric quantities from the SOLPS mesh and extrapolate them inwards to get profiles of elongation, kappa, Rshift, aminor, etc. vs. Psin.
7. Calculate the boundary of each of the extrapolated surfaces, such as with the Miller equilibrium definition (which isn't great at the edge but is fine in the core, where we care).
   • Might have to do something fancy about upper/lower triangularity if they're really dissimilar. If they're close, call them the same to make it easier.
8. Think about how to handle the diamagnetic effect of the plasma on the toroidal field.

This should give the flux surfaces needed in order to integrate B to get rho. Then profiles of pressure or whatever vs. rho can be turned into quantities vs R,Z and then input to synthetic diagnostic instrument functions.

Here are trends in some key quantities in a sample equilibrium:

• In the top left, I plotted the toroidal field to compare deviation from the vacuum field.
• Left center worries more about vacuum field vs total field; I'm looking for the diamagnetic effect from the plasma pressure.
• Bottom left has normalized psi vs R for a slice through the midplane, and also shows psi of each flux surface plotted at that surface's geometric center
• Top right is Rgeo vs psi; basically the green line from bottom left but with axes swapped. Easier to read this way.
• Middle right is the Shafranov shift and the minor radius of each surface.
• Bottom right is elongation (kappa) and triangularity (delta).

[eldond]

Here's the code for making this plot:

# One-time prep
OMFIT['eqstuff'] = from_mds_plus(device='DIII-D', shot=180257, times=[3000], get_afile=False)

# Every time you want to repeat the plot
r = OMFIT['eqstuff']['gEQDSK'][3000]['AuxQuantities']['R']
bt = OMFIT['eqstuff']['gEQDSK'][3000]['AuxQuantities']['Bt'][33, :]
OMFIT['eqstuff']['gEQDSK'][3000]['RMAXIS']
bcentr = OMFIT['eqstuff']['gEQDSK'][3000]['BCENTR']
rcentr = OMFIT['eqstuff']['gEQDSK'][3000]['RCENTR']
btvac = OMFIT['eqstuff']['gEQDSK'][3000]['AuxQuantities']['Bt'][-1, :]
btvac2 = bcentr * rcentr / r
psin = OMFIT['eqstuff']['gEQDSK'][3000]['AuxQuantities']['PSIRZ_NORM'][33, :]
middle = psin.argmin()
psinr_left = interp1d(psin[0:middle+1], r[0:middle+1], bounds_error=False)
psinr_right = interp1d(psin[middle:], r[middle:], bounds_error=False)
lcfs_in = psinr_left(1)
lcfs_out = psinr_right(1)
db = bt - btvac2

fig, axs = plt.subplots(3, 2, sharex='col')
axs[-1, 0].set_xlabel('$R$ / m')
axs[-1, 1].set_xlabel('$\psi_N$')

ax = axs[0, 0]
ax.plot(r, bt, label='Full field @ midplane')
ax.plot(r, btvac, label='Vac field = field @ slice w/ no plasma')
ax.plot(r, btvac2, label=r'Vacuum field = $\frac{B_0 R_0}{R}$')
ax.set_ylabel('$B_\phi$ / T')
ax.legend(loc=0)

ax = axs[1, 0]
ax.plot(r, bt-btvac, label='$B_{midplane} - B_{no plasma}$')
ax.plot(r, db, label = r'$B_{midplane} - \frac{R_0 B_0}{R}$')
ax.plot(r, btvac2 - btvac, label=r'$\frac{R_0 B_0}{R} - B_{no plasma}$')
ax.set_ylabel('$\Delta B_\phi$ / T')
ax.legend(loc=0)

ax = axs[2, 0]
ax.plot(r, psin, label='slice through midplane')
for ax in axs[:, 0]:
    for v in [lcfs_in, lcfs_out]:
        ax.axvline(v, ls='--', color='r')
ax.set_ylabel('$\psi_N$')

psin1 = np.linspace(0, 1, 101)
rgeo = np.zeros(len(psin1))
aminor = np.zeros(len(psin1))
for i in range(len(psin1)):
    rl = psinr_left(psin1[i])
    rr = psinr_right(psin1[i])
    rgeo[i] = (rl + rr)/2.0
    aminor[i] = (rr - rl)/2.0

shift = rgeo - rgeo[-1]

axs[2, 0].plot(rgeo, psin1, label='flux surface vs geometric center')
axs[2, 0].legend(loc=0)

ax = axs[0, 1]
ax.set_ylabel('$R_{geo}$ / m')
ax.plot(psin1, rgeo)
ax = axs[1, 1]
ax.plot(psin1, shift, label='$R_{shift}$')
ax.set_ylabel('m')
ax.plot(psin1, aminor, label='$a_{minor}$')
ax.legend(loc=0)
ax = axs[2, 1]
ax.plot(OMFIT['eqstuff']['gEQDSK'][3000]['fluxSurfaces']['geo']['psin'], OMFIT['eqstuff']['gEQDSK'][3000]['fluxSurfaces']['geo']['kap'], label='$\kappa$')
ax.plot(OMFIT['eqstuff']['gEQDSK'][3000]['fluxSurfaces']['geo']['psin'], OMFIT['eqstuff']['gEQDSK'][3000]['fluxSurfaces']['geo']['delta'], label='$\delta$')
ax.legend(loc=0)