Better nonorthogonal poloidal spacing algorithm

[johnomotani]

The current algorithm for defining poloidal spacing for nonorthogonal grids is not particularly robust, and requires a significant amount of manual tweaking to produce a good grid. An improved method would be very nice!

The current algorithm defines a 'spacing function' for each PsiContour independently. The parameters of the spacing functions vary in a smooth way radially in order to create a sensible grid. This method was chosen in large part because it was convenient to implement.

Now that there are the FineContour objects (which are indepedent of the PsiContour poloidal spacing) providing the underlying skeleton for the grid, and there is the regrid() functionality, it would be possible to implement other algorithms. PsiContour.getRegridded() could easily be modified to accept an array of poloidal positions instead of a spacing function, by modifying this block

hypnotoad/hypnotoad/core/equilibrium.py

Lines 1678 to 1686 in b7ec7d0

	if sfunc is not None:
	s = sfunc(indices)
	# offset fine_contour.interpFunction in case sfunc(0.)!=0.
	sbegin = sfunc(0.0)
	else:
	d = self.get_distance(psi=psi)
	s = (d[self.endInd] - d[self.startInd]) / (npoints - 1) * indices
	sbegin = 0.0

so s is taken from the array passed in instead of calculated from sfunc. Then it would be possible to do some sort of iteration (pseudo-timestepping?) on the array of poloidal positions.

For example, given an array of poloidal positions:

1. call a suitably modified MeshRegion.distributePointsNonorthogonal() to get the R-Z positions
2. work out how 'nonorthogonal' the grid is at each grid point, and which direction it would have to move to make the grid more orthogonal - turn this into an effective 'force' on the grid point
3. define another 'force' that repels each point from other points on the same contour
4. combine the forces to work out an update for the position of each point

...and iterate until a relaxed state is found. To control the grid spacing and deal with non-orthogonal features of the grid (like target plates) it would probably be necessary to have some parameters that adjust the 'forces' depending on position in the grid:

• enhanced or decreased 'spring constant' for the repulsion between points near the targets or X-point, set by nonorthogonal_*_poloidal_spacing_length.
• different strengths for the force that pushes the grid towards orthogonal - this should probably be stronger than the repulsion between grid points when far from X-points or targets (so that a set of grid points is pushed onto a line that follows Grad(psi) perpendicular to flux surfaces, and then that line is moved around as a rigid object relative to the nearby lines), but weaker than the repulsion between grid points when near the targets or X-point, so that the spacing does become non-orthogonal there. Controlled by (some of) the nonorthogonal_*_poloidal_spacing_range* parameters?

If implementing this suggestion, it would probably be a good idea to parallelise MeshRegion.distributePointsNonorthogonal(). The reason this method was not parallelised before was to avoid having to construct sfunc_orthogonal_list on the parallel workers (have to construct on the workers because functions cannot be serialised by dill to be communicated), but the 'orthogonal spacing function' would not be required in the new method.

As another possible speed-up, it should probably be OK to skip the PsiContour.refine() calls until the iteration is nearly converged, as the interpolation from FineContour objects should be pretty accurate.

As an extension of the above suggestion, with a more substantial refactoring of the code, it would be possible to allow the PsiContour objects (for nonorthogonal grids) to extend between (what are currently) different regions - i.e. a single PsiContour object would extend from target to target in the SOL or PFRs, or all the way around the core. Then even the poloidal boundaries between regions could be allowed to move, which might be useful for some configurations (e.g. where the divertor legs are very short so the X-point is close to the wall).

Possibly the region boundaries could be updated in a different way, with the structure kept more like it is at present but with the boundaries being moved. The advantage of the suggestion above (compared to this one) is that the FineContour objects would extend from target to target, etc. and so would not need to be changed when the region boundaries move.

[bendudson]

Hey @johnomotani ! This sounds great; I think the idea of starting with an orthogonal grid, and then deforming it to better match the wall is a good one. I'd started playing around with using JAX (https://github.com/google/jax) for doing this (due to ability to automatically differentiate), but ran into limitations on accuracy of the available interpolation methods. If this deformation could be done inside hypnotoad then that might be an easier path.

It sounds like your proposal restricts points to moving around on a flux surface, so only moving in 1D. That should simplify things considerably, and avoid our mesh becoming too tangled.

Thanks! Happy to help!

[johnomotani]

Hi @bendudson - I have no time to work on this, so if you want to take a shot at it, that would be awesome! I think the sort of 'enabling changes' in hypnotoad should be simple (famous last words!) - not sure how clear I was above... The actual algorithm for deciding how to move the points I've given no thought to!