Nonorthogonal grid generation refactor proposal

[johnomotani]

Current nonorthogonal poloidal spacing algorithm

The current grid generation process for nonorthogonal grids (https://hypnotoad.readthedocs.io/en/latest/nonorthogonal-grid.html), in particular defining the poloidal positions of the points along each flux surface, is not very robust. It is based on defining (analytically) a 'spacing function' $s(i/N)$ that maps grid index $i$ (normalised by total number of poloidal grid points $N$) to poloidal distance along a flux surface contour. Constructing an analytic function that conforms to the vessel wall but also tends towards an orthogonal grid far from the wall or X-points is difficult, and requires rather too many parameters to get even as far as it has in the current implementation. The one advantage of the current implementation is that a series of grids can be generated for a scan in poloidal resolution that are consistent with each other because they use the same spacing functions - for example that if the resolution is doubled, each grid cell is divided exactly into two with the original cell faces remaining in the same positions. For this reason it is probably useful to keep the original algorithm as an option, but the default for most use cases should be a more robust algorithm such as the one described below.

Proposed new algorithm

The code generating sfuncs for poloidal spacing is a spaghetti-fied mess - I would suggest not touching it or trying to understand it. Instead, add a new algorithm that can be called by PsiContour.getRegridded() here

hypnotoad/hypnotoad/core/equilibrium.py

Lines 1774 to 1782 in 0365640

	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

with a new option added to control which algorithm to use.

Algorithm:

1. As an initial guess, place the required number of grid points by some simple, robust algorithm (e.g. uniform poloidal spacing, as is currently done if no sfunc is passed to PsiContour.getRegridded()).

2. Think of the grid points as beads on wires, with the beads connected by springs, where the 'wires' are the flux surfaces (represented by PsiContour/FineContour objects), and the springs are various forcing terms added to create the desired grid:

   1. Calculate how 'orthogonal' the grid is at the current iteration, by measuring the deviation of the lines connecting one grid point to its inner and outer radial neighbours from the normal to the flux surface (the vector [eq.f_R(p.R,p.Z), eq.f_Z(p.R,p.Z)], where p is the Point2D corresponding to the grid point and eq is the Equilibrium object). Push the grid point in the direction that reduces the average of the two deviations from orthogonal, by an amount proportional to the average deviation.
   2. Through most of the grid, all other things being equal, we would like a constant poloidal spacing of grid points. Making the grid approximately orthogonal should (probably?) be dominant, but a sub-dominant term that repels grid points from each other would make the radial grid lines repel each other and arrange themselves reasonably nicely. So for each of the nearest grid points on either side, have a force inversely proportional to the distance to that grid point (or maybe the square of the distance, by analogy with the Coulomb force?) and directed away from it.
   3. We want the option to have (poloidally) compressed grids near the targets, with a spacing that has a fixed value at the target and some kind of range that blends into the spacing far from the target. Suggestion: say the user-requested grid spacing at the targets is d_t. For the i'th grid point away from the target, have a force proportional to (s - i*d_t). This force should be large enough to be dominant for i=1, and then drop off (probably with a user-definable range r_t). I'm not sure whether it's better that the range depends on i or on the poloidal distance from the target - I would guess distance from the target. Maybe the prefactor could be something like (r_t / |s - s_t|)**p_t where s_t is the position of the target, and p_t is some user-definable power. This function diverges as |s - s_t| goes to zero, so that if r_t > d_t the force will be very large at the first grid point, and p_t can be tuned to set how quickly the target-compression switches off.
   4. Probably want a similar option to iii. to control spacing around the X-points. This could be skipped for the initial implementation though, and added as a later stage, as it should be possible to generate grids without it. In this case it might be useful to allow the user to choose whether the spacing should be fixed in poloidal distance, or in parallel distance along the field line. To do spacing in parallel distance, probably need to write an interpolation function that calculates parallel distance (which will be calculated on the FineContour objects once Calculate/save parallel distance along flux surfaces #193 is merged) from the poloidal distance.
   5. Another possibly useful extra would be a minimum grid spacing, which might help prevent grids getting too compressed when there is a very inclined or funny-shaped target. Could take the form of an extra repulsive force, similar to ii. but proportional to a higher power, e.g. 1 / separation**8.

   The actual grid would be given by the solution where all these forces are in equilibrium.

3. To find the equilibrium, likely a variety of approaches are possible. I'm not familiar with this kind of problem, so do not know what the standard or optimal solution would be. One possibility that seems likely to be reasonably simple and reasonably robust would be a crude time-stepping algorithm. Define a (poloidal) velocity v for each grid point. Start them all at rest, then apply an acceleration with the net force from 2. Add an extra damping term proportional to -v for stability. So the ODEs to be solved are
   $$\dot{s} = v$$
   $$\dot{v} = F - D v$$
   where $F$ is the force from 2. and $D$ is some damping coefficient. Advance, using forward-Euler timestepping (for simplicity) until all the magnitudes of all the velocities are below some tolerance. When updating grid point positions

   hypnotoad/hypnotoad/core/equilibrium.py

   Lines 835 to 852 in 0365640

   	def interpFunction(self, *, kind="linear"):
   	distance = self.distance - self.distance[self.startInd]
   	interpR = interpolate.interp1d(
   	distance,
   	self.positions[:, 0],
   	kind=kind,
   	assume_sorted=True,
   	fill_value="extrapolate",
   	)
   	interpZ = interpolate.interp1d(
   	distance,
   	self.positions[:, 1],
   	kind=kind,
   	assume_sorted=True,
   	fill_value="extrapolate",
   	)
   	return lambda s: Point2D(float(interpR(s)), float(interpZ(s)))

   will be useful to convert poloidal location s into (R,Z) positions.

4. For performance, I would suggest not using c.refine() where c is a PsiContour object during the timestepping. The pattern in the existing code is to do a refine() after just about every operation that moves the grid points - partially this as a hang-over from before the FineContour objects were introduced, which made the interpolated positions a lot more accurate, and partially it is just done for accuracy because the cost was not too high. For this algorithm, I think it would be sufficient to call refine() once, after the time-loop has converged.

Performance?

The current structure of PsiContour objects as essentially lists of Point2D objects is probably very far from optimal in terms of performance. For this algorithm you'd prefer a numpy array of all the points in the MeshRegion that lumped together a bunch of different PsiContours. Performance during grid generation is not the biggest issue though, so I would suggest not thinking about refactoring, etc. unless the performance is prohibitive. I've wondered if it might be nice to create temporary full-MeshRegion numpy arrays of some quantities just for this time-loop, but I think it might then be difficult/impossible to use the existing functions from PsiContour and FineContour that we'd want for interpolation, etc. so not sure this is a viable option.

Other notes

If there are a lot of errors when generating the initial set of PsiContour objects, before starting step 1., it might be worth implementing the suggestion in #197. I don't expect that to be a problem though, so wouldn't bother with #197 unless it turns out to be one.