Small fixes to enable grids with nonorthogonal X-point but orthogonal targets

doc/nonorthogonal-grid.rst

@@ -29,6 +29,10 @@ detail in the subsections below:
 * Linear spacing: Robust but inflexible, and produces grids with discontinuous
   spacing. Not recommended for simulations, but may be useful for
   exploration/scoping or to provide inputs to other tools.
+* Nonorthogonal X-point: with approximately orthogonal divertor targets. Should
+  result an an almost orthogonal grid, but alleviate the very small poloidal
+  spacing of grid cells on the 'seams' around the X-points that usually results
+  from fully orthogonal grids, and which often limits the timestep size.
 
 Default method
 --------------
@@ -200,3 +204,23 @@ To use set `nonorthogonal_spacing_method = "linear"` and
 
 In each region (target to X-point, or X-point to X-point) distributes the grid
 points poloidally with a uniform spacing in poloidal (not parallel!) distance.
+
+Nonorthogonal X-point
+---------------------
+
+To use set `nonorthogonal_target_all_poloidal_spacing_range = None` and either
+leave `nonorthogonal_target_all_poloidal_spacing_range_inner` and
+`nonorthogonal_target_all_poloidal_spacing_range_outer` unset, or set them to
+`None`.
+
+In most of the grid, uses the orthogonal spacing function so that the grid is
+approximately orthogonal. Near to the X-point uses a fixed perpendicular
+spacing from the region boundaries that start from the X-point and follow the
+directions perpendicular to flux surfaces. This should not make the grid
+strongly nonorthogonal, but reduces the poloidal compression of grid points in
+these areas, which can force timestep reductions (or badly-conditioned matrices
+for implicit solvers, etc.).
+
+Ideally, codes using grids generated with this option should support
+nonorthogonal grids. The sensitivity of results to slight non-orthogonality in
+codes that assume orthogonal grids has not been tested yet...

doc/whats-new.md

