Allow for profiles created with map3d to extend outside the last closed flux surface

[TKhare]

I am currently using CHERAB to create a forward model for a diagnostic that views the edge region of a tokamak, outside the last closed flux surface. I have been using map3d() to map one dimensional profiles that extend past the lcfs in 3d space. However, I recently noticed that the map3d() function has a value_outside_lcfs parameter, and prevents the user from defining non-constant values for profiles outside the lcfs.

I've gotten around this with a simple modification to the source code in efit.pyx, changing the map2d() function in line 252 from return ScalarBlend2D(value_outside_lcfs, f, self.inside_lcfs) to return ScalarBlend2D(value_outside_lcfs, f, 1) so that values outside the lcfs are not masked off. Could the map3d function be changed altogether to be more general and allow for profiles outside of psin=1? Perhaps taking a psin cutoff parameter for profiles that extend to a certain psin value?

[vsnever]

Hi @TKhare,

I think that removing the inside_lcfs mask will cause trouble in the area below the X-point, where ψ_norm is also < 1, but where the 1D profiles are definitely not applicable. The intended way to build a core + edge plasma is to blend the 1D core profiles with 2D edge profiles simulated with a dedicated edge code.

However, the map functions can be generalised by allowing the user to provide a custom mask:

map2d(self, object profile, Function2D mask=None, double value_outside_mask=0.0)

If no mask is specified, the inside_lcfs mask is used.
We'll have to extrapolate the profile, to ensure that the resulting function does not throw an error within the domain where the equilibrium is defined:

    profile = Interpolator1DArray(profile[0, :], profile[1, :], 'cubic', ''quadratic'', INFINITY)

I would like to hear others' opinions. Is it worth generalising the map functions, or, given that they are quite simple and not cythonised, users can define their own mapping functions if necessary?

[Mateasek]

Hi @TKhare and @vsnever

A while back I was also proposing to extend the interpolation ouside LCFS, but the decision was not to take this step (sorry, I can't find the issue). Although there can be some asymmetries within LCFS, the assumption of constant values on surfaces is a good approximation and is being used quite often. This does not apply for SOL, where the quantities are far from being constant on surfaces. I personally think that we should not allow the SOL mapping within the EFITEquilibrium class. The Raysect function framework should be capable of covering everything needed to make some kind of edge profiles, shouldn't it? You can use blending to distinguish between in/out of lcfs (as is done in the map2d), it should be possible to use blending and clamp on psi normallised to distinguish between PFR and SOL, shouldn't it? This would still not add the evolution along surfaces and should be used in limited cases. Maybe we could add a demo showing some tricks which can be used to get some arbitrary edge profiles instead?

[jacklovell]

As above, it doesn't make sense to map SOL/PFR quantities to the equilibrium: the relavant quantities are not flux functions outside the LCFS. If you don't have a sensible 2D or 3D profile from a code output to build an interpolator from, you could consider parametrising a SOL quantity as $f(\psi, L_\parallel$) which would be more realistic than $f(\psi)$, but would still need separate treatment in the PFR and some sort of contour-following algorithm to calculate $L_\parallel$. I'm not aware of any quantities in the EFIT output which enable differentiating between SOL and PFR, so you'd likely have to make your own too or find an external tool.

From memory, the main reason we didn't add mapping outside the LCFS in Cherab core is because there isn't a generic, machine-independent way to do it well (and core should be general and machine independent).

[rmchurch]

There are situations where specifying an $f(\psi)$ even in the SOL could be useful (e.g. general testing, or also our specific case where the diagnostic line-of-sight doesn't vary much in Z, and so doesn't need $f(\psi,L_\parallel)$). For right now, we can simply reimplement map2d. Its probably good not to have extrapolation past the LCFS as default behavior, but it doesn't seem to be a completely useless ability to have.
I think you're right by default EFIT doesn't have quantities to identify PFR, would have to use something like OMFIT's X-point finder, and then check psin and the Z value compared to the Zx https://omfit.io/_modules/omfit_classes/omfit_eqdsk.html

[jacklovell]

For the use case of wanting to evaluate a function $f(\psi(r, z))$ at an arbitrary point in space, with no special treatment of the separatrix as a boundary between closed and open field lines, the IsoMapper2D class is more appropriate than EFITEquilibrium.map2d. Indeed, the major significant difference between IsoMapper2D and EFITEquilibrium.map2d is that the latter is specifically designed not to extend the profile outside the separatrix because it's so complicated to do properly and generically.

map2d is really designed with core profiles in mind, though I don't think this is properly explained in the documentation. We could look at providing a separate function or method for mapping edge profiles, but this will need some careful thought given to its design as it would need a different API to the core method.

We should definitely add a demo illustrating the simple (and common) use cases @rmchurch describes, where we make clear what simplifications are being made and what caveats about the overall validity of the model apply.