@@ -31,6 +31,8 @@ Release history
   By [Ben Dudson](https://github.com/bendudson)
 - Linear poloidal spacing option for nonorthogonal grids (#190).
   By [John Omotani](https://github.com/johnomotani)
+- Grids that are nonorthogonal only at the X-point (#180).
+  By [John Omotani](https://github.com/johnomotani)
 
 0.5.2 (13th March 2023)
 -------------------------

hypnotoad/core/equilibrium.py

@@ -1721,7 +1721,7 @@ def _coarseExtrapUpper(self, reference_ind, *, kind="cubic"):
         )
         return lambda s: Point2D(interpR(s), interpZ(s))
 
-    def contourSfunc(self, *, psi, equilibrium, kind="cubic"):
+    def contourSfunc(self, *, psi, equilibrium, spacings, kind="cubic"):
         """
         Function interpolating distance as a function of index for the current state of
         this contour. When outside [startInd, endInd], set to constant so the results
@@ -1743,12 +1743,44 @@ def contourSfunc(self, *, psi, equilibrium, kind="cubic"):
             thisEndInd += len(self)
         startDistance = distance[thisStartInd]
         endDistance = distance[thisEndInd]
+
+        # When doing nonorthogonal X-point but 'orthogonal' targets, do want to be able
+        # to extrapolate past the end distance, otherwise do not extrapolate.
+        if (
+            (spacings["nonorthogonal_range_lower"] is None)
+            and (spacings["nonorthogonal_range_lower_inner"] is None)
+            and (spacings["nonorthogonal_range_lower_outer"] is None)
+        ):
+
+            def lower_extrap(i):
+                return i * (distance[thisStartInd + 1] - distance[thisStartInd])
+
+        else:
+            lower_extrap = 0.0
+
+        if (
+            (spacings["nonorthogonal_range_upper"] is None)
+            and (spacings["nonorthogonal_range_upper_inner"] is None)
+            and (spacings["nonorthogonal_range_upper_outer"] is None)
+        ):
+
+            def upper_extrap(i):
+                return (
+                    endDistance
+                    - startDistance
+                    + (i - thisEndInd + thisStartInd)
+                    * (distance[thisEndInd] - distance[thisEndInd - 1])
+                )
+
+        else:
+            upper_extrap = endDistance - startDistance
+
         return lambda i: numpy.piecewise(
             i,
             [i <= 0.0, i >= thisEndInd - thisStartInd],
             [
-                0.0,
-                endDistance - startDistance,
+                lower_extrap,
+                upper_extrap,
                 lambda i: interpS(i + thisStartInd) - startDistance,
             ],
         )
@@ -2171,7 +2203,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_positive,
         ),
         nonorthogonal_xpoint_poloidal_spacing_range=WithMeta(
-            lambda options: 0.02 * options.nonorthogonal_xpoint_poloidal_spacing_length,
+            lambda options: (
+                None
+                if options.nonorthogonal_xpoint_poloidal_spacing_length is None
+                else 0.02 * options.nonorthogonal_xpoint_poloidal_spacing_length
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "X-point. This range is used at the radial location of separatrices"
@@ -2180,7 +2216,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_non_negative_or_None,
         ),
         nonorthogonal_xpoint_poloidal_spacing_range_inner=WithMeta(
-            lambda options: 5.0 * options.nonorthogonal_xpoint_poloidal_spacing_range,
+            lambda options: (
+                None
+                if options.nonorthogonal_xpoint_poloidal_spacing_range is None
+                else 5.0 * options.nonorthogonal_xpoint_poloidal_spacing_range
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "X-point. This range is used at 'inner' radial boundaries (core and PFR)"
@@ -2189,7 +2229,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_non_negative_or_None,
         ),
         nonorthogonal_xpoint_poloidal_spacing_range_outer=WithMeta(
-            lambda options: 5.0 * options.nonorthogonal_xpoint_poloidal_spacing_range,
+            lambda options: (
+                None
+                if options.nonorthogonal_xpoint_poloidal_spacing_range is None
+                else 5.0 * options.nonorthogonal_xpoint_poloidal_spacing_range
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "X-point. This range is used at 'outer' radial boundaries (SOL)"
@@ -2209,8 +2253,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_positive,
         ),
         nonorthogonal_target_all_poloidal_spacing_range=WithMeta(
-            lambda options: 0.5
-            * options.nonorthogonal_target_all_poloidal_spacing_length,
+            lambda options: (
+                None
+                if options.nonorthogonal_target_all_poloidal_spacing_length is None
+                else 0.5 * options.nonorthogonal_target_all_poloidal_spacing_length
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "target. This range is used at the radial location of separatrices"
@@ -2219,8 +2266,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_non_negative_or_None,
         ),
         nonorthogonal_target_all_poloidal_spacing_range_inner=WithMeta(
-            lambda options: 2.0
-            * options.nonorthogonal_target_all_poloidal_spacing_range,
+            lambda options: (
+                None
+                if options.nonorthogonal_target_all_poloidal_spacing_range is None
+                else 2.0 * options.nonorthogonal_target_all_poloidal_spacing_range
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "target. This range is used at 'inner' radial boundaries (PFR)"
@@ -2229,8 +2279,11 @@ class EquilibriumRegion(PsiContour):
             check_all=is_non_negative_or_None,
         ),
         nonorthogonal_target_all_poloidal_spacing_range_outer=WithMeta(
-            lambda options: 2.0
-            * options.nonorthogonal_target_all_poloidal_spacing_range,
+            lambda options: (
+                None
+                if options.nonorthogonal_target_all_poloidal_spacing_range is None
+                else 2.0 * options.nonorthogonal_target_all_poloidal_spacing_range
+            ),
             doc=(
                 "Poloidal range over which to use perpendicular spacing near the "
                 "target. This range is used at 'outer' radial boundaries (SOL)"
@@ -2925,9 +2978,7 @@ def getSfuncFixedSpacing(
                     spacing_lower=spacings["nonorthogonal_orthogonal_d_lower"],
                     spacing_upper=spacings["nonorthogonal_orthogonal_d_upper"],
                 )
-            elif (
-                self.nonorthogonal_options.nonorthogonal_spacing_method == "orthogonal"
-            ):
+            elif self.nonorthogonal_options.nonorthogonal_spacing_method == "linear":
                 sfunc = self.getLinearPoloidalDistanceFunc(
                     distance,
                     npoints - 1,
@@ -3250,14 +3301,14 @@ def new_sfunc(i):
                 spacings["nonorthogonal_range_lower_inner"],
                 spacings["nonorthogonal_range_lower"],
                 spacings["nonorthogonal_range_lower_outer"],
-                this_range_lower,
+                # this_range_lower,
             )
             print(
                 "check upper ranges",
                 spacings["nonorthogonal_range_upper_inner"],
                 spacings["nonorthogonal_range_upper"],
                 spacings["nonorthogonal_range_upper_outer"],
-                this_range_upper,
+                # this_range_upper,
             )
             raise
 

hypnotoad/core/mesh.py

@@ -363,10 +363,26 @@ def addPointAtWallToContours(self):
         max_extend = 100
 
         # should the contour intersect a wall at the lower end?
-        lower_wall = self.connections["lower"] is None
+        # if nonorthogonal_target_poloidal_spacing_range_* for this region is `None`,
+        # then ignore the wall intersections.
+        lower_wall = self.connections["lower"] is None and not (
+            not self.user_options.orthogonal
+            and self.equilibriumRegion.getTargetParameter(
+                "nonorthogonal_target_poloidal_spacing_range"
+            )
+            is None
+        )
 
         # should the contour intersect a wall at the upper end?
-        upper_wall = self.connections["upper"] is None
+        # if nonorthogonal_target_poloidal_spacing_range_* for this region is `None`,
+        # then ignore the wall intersections.
+        upper_wall = self.connections["upper"] is None and not (
+            not self.user_options.orthogonal
+            and self.equilibriumRegion.getTargetParameter(
+                "nonorthogonal_target_poloidal_spacing_range"
+            )
+            is None
+        )
 
         # Note: parallel_map will pass additional keywords to the function,
         #       including `equilibrium` and `psi`.
@@ -387,7 +403,9 @@ def addPointAtWallToContours(self):
         # the change in distance after redefining startInd to be at the wall
         self.sfunc_orthogonal_list = [
             contour.contourSfunc(
-                psi=self.equilibriumRegion.psi, equilibrium=self.meshParent.equilibrium
+                psi=self.equilibriumRegion.psi,
+                equilibrium=self.meshParent.equilibrium,
+                spacings=self.equilibriumRegion.getSpacings(),
             )
             for contour in self.contours
         ]

hypnotoad/test_suite/test_equilibrium.py

@@ -604,7 +604,17 @@ def test_contourSfunc(self, testcontour, eq):
         c.extend_lower = 2
         c.extend_upper = 2
 
-        f = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
+        # Put in dummy entry for 'spacings' - doesn't matter what it is as long as it is
+        # not `None`.
+        fake_spacings = {
+            "nonorthogonal_range_lower": 1.0,
+            "nonorthogonal_range_lower_inner": 1.0,
+            "nonorthogonal_range_lower_outer": 1.0,
+            "nonorthogonal_range_upper": 1.0,
+            "nonorthogonal_range_upper_inner": 1.0,
+            "nonorthogonal_range_upper_outer": 1.0,
+        }
+        f = c.contourSfunc(psi=testcontour.psi, equilibrium=eq, spacings=fake_spacings)
 
         indices = numpy.arange(n, dtype=float)
         assert f(indices) == tight_approx(
@@ -659,8 +669,18 @@ def test_contourSfunc_list(self, testcontour, eq):
         # evaluated, it takes 'sfunc_orig' from the enclosing scope, which means it uses
         # the last value from the loop instead of the value when 'sfunc' was added to
         # 'sfunc_list'
+        fake_spacings = {
+            "nonorthogonal_range_lower": 1.0,
+            "nonorthogonal_range_lower_inner": 1.0,
+            "nonorthogonal_range_lower_outer": 1.0,
+            "nonorthogonal_range_upper": 1.0,
+            "nonorthogonal_range_upper_inner": 1.0,
+            "nonorthogonal_range_upper_outer": 1.0,
+        }
         for c in c_list:
-            sfunc_orig = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
+            sfunc_orig = c.contourSfunc(
+                psi=testcontour.psi, equilibrium=eq, spacings=fake_spacings
+            )
 
             sfunc_list.append(lambda i: sfunc_orig(i) - 3.0)
 
@@ -677,7 +697,9 @@ def test_contourSfunc_list(self, testcontour, eq):
         # This version does work, because when the lambda is evaluated it uses
         # 'sfunc_orig' from the scope of 'shift_sfunc' in which it was created.
         def shift_sfunc(c):
-            sfunc_orig = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
+            sfunc_orig = c.contourSfunc(
+                psi=testcontour.psi, equilibrium=eq, spacings=fake_spacings
+            )
             return lambda i: sfunc_orig(i) - 3.0
 
         for c in c_list:
@@ -1422,8 +1444,16 @@ def test_combineSfuncsPerpSpacingIntegrated(self, eqReg):
         eqReg.sin_angle_at_end = 1.0
         eqReg.name = "inner_lower"
 
+        fake_spacings = {
+            "nonorthogonal_range_lower": 1.0,
+            "nonorthogonal_range_lower_inner": 1.0,
+            "nonorthogonal_range_lower_outer": 1.0,
+            "nonorthogonal_range_upper": 1.0,
+            "nonorthogonal_range_upper_inner": 1.0,
+            "nonorthogonal_range_upper_outer": 1.0,
+        }
         sfunc_orthogonal_original = eqReg.contourSfunc(
-            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium, spacings=fake_spacings
         )
 
         # as if lower_wall
@@ -1476,8 +1506,16 @@ def test_combineSfuncsPoloidalSpacingIntegrated(self, eqReg):
         eqReg.ny_total = 40
         eqReg.name = "outer_lower"
 
+        fake_spacings = {
+            "nonorthogonal_range_lower": 1.0,
+            "nonorthogonal_range_lower_inner": 1.0,
+            "nonorthogonal_range_lower_outer": 1.0,
+            "nonorthogonal_range_upper": 1.0,
+            "nonorthogonal_range_upper_inner": 1.0,
+            "nonorthogonal_range_upper_outer": 1.0,
+        }
         sfunc_orthogonal_original = eqReg.contourSfunc(
-            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium, spacings=fake_spacings
         )
 
         # as if lower_wall