From 7fc3e436dabac85a52ef00217cc8465674818f29 Mon Sep 17 00:00:00 2001
From: John Omotani <john.omotani@ukaea.uk>
Date: Mon, 9 Jun 2025 22:53:31 +0100
Subject: [PATCH] Calculate/save parallel distance along flux surfaces

Expected results of tests updated to add the new `parallel_distance`
variables.

Allow tolerances to be increased for certain variables in the integrated
tests when adding random noise, to deal with certain variables seeming
to be more sensitive. For now, applies only to "ShiftTorsion_ylow".
---
 doc/grid-file.rst                             |  12 +
 doc/whats-new.md                              |   1 +
 hypnotoad/core/equilibrium.py                 | 153 +++++++++----
 hypnotoad/core/mesh.py                        | 210 +++++++++++++++---
 hypnotoad/test_suite/test_equilibrium.py      | 171 ++++++++------
 .../expected_nonorthogonal.grd.nc             |   4 +-
 .../expected_orthogonal.grd.nc                |   4 +-
 integrated_tests/utils.py                     |  19 +-
 8 files changed, 427 insertions(+), 147 deletions(-)

diff --git a/doc/grid-file.rst b/doc/grid-file.rst
index a0f38836..5b104039 100644
--- a/doc/grid-file.rst
+++ b/doc/grid-file.rst
@@ -132,6 +132,11 @@ BoutMesh writes three poloidal coordinates to the grid file:
      - The total poloidal distance around a closed flux surface in the core.
        Not calculated on open flux surfaces.
 
+   * - ``total_parallel_distance``
+
+     - The total distance parallel to a magnetic field line around a closed
+       flux surface in the core.  Not calculated on open flux surfaces.
+
    * - ``ShiftAngle``
 
      - The total toroidal angular displacement when following a field line one
@@ -259,6 +264,13 @@ Integral quantities
        surface to the grid point (on open field lines), or from the poloidal
        location of the lower X-point (on closed field lines).
 
+   * - ``parallel_distance``
+
+     - Distance (in metres) parallel to a magnetic field line from the lower
+       divertor target of each flux surface to the grid point (on open field
+       lines), or from the poloidal location of the lower X-point (on closed
+       field lines).
+
    * - ``zShift``
 
      - Toroidal displacement of a field line followed from some reference
diff --git a/doc/whats-new.md b/doc/whats-new.md
index b2e376a4..6f8dbdcb 100644
--- a/doc/whats-new.md
+++ b/doc/whats-new.md
@@ -8,6 +8,7 @@ Release history
 - Set the dimension for R_closed_wall and Z_closed_wall to 'closed_wall'. Fixes
   loading of grid files by xBOUT (#191).
   By [John Omotani](https://github.com/johnomotani)
+- Calculate the parallel distance along field lines, and save to the output (#193).
 
 
 ### New features
diff --git a/hypnotoad/core/equilibrium.py b/hypnotoad/core/equilibrium.py
index d9171da0..f9885f19 100644
--- a/hypnotoad/core/equilibrium.py
+++ b/hypnotoad/core/equilibrium.py
@@ -494,10 +494,11 @@ class FineContour:
         ),
     )
 
-    def __init__(self, parentContour, settings, *, psi):
+    def __init__(self, parentContour, settings, *, psi, equilibrium):
         self.parentContour = parentContour
         self.user_options = self.user_options_factory.create(settings)
         self.distance = None
+        self.parallel_distance = None
         Nfine = self.user_options.finecontour_Nfine
 
         endInd = self.parentContour.endInd
@@ -561,9 +562,9 @@ def __init__(self, parentContour, settings, *, psi):
             self.parentContour.endInd
         ].as_ndarray()
 
-        self.equaliseSpacing(psi=psi)
+        self.equaliseSpacing(psi=psi, equilibrium=equilibrium)
 
-    def extend(self, *, psi, extend_lower=0, extend_upper=0):
+    def extend(self, *, psi, equilibrium, extend_lower=0, extend_upper=0):
         Nfine = self.user_options.finecontour_Nfine
 
         parentCopy = self.parentContour.newContourFromSelf()
@@ -646,9 +647,9 @@ def extend(self, *, psi, extend_lower=0, extend_upper=0):
         self.startInd = self.extend_lower_fine
         self.endInd = Nfine - 1 + self.extend_lower_fine
 
-        self.equaliseSpacing(psi=psi, reallocate=True)
+        self.equaliseSpacing(psi=psi, equilibrium=equilibrium, reallocate=True)
 
-    def equaliseSpacing(self, *, psi, reallocate=False):
+    def equaliseSpacing(self, *, psi, equilibrium, reallocate=False):
         """
         Adjust the positions of points in this :class:`FineContour
         <hypnotoad.core.equilibrium.FineContour>` so they have a constant distance
@@ -686,7 +687,7 @@ def equaliseSpacing(self, *, psi, reallocate=False):
 
         self.refine(psi=psi, skip_endpoints=True)
 
-        self.calcDistance(reallocate=reallocate)
+        self.calcDistance(reallocate=reallocate, equilibrium=equilibrium)
 
         ds = self.distance[1:] - self.distance[:-1]
         # want constant spacing, so ds has a constant value
@@ -773,7 +774,7 @@ def equaliseSpacing(self, *, psi, reallocate=False):
 
             self.refine(psi=psi, skip_endpoints=True)
 
-            self.calcDistance()
+            self.calcDistance(equilibrium=equilibrium)
 
             ds = self.distance[1:] - self.distance[:-1]
             # want constant spacing, so ds has a constant value
@@ -818,7 +819,7 @@ def equaliseSpacing(self, *, psi, reallocate=False):
     def totalDistance(self):
         return self.distance[self.endInd] - self.distance[self.startInd]
 
-    def calcDistance(self, *, reallocate=False):
+    def calcDistance(self, *, equilibrium, reallocate=False):
         """
         Calculate poloidal distance from the start of this :class:`FineContour
         <hypnotoad.core.equilibrium.FineContour>`.
@@ -829,8 +830,35 @@ def calcDistance(self, *, reallocate=False):
         """
         if self.distance is None or reallocate:
             self.distance = numpy.zeros(self.positions.shape[0])
+        if self.parallel_distance is None or reallocate:
+            self.parallel_distance = numpy.zeros(self.positions.shape[0])
         deltaSquared = (self.positions[1:] - self.positions[:-1]) ** 2
-        self.distance[1:] = numpy.cumsum(numpy.sqrt(numpy.sum(deltaSquared, axis=1)))
+        delta_poloidal = numpy.sqrt(numpy.sum(deltaSquared, axis=1))
+        self.distance[1:] = numpy.cumsum(delta_poloidal)
+
+        if self.parentContour.psival is None:
+            # This is an EquilibriumRegion that is not necessarily following a single
+            # flux surface, so does not have a `psival`. We do not need the distance
+            # along this contour, so just set to NaN.
+            self.parallel_distance[:] = numpy.nan
+        else:
+            fine_contour_Bt = (
+                equilibrium.fpol(self.parentContour.psival) / self.positions[:, 0]
+            )
+            fine_contour_Br = equilibrium.Bp_R(
+                self.positions[:, 0], self.positions[:, 1]
+            )
+            fine_contour_Bz = equilibrium.Bp_Z(
+                self.positions[:, 0], self.positions[:, 1]
+            )
+            fine_contour_mod_Bp = numpy.sqrt(fine_contour_Br**2 + fine_contour_Bz**2)
+            fine_contour_B = numpy.sqrt(fine_contour_Bt**2 + fine_contour_mod_Bp**2)
+            midpoints_B = 0.5 * (fine_contour_B[:-1] + fine_contour_B[1:])
+            midpoints_mod_Bp = 0.5 * (
+                fine_contour_mod_Bp[:-1] + fine_contour_mod_Bp[1:]
+            )
+            delta_parallel = midpoints_B / midpoints_mod_Bp * delta_poloidal
+            self.parallel_distance[1:] = numpy.cumsum(delta_parallel)
 
     def interpFunction(self, *, kind="linear"):
         distance = self.distance - self.distance[self.startInd]
@@ -973,17 +1001,18 @@ def interpSSperp(self, vec, kind="linear"):
 
         return s_of_sperp, s_perp_total
 
-    def getDistance(self, p):
+    def getDistance(self, p, *, parallel=False):
         """
-        Find the poloidal distance from the start of this contour of a point ``p``.
+        Find the poloidal and parallel distances from the start of this contour of a
+        point ``p``.
 
         Assume ``p`` is a point on the contour so has the correct psi-value.
 
-        Result is calculated as the weighted mean of the poloidal distances of the two
-        nearest points on the :class:`FineContour
-        <hypnotoad.core.equilibrium.FineContour>` (weighted by the relative distance
-        from ``p`` to each :class:`FineContour <hypnotoad.core.equilibrium.FineContour>`
-        point).
+        Result is calculated as the weighted mean of the poloidal (or parallel)
+        distances of the two nearest points on the :class:`FineContour
+        <hypnotoad.core.equilibrium.FineContour>` (weighted by the relative poloidal (in
+        both cases) distance from ``p`` to each :class:`FineContour
+        <hypnotoad.core.equilibrium.FineContour>` point).
         """
         p = p.as_ndarray()
 
@@ -1012,7 +1041,15 @@ def getDistance(self, p):
         # as their distances from the point
         r = d2 / (d1 + d2)
 
-        return r * self.distance[i1] + (1.0 - r) * self.distance[i2]
+        distance = r * self.distance[i1] + (1.0 - r) * self.distance[i2]
+
+        # Weight by poloidal distances even when calculating parallel distance,
+        # because we cannot get the parallel displacement of `p`, only its position.
+        parallel_distance = (
+            r * self.parallel_distance[i1] + (1.0 - r) * self.parallel_distance[i2]
+        )
+
+        return distance, parallel_distance
 
     def plot(self, *args, psi=None, ax=None, **kwargs):
         """
@@ -1084,6 +1121,7 @@ def __init__(self, *, points, psival, settings, Rrange, Zrange):
         self._fine_contour = None
 
         self._distance = None
+        self._parallel_distance = None
 
         # Value of vector potential on this contour
         self.psival = psival
@@ -1106,6 +1144,7 @@ def _reset_cached(self):
         # Reset all cached objects/values because the contour has been changed
         self._fine_contour = None
         self._distance = None
+        self._parallel_distance = None
 
     @property
     def startInd(self):
@@ -1153,7 +1192,7 @@ def extend_upper(self, val):
             self._reset_cached()
             self._extend_upper = val
 
-    def get_fine_contour(self, *, psi=None):
+    def get_fine_contour(self, *, psi=None, equilibrium=None):
         """
         Get the FineContour associated with this PsiContour
 
@@ -1163,15 +1202,21 @@ def get_fine_contour(self, *, psi=None):
         if self._fine_contour is None:
             if psi is None:
                 raise ValueError("Poloidal flux psi needed to create FineContour")
-            self._fine_contour = FineContour(self, dict(self.user_options), psi=psi)
+            if equilibrium is None:
+                raise ValueError("equilibrium object needed to create FineContour")
+            self._fine_contour = FineContour(
+                self, dict(self.user_options), psi=psi, equilibrium=equilibrium
+            )
             # Ensure that the fine contour is long enough
-            self.checkFineContourExtend(psi=psi)
+            self.checkFineContourExtend(psi=psi, equilibrium=equilibrium)
         return self._fine_contour
 
-    def get_distance(self, *, psi):
+    def get_distance(self, *, psi, equilibrium):
         if self._distance is None:
-            fine_contour = self.get_fine_contour(psi=psi)
-            self._distance = [fine_contour.getDistance(p) for p in self]
+            fine_contour = self.get_fine_contour(psi=psi, equilibrium=equilibrium)
+            both_distances = [fine_contour.getDistance(p) for p in self]
+            self._distance = [d[0] for d in both_distances]
+            self._parallel_distance = [d[1] for d in both_distances]
             d = numpy.array(self._distance)
             if not numpy.all(d[1:] - d[:-1] > 0.0):
                 print("\nPsiContour distance", self._distance)
@@ -1192,6 +1237,11 @@ def get_distance(self, *, psi):
                 )
         return self._distance
 
+    def get_parallel_distance(self, *, psi, equilibrium):
+        if self._parallel_distance is None:
+            self.get_distance(psi=psi, equilibrium=equilibrium)
+        return self._parallel_distance
+
     def __iter__(self):
         return self.points.__iter__()
 
@@ -1312,8 +1362,8 @@ def insertFindPosition(self, point):
             self.insert(minind, point)
             return minind
 
-    def totalDistance(self, *, psi):
-        distance = self.get_distance(psi=psi)
+    def totalDistance(self, *, psi, equilibrium):
+        distance = self.get_distance(psi=psi, equilibrium=equilibrium)
         return distance[self.endInd] - distance[self.startInd]
 
     def reverse(self):
@@ -1600,8 +1650,8 @@ def getRefined(self, skip_endpoints=False, *, width=None, atol=None, **kwargs):
 
         return self.newContourFromSelf(points=newpoints)
 
-    def interpFunction(self, *, psi):
-        return self.get_fine_contour(psi=psi).interpFunction()
+    def interpFunction(self, *, psi, equilibrium):
+        return self.get_fine_contour(psi=psi, equilibrium=equilibrium).interpFunction()
 
     def _coarseInterp(self, *, kind="cubic"):
         distance = [0.0]
@@ -1671,13 +1721,13 @@ def _coarseExtrapUpper(self, reference_ind, *, kind="cubic"):
         )
         return lambda s: Point2D(interpR(s), interpZ(s))
 
-    def contourSfunc(self, *, psi, kind="cubic"):
+    def contourSfunc(self, *, psi, equilibrium, 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
         aren't affected by extrapolation errors.
         """
-        distance = self.get_distance(psi=psi)
+        distance = self.get_distance(psi=psi, equilibrium=equilibrium)
         interpS = interpolate.interp1d(
             numpy.arange(len(self), dtype=float),
             distance,
@@ -1703,7 +1753,7 @@ def contourSfunc(self, *, psi, kind="cubic"):
             ],
         )
 
-    def interpSSperp(self, vec, *, psi):
+    def interpSSperp(self, vec, *, psi, equilibrium):
         """
         Returns
         -------
@@ -1714,7 +1764,7 @@ def interpSSperp(self, vec, *, psi):
            contour.
         2. the total perpendicular distance between startInd and endInd of the contour.
         """
-        return self.get_fine_contour(psi=psi).interpSSperp(vec)
+        return self.get_fine_contour(psi=psi, equilibrium=equilibrium).interpSSperp(vec)
 
     def regrid(self, *args, **kwargs):
         """
@@ -1728,6 +1778,7 @@ def getRegridded(
         npoints,
         *,
         psi,
+        equilibrium,
         width=None,
         atol=None,
         sfunc=None,
@@ -1777,7 +1828,7 @@ def getRegridded(
             # offset fine_contour.interpFunction in case sfunc(0.)!=0.
             sbegin = sfunc(0.0)
         else:
-            d = self.get_distance(psi=psi)
+            d = self.get_distance(psi=psi, equilibrium=equilibrium)
             s = (d[self.endInd] - d[self.startInd]) / (npoints - 1) * indices
             sbegin = 0.0
 
@@ -1787,7 +1838,7 @@ def getRegridded(
         # re-gridding should be the last change that is made to a PsiContour
 
         # Calling get_fine_contour() ensures that self._fine_contour has been initialised
-        self.get_fine_contour(psi=psi)
+        self.get_fine_contour(psi=psi, equilibrium=equilibrium)
 
         orig_extend_lower = self._fine_contour.extend_lower_fine
         orig_extend_upper = self._fine_contour.extend_upper_fine
@@ -1798,7 +1849,9 @@ def getRegridded(
         while (
             s[0] < -self._fine_contour.distance[self._fine_contour.startInd] - tol_lower
         ):
-            self._fine_contour.extend(extend_lower=max(orig_extend_lower, 1), psi=psi)
+            self._fine_contour.extend(
+                extend_lower=max(orig_extend_lower, 1), psi=psi, equilibrium=equilibrium
+            )
 
         tol_upper = 0.25 * (
             self._fine_contour.distance[-1] - self._fine_contour.distance[-2]
@@ -1809,7 +1862,9 @@ def getRegridded(
             - self._fine_contour.distance[self._fine_contour.startInd]
             + tol_upper
         ):
-            self._fine_contour.extend(extend_upper=max(orig_extend_upper, 1), psi=psi)
+            self._fine_contour.extend(
+                extend_upper=max(orig_extend_upper, 1), psi=psi, equilibrium=equilibrium
+            )
 
         interp_unadjusted = self._fine_contour.interpFunction()
 
@@ -1838,13 +1893,13 @@ def interp(s):
 
         return new_contour
 
-    def checkFineContourExtend(self, *, psi):
+    def checkFineContourExtend(self, *, psi, equilibrium):
         """
         Ensure that self._fine_contour extends past the first and last points of this
         PsiContour
         """
 
-        fine_contour = self.get_fine_contour(psi=psi)
+        fine_contour = self.get_fine_contour(psi=psi, equilibrium=equilibrium)
 
         # check first point
         p = numpy.array([*self[0]])
@@ -1902,10 +1957,13 @@ def checkFineContourExtend(self, *, psi):
             return
         else:
             fine_contour.extend(
-                psi=psi, extend_lower=n_extend_lower, extend_upper=n_extend_upper
+                psi=psi,
+                equilibrium=equilibrium,
+                extend_lower=n_extend_lower,
+                extend_upper=n_extend_upper,
             )
             # Call recursively to check extending has gone far enough
-            self.checkFineContourExtend(psi=psi)
+            self.checkFineContourExtend(psi=psi, equilibrium=equilibrium)
 
     def temporaryExtend(
         self, *, psi, extend_lower=0, extend_upper=0, ds_lower=None, ds_upper=None
@@ -2689,7 +2747,7 @@ def nxInsideSeparatrix(self):
     def getRefined(self, *args, **kwargs):
         return self.newRegionFromPsiContour(super().getRefined(*args, **kwargs))
 
-    def getRegridded(self, radialIndex, *, psi, **kwargs):
+    def getRegridded(self, radialIndex, *, psi, equilibrium, **kwargs):
         """ """
         for wrong_argument in ["npoints", "extend_lower", "extend_upper", "sfunc"]:
             # these are valid arguments to PsiContour.getRegridded, but not to
@@ -2708,7 +2766,7 @@ def getRegridded(self, radialIndex, *, psi, **kwargs):
             extend_upper = 2 * self.user_options.y_boundary_guards
         else:
             extend_upper = 0
-        distance = self.get_distance(psi=psi)
+        distance = self.get_distance(psi=psi, equilibrium=equilibrium)
         sfunc = self.getSfuncFixedSpacing(
             2 * self.ny_noguards + 1,
             distance[self.endInd] - distance[self.startInd],
@@ -2717,6 +2775,7 @@ def getRegridded(self, radialIndex, *, psi, **kwargs):
             super().getRegridded(
                 2 * self.ny_noguards + 1,
                 psi=psi,
+                equilibrium=equilibrium,
                 extend_lower=extend_lower,
                 extend_upper=extend_upper,
                 sfunc=sfunc,
@@ -2920,7 +2979,7 @@ def combineSfuncs(
         if vec_lower is None:
             sfunc_fixed_lower = self.getSfuncFixedSpacing(
                 2 * self.ny_noguards + 1,
-                contour.totalDistance(psi=self.psi),
+                contour.totalDistance(psi=self.psi, equilibrium=self.equilibrium),
                 method="monotonic",
                 spacing_lower=spacing_lower,
                 spacing_upper=spacing_upper,
@@ -2938,7 +2997,7 @@ def combineSfuncs(
         if vec_upper is None:
             sfunc_fixed_upper = self.getSfuncFixedSpacing(
                 2 * self.ny_noguards + 1,
-                contour.totalDistance(psi=self.psi),
+                contour.totalDistance(psi=self.psi, equilibrium=self.equilibrium),
                 method="monotonic",
                 spacing_lower=spacing_lower,
                 spacing_upper=spacing_upper,
@@ -3171,7 +3230,9 @@ def new_sfunc(i):
                     (sfunc_fixed_upper, "fixed perp upper"),
                 ],
                 xind=contour.global_xind,
-                total_distance=contour.totalDistance(psi=self.psi),
+                total_distance=contour.totalDistance(
+                    psi=self.psi, equilibrium=self.equilibrium
+                ),
                 prefix="nonorthogonal_",
             )
         except ValueError:
@@ -3226,7 +3287,9 @@ def getSfuncFixedPerpSpacing(
         else:
             d_upper = spacing_upper
 
-        s_of_sperp, s_perp_total = contour.interpSSperp(surface_direction, psi=self.psi)
+        s_of_sperp, s_perp_total = contour.interpSSperp(
+            surface_direction, psi=self.psi, equilibrium=self.equilibrium
+        )
         sperp_func = self.getMonotonicPoloidalDistanceFunc(
             s_perp_total, N - 1, N_norm, d_lower=d_lower, d_upper=d_upper
         )
diff --git a/hypnotoad/core/mesh.py b/hypnotoad/core/mesh.py
index d303a350..b5221f60 100644
--- a/hypnotoad/core/mesh.py
+++ b/hypnotoad/core/mesh.py
@@ -289,7 +289,8 @@ def __init__(
         # Use self.equilibriumRegion.fine_contour for the vector along the separatrix
         # because then the vector will not change when the grid resolution changes
         fine_contour = self.equilibriumRegion.get_fine_contour(
-            psi=self.equilibriumRegion.psi
+            psi=self.equilibriumRegion.psi,
+            equilibrium=self.meshParent.equilibrium,
         )
         if self.equilibriumRegion.wallSurfaceAtStart is None:
             # lower end
@@ -385,7 +386,9 @@ def addPointAtWallToContours(self):
         # point makes the spacing of points on the contour not-smooth) and adjusted for
         # the change in distance after redefining startInd to be at the wall
         self.sfunc_orthogonal_list = [
-            contour.contourSfunc(psi=self.equilibriumRegion.psi)
+            contour.contourSfunc(
+                psi=self.equilibriumRegion.psi, equilibrium=self.meshParent.equilibrium
+            )
             for contour in self.contours
         ]
 
@@ -447,7 +450,8 @@ def correct_sfunc_orthogonal_and_set_startInd(
                     # Calculate total distance using a new FineContour, created at this
                     # point
                     new_total_distance = contour.totalDistance(
-                        psi=self.equilibriumRegion.psi
+                        psi=self.equilibriumRegion.psi,
+                        equilibrium=self.meshParent.equilibrium,
                     )
 
                     return (
@@ -647,6 +651,7 @@ def get_sfunc(i_contour, contour, sfunc_orthogonal):
             c.regrid(
                 2 * self.ny_noguards + 1,
                 psi=self.equilibriumRegion.psi,
+                equilibrium=self.meshParent.equilibrium,
                 sfunc=get_sfunc(i, c, sfunc_orth),
                 extend_lower=self.equilibriumRegion.extend_lower,
                 extend_upper=self.equilibriumRegion.extend_upper,
@@ -856,6 +861,7 @@ def geometry1(self):
         self.Bpxy = numpy.sqrt(self.Brxy**2 + self.Bzxy**2)
 
         self.calcPoloidalDistance()
+        self.calcParallelDistance()
 
         if hasattr(
             self.meshParent.equilibrium.regions[self.equilibriumRegion.name], "pressure"
@@ -1363,13 +1369,19 @@ def calcHy(self):
         for i in range(self.nx):
             print(f"{self.name} calcHy {i} / {2 * self.nx + 1}", end="\r", flush=True)
             d = numpy.array(
-                self.contours[2 * i + 1].get_distance(psi=self.equilibriumRegion.psi)
+                self.contours[2 * i + 1].get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
             )
             hy.centre[i, :] = d[2::2] - d[:-2:2]
             hy.ylow[i, 1:-1] = d[3:-1:2] - d[1:-3:2]
             if self.connections["lower"] is not None:
                 cbelow = self.getNeighbour("lower").contours[2 * i + 1]
-                dbelow = cbelow.get_distance(psi=self.equilibriumRegion.psi)
+                dbelow = cbelow.get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
                 hy.ylow[i, 0] = d[1] - d[0] + dbelow[-1] - dbelow[-2]
             else:
                 # no region below, so estimate distance to point before '0' as the same
@@ -1377,7 +1389,10 @@ def calcHy(self):
                 hy.ylow[i, 0] = 2.0 * (d[1] - d[0])
             if self.connections["upper"] is not None:
                 cabove = self.getNeighbour("upper").contours[2 * i + 1]
-                dabove = cabove.get_distance(psi=self.equilibriumRegion.psi)
+                dabove = cabove.get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
                 hy.ylow[i, -1] = d[-1] - d[-2] + dabove[1] - dabove[0]
             else:
                 # no region below, so estimate distance to point before '0' as the same
@@ -1391,13 +1406,19 @@ def calcHy(self):
                 flush=True,
             )
             d = numpy.array(
-                self.contours[2 * i].get_distance(psi=self.equilibriumRegion.psi)
+                self.contours[2 * i].get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
             )
             hy.xlow[i, :] = d[2::2] - d[:-2:2]
             hy.corners[i, 1:-1] = d[3:-1:2] - d[1:-3:2]
             if self.connections["lower"] is not None:
                 cbelow = self.getNeighbour("lower").contours[2 * i]
-                dbelow = cbelow.get_distance(psi=self.equilibriumRegion.psi)
+                dbelow = cbelow.get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
                 hy.corners[i, 0] = d[1] - d[0] + dbelow[-1] - dbelow[-2]
             else:
                 # no region below, so estimate distance to point before '0' as the same
@@ -1405,7 +1426,10 @@ def calcHy(self):
                 hy.corners[i, 0] = 2.0 * (d[1] - d[0])
             if self.connections["upper"] is not None:
                 cabove = self.getNeighbour("upper").contours[2 * i]
-                dabove = cabove.get_distance(psi=self.equilibriumRegion.psi)
+                dabove = cabove.get_distance(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
                 hy.corners[i, -1] = d[-1] - d[-2] + dabove[1] - dabove[0]
             else:
                 # no region below, so estimate distance to point before '0' as the same
@@ -1548,7 +1572,10 @@ def calcZShift(self):
         while True:
             print("calcZShift", region.name, end="\r", flush=True)
             for i, contour in enumerate(region.contours):
-                fine_contour = contour.get_fine_contour(psi=self.equilibriumRegion.psi)
+                fine_contour = contour.get_fine_contour(
+                    psi=self.equilibriumRegion.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )
                 fine_distance = fine_contour.distance
 
                 def integrand_func(R, Z):
@@ -1581,7 +1608,10 @@ def integrand_func(R, Z):
                     fine_distance, zShift_fine, kind="linear", assume_sorted=True
                 )
                 zShift_contour = zShift_interpolator(
-                    contour.get_distance(psi=self.equilibriumRegion.psi)
+                    contour.get_distance(
+                        psi=self.equilibriumRegion.psi,
+                        equilibrium=self.meshParent.equilibrium,
+                    )
                 )
 
                 if i % 2 == 0:
@@ -1650,21 +1680,25 @@ def calcPoloidalDistance(self):
             c = region.contours[2 * i + 1]
             # Cell-centre points
             region.poloidal_distance.centre[i, :] -= c.get_distance(
-                psi=self.meshParent.equilibrium.psi
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
             )[c.startInd]
             # ylow points
             region.poloidal_distance.ylow[i, :] -= c.get_distance(
-                psi=self.meshParent.equilibrium.psi
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
             )[c.startInd]
         for i in range(self.nx + 1):
             c = region.contours[2 * i]
             # Cell-centre points
             region.poloidal_distance.xlow[i, :] -= c.get_distance(
-                psi=self.meshParent.equilibrium.psi
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
             )[c.startInd]
             # ylow points
             region.poloidal_distance.corners[i, :] -= c.get_distance(
-                psi=self.meshParent.equilibrium.psi
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
             )[c.startInd]
 
         # Get distances from contours
@@ -1673,21 +1707,25 @@ def calcPoloidalDistance(self):
                 c = region.contours[2 * i + 1]
                 # Cell-centre points
                 region.poloidal_distance.centre[i, :] += c.get_distance(
-                    psi=self.meshParent.equilibrium.psi
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
                 )[1::2]
                 # ylow points
                 region.poloidal_distance.ylow[i, :] += c.get_distance(
-                    psi=self.meshParent.equilibrium.psi
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
                 )[::2]
             for i in range(self.nx + 1):
                 c = region.contours[2 * i]
                 # Cell-centre points
                 region.poloidal_distance.xlow[i, :] += c.get_distance(
-                    psi=self.meshParent.equilibrium.psi
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
                 )[1::2]
                 # ylow points
                 region.poloidal_distance.corners[i, :] += c.get_distance(
-                    psi=self.meshParent.equilibrium.psi
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
                 )[::2]
 
             next_region = region.getNeighbour("upper")
@@ -1726,6 +1764,115 @@ def calcPoloidalDistance(self):
                 :, -1
             ]
 
+    def calcParallelDistance(self):
+        """
+        Calculate parallel distance by following contours between regions.
+        """
+        # Cannot just test 'connections['lower'] is not None' because periodic regions
+        # always have a lower connection - requires us to give a yGroupIndex to each
+        # region when creating the groups.
+        if self.yGroupIndex != 0:
+            return None
+
+        region = self
+        region.parallel_distance = MultiLocationArray(region.nx, region.ny)
+        region.parallel_distance.centre = 0.0
+        region.parallel_distance.ylow = 0.0
+        region.parallel_distance.xlow = 0.0
+        region.parallel_distance.corners = 0.0
+
+        # Initialise so that distance counts from the lower wall (for SOL/PFR) or wall
+        # (for core)
+        for i in range(self.nx):
+            c = region.contours[2 * i + 1]
+            # Cell-centre points
+            region.parallel_distance.centre[i, :] -= c.get_parallel_distance(
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
+            )[c.startInd]
+            # ylow points
+            region.parallel_distance.ylow[i, :] -= c.get_parallel_distance(
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
+            )[c.startInd]
+        for i in range(self.nx + 1):
+            c = region.contours[2 * i]
+            # Cell-centre points
+            region.parallel_distance.xlow[i, :] -= c.get_parallel_distance(
+                psi=self.meshParent.equilibrium.psi,
+                equilibrium=self.meshParent.equilibrium,
+            )[c.startInd]
+            # ylow points
+            region.parallel_distance.corners[i, :] -= (
+                c.get_parallel_distance(
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )[c.startInd],
+            )
+
+        # Get distances from contours
+        while True:
+            for i in range(self.nx):
+                c = region.contours[2 * i + 1]
+                # Cell-centre points
+                region.parallel_distance.centre[i, :] += c.get_parallel_distance(
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )[1::2]
+                # ylow points
+                region.parallel_distance.ylow[i, :] += c.get_parallel_distance(
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )[::2]
+            for i in range(self.nx + 1):
+                c = region.contours[2 * i]
+                # Cell-centre points
+                region.parallel_distance.xlow[i, :] += c.get_parallel_distance(
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )[1::2]
+                # ylow points
+                region.parallel_distance.corners[i, :] += c.get_parallel_distance(
+                    psi=self.meshParent.equilibrium.psi,
+                    equilibrium=self.meshParent.equilibrium,
+                )[::2]
+
+            next_region = region.getNeighbour("upper")
+            if (next_region is None) or (next_region is self):
+                # Note: If periodic, next_region is self (back to start)
+                break
+            else:
+                # Initialise with values at the lower y-boundary of next_region
+                next_region.parallel_distance = MultiLocationArray(
+                    next_region.nx, next_region.ny
+                )
+                next_region.parallel_distance.centre[:, :] = (
+                    region.parallel_distance.ylow[:, -1, numpy.newaxis]
+                )
+                next_region.parallel_distance.ylow[:, :] = (
+                    region.parallel_distance.ylow[:, -1, numpy.newaxis]
+                )
+                next_region.parallel_distance.xlow[:, :] = (
+                    region.parallel_distance.corners[:, -1, numpy.newaxis]
+                )
+                next_region.parallel_distance.corners[:, :] = (
+                    region.parallel_distance.corners[:, -1, numpy.newaxis]
+                )
+                region = next_region
+
+        # Save total parallel distance in core
+        self.total_parallel_distance = MultiLocationArray(region.nx, 1)
+        if self.connections["lower"] is not None:
+            # This is a periodic region (we already checked that the self.yGroupIndex is
+            # 0).
+            # 'region' is the last region in the y-group
+            self.total_parallel_distance.centre[:, 0] = region.parallel_distance.ylow[
+                :, -1
+            ]
+            self.total_parallel_distance.xlow[:, 0] = region.parallel_distance.corners[
+                :, -1
+            ]
+
     def getNeighbour(self, face):
         if self.connections[face] is None:
             return None
@@ -2306,13 +2453,13 @@ def smooth_mask(mark):
         setattr(self, varname, tmp)
 
 
-def _calc_contour_distance(i, c, *, psi, **kwargs):
+def _calc_contour_distance(i, c, *, psi, equilibrium, **kwargs):
     print(
         f"Calculating contour distances: {i + 1}",
         end="\r",
         flush=True,
     )
-    c.get_distance(psi=psi)
+    c.get_distance(psi=psi, equilibrium=equilibrium)
     return c
 
 
@@ -2369,7 +2516,7 @@ def _find_intersection(
                 break
 
         count = 0
-        distance = contour.get_distance(psi=psi)
+        distance = contour.get_distance(psi=psi, equilibrium=equilibrium)
         ds_extend = distance[1] - distance[0]
         while coarse_lower_intersect is None:
             # contour has not yet intersected with wall, so make it longer and
@@ -2401,8 +2548,8 @@ def _find_intersection(
         # points
         #
         # first find nearest FineContour points
-        fine_contour = contour.get_fine_contour(psi=psi)
-        d = fine_contour.getDistance(coarse_lower_intersect)
+        fine_contour = contour.get_fine_contour(psi=psi, equilibrium=equilibrium)
+        d, _ = fine_contour.getDistance(coarse_lower_intersect)
         i_fine = numpy.searchsorted(fine_contour.distance, d)
         # Intersection should be between i_fine-1 and i_fine, but check
         # intervals on either side if necessary
@@ -2455,7 +2602,7 @@ def _find_intersection(
                 break
 
         count = 0
-        distance = contour.get_distance(psi=psi)
+        distance = contour.get_distance(psi=psi, equilibrium=equilibrium)
         ds_extend = distance[-1] - distance[-2]
         while coarse_upper_intersect is None:
             # contour has not yet intersected with wall, so make it longer and
@@ -2487,8 +2634,8 @@ def _find_intersection(
         # points
         #
         # first find nearest FineContour points
-        fine_contour = contour.get_fine_contour(psi=psi)
-        d = fine_contour.getDistance(coarse_upper_intersect)
+        fine_contour = contour.get_fine_contour(psi=psi, equilibrium=equilibrium)
+        d, _ = fine_contour.getDistance(coarse_upper_intersect)
         i_fine = numpy.searchsorted(fine_contour.distance, d)
         # Intersection should be between i_fine-1 and i_fine, but check
         # intervals on either side if necessary
@@ -2525,7 +2672,7 @@ def _find_intersection(
         )
 
     # Create FineContour for result, so that this is done in parallel
-    contour.get_fine_contour(psi=psi)
+    contour.get_fine_contour(psi=psi, equilibrium=equilibrium)
 
     return (
         contour,
@@ -2536,7 +2683,7 @@ def _find_intersection(
     )
 
 
-def _refine_extend(i, contour, *, psi, **kwargs):
+def _refine_extend(i, contour, *, psi, equilibrium, **kwargs):
     print(
         "refine and extend",
         i,
@@ -2544,7 +2691,7 @@ def _refine_extend(i, contour, *, psi, **kwargs):
         flush=True,
     )
     contour.refine(psi=psi)
-    contour.checkFineContourExtend(psi=psi)
+    contour.checkFineContourExtend(psi=psi, equilibrium=equilibrium)
     return contour
 
 
@@ -2695,6 +2842,7 @@ def makeRegions(self, parallel_map):
                 eq_region_with_boundaries = eq_region.getRegridded(
                     radialIndex=i,
                     psi=self.equilibrium.psi,
+                    equilibrium=self.equilibrium,
                     width=self.user_options.refine_width,
                 )
                 self.regions[region_id] = MeshRegion(
@@ -3614,6 +3762,8 @@ def addFromRegionsXArray(name):
         addFromRegions("dy")
         addFromRegions("poloidal_distance")
         addFromRegionsXArray("total_poloidal_distance")
+        addFromRegions("parallel_distance")
+        addFromRegionsXArray("total_parallel_distance")
         addFromRegions("Brxy")
         addFromRegions("Bzxy")
         addFromRegions("Bpxy")
diff --git a/hypnotoad/test_suite/test_equilibrium.py b/hypnotoad/test_suite/test_equilibrium.py
index ba8e1bb1..d8de79bb 100644
--- a/hypnotoad/test_suite/test_equilibrium.py
+++ b/hypnotoad/test_suite/test_equilibrium.py
@@ -362,6 +362,24 @@ def test_closest_approach_right():
     assert numpy.isclose(cpa, 0.5)
 
 
+class ThisEquilibrium(Equilibrium):
+    def __init__(self, settings=None, wall=None):
+        if settings is None:
+            settings = {}
+        self.user_options = Equilibrium.user_options_factory.add(
+            refine_width=1.0e-5, refine_atol=2.0e-8
+        ).create(settings)
+
+        if wall is not None:
+            self.wall = wall
+
+        self.fpol = lambda psi: 1.0 + 0.0 * psi
+        self.Bp_R = lambda R, Z: 1.0 + 0.0 * R
+        self.Bp_Z = lambda R, Z: 1.0 + 0.0 * R
+
+        super().__init__({})
+
+
 class TestContour:
     @pytest.fixture
     def testcontour(self):
@@ -394,9 +412,22 @@ def psifunc(R, Z):
 
         return returnObject()
 
-    def test_distance(self, testcontour):
+    @pytest.fixture
+    def eq(self):
+        wall = [
+            Point2D(-1.0, -1.0),
+            Point2D(1.0, -1.0),
+            Point2D(1.0, 1.0),
+            Point2D(-1.0, 1.0),
+        ]
+        eq = ThisEquilibrium(wall=wall)
+        return eq
+
+    def test_distance(self, testcontour, eq):
         segment_length = testcontour.r * numpy.pi / (testcontour.npoints - 1)
-        assert testcontour.c.get_distance(psi=testcontour.psi) == pytest.approx(
+        assert testcontour.c.get_distance(
+            psi=testcontour.psi, equilibrium=eq
+        ) == pytest.approx(
             segment_length * numpy.arange(testcontour.npoints), abs=1.0e-4
         )
 
@@ -412,29 +443,29 @@ def test_getitem(self, testcontour):
         assert p.R == tight_approx(testcontour.R[5])
         assert p.Z == tight_approx(testcontour.Z[5])
 
-    def test_append(self, testcontour):
+    def test_append(self, testcontour, eq):
         c = testcontour.c
         point_to_add = c[-2]
-        expected_distance = c.get_distance(psi=testcontour.psi)[-2]
+        expected_distance = c.get_distance(psi=testcontour.psi, equilibrium=eq)[-2]
         del c.points[-1]
         del c.points[-1]
         c.endInd = -1
         c.append(point_to_add)
-        assert c.get_distance(psi=testcontour.psi)[-1] == pytest.approx(
+        assert c.get_distance(psi=testcontour.psi, equilibrium=eq)[-1] == pytest.approx(
             expected_distance, abs=4.0e-5
         )
 
-    def test_prepend(self, testcontour):
+    def test_prepend(self, testcontour, eq):
         c = testcontour.c
         c.startInd = 3
         c.endInd = -1
         point_to_add = c[1]
-        expected_distance = c.get_distance(psi=testcontour.psi)[1]
+        expected_distance = c.get_distance(psi=testcontour.psi, equilibrium=eq)[1]
         del c.points[0]
         del c.points[0]
         c.startInd = 1
         c.prepend(point_to_add)
-        assert c.get_distance(psi=testcontour.psi)[1] == pytest.approx(
+        assert c.get_distance(psi=testcontour.psi, equilibrium=eq)[1] == pytest.approx(
             expected_distance, abs=6.0e-4
         )
         assert c.startInd == 1
@@ -455,20 +486,20 @@ def test_replace(self, testcontour):
         assert c[5] == point_to_add
         assert c.endInd == 22
 
-    def test_reverse(self, testcontour):
+    def test_reverse(self, testcontour, eq):
         c = testcontour.c
         orig = deepcopy(c)
 
         c.reverse()
 
         n = len(orig)
-        total_d = orig.get_distance(psi=testcontour.psi)[-1]
+        total_d = orig.get_distance(psi=testcontour.psi, equilibrium=eq)[-1]
         for i in range(n):
             assert orig[n - 1 - i].R == tight_approx(c[i].R)
             assert orig[n - 1 - i].Z == tight_approx(c[i].Z)
-            assert total_d - orig.get_distance(psi=testcontour.psi)[
+            assert total_d - orig.get_distance(psi=testcontour.psi, equilibrium=eq)[
                 n - 1 - i
-            ] == tight_approx(c.get_distance(psi=testcontour.psi)[i])
+            ] == tight_approx(c.get_distance(psi=testcontour.psi, equilibrium=eq)[i])
 
     def test_insertFindPosition(self, testcontour):
         c = testcontour.c
@@ -500,8 +531,8 @@ def test_coarseInterp(self, testcontour):
         assert pend.R == pytest.approx(testcontour.R0 - testcontour.r, abs=1.0e-9)
         assert pend.Z == pytest.approx(testcontour.Z0, abs=1.0e-5)
 
-    def test_interpFunction(self, testcontour):
-        f = testcontour.c.interpFunction(psi=testcontour.psi)
+    def test_interpFunction(self, testcontour, eq):
+        f = testcontour.c.interpFunction(psi=testcontour.psi, equilibrium=eq)
         pstart = f(0.0)
         pend = f(numpy.pi * testcontour.r)
         assert pstart.R == pytest.approx(testcontour.R0 + testcontour.r, abs=1.0e-9)
@@ -509,7 +540,7 @@ def test_interpFunction(self, testcontour):
         assert pend.R == pytest.approx(testcontour.R0 - testcontour.r, abs=1.0e-8)
         assert pend.Z == pytest.approx(testcontour.Z0, abs=1.0e-5)
 
-    def test_getRegridded(self, testcontour):
+    def test_getRegridded(self, testcontour, eq):
         orig = testcontour.c
         r = testcontour.r
 
@@ -521,7 +552,7 @@ def sfunc_true(i):
         def sfunc(i):
             return (
                 numpy.sqrt(i / (newNpoints - 1))
-                * orig.get_distance(psi=testcontour.psi)[-1]
+                * orig.get_distance(psi=testcontour.psi, equilibrium=eq)[-1]
             )
 
         newTheta = sfunc_true(numpy.arange(newNpoints)) / r
@@ -529,19 +560,20 @@ def sfunc(i):
         newZ = testcontour.Z0 + r * numpy.sin(newTheta)
 
         new = orig.getRegridded(
-            newNpoints, psi=testcontour.psi, sfunc=sfunc, width=1.0e-3
+            newNpoints, psi=testcontour.psi, equilibrium=eq, sfunc=sfunc, width=1.0e-3
         )
 
         assert [p.R for p in new] == pytest.approx(newR, abs=4.0e-4)
         assert [p.Z for p in new] == pytest.approx(newZ, abs=4.0e-4)
 
-    def test_getRegridded_extend(self, testcontour):
+    def test_getRegridded_extend(self, testcontour, eq):
         c = testcontour.c
         orig = c.newContourFromSelf()
 
         new = c.getRegridded(
             testcontour.npoints,
             psi=testcontour.psi,
+            equilibrium=eq,
             width=0.1,
             extend_lower=1,
             extend_upper=2,
@@ -564,7 +596,7 @@ def test_getRegridded_extend(self, testcontour):
             [orig[-3].R, 2.0 * testcontour.Z0 - orig[-3].Z], abs=1.0e-7
         )
 
-    def test_contourSfunc(self, testcontour):
+    def test_contourSfunc(self, testcontour, eq):
         c = testcontour.c
         c.startInd = 2
         c.endInd = len(c) - 2
@@ -572,22 +604,24 @@ def test_contourSfunc(self, testcontour):
         c.extend_lower = 2
         c.extend_upper = 2
 
-        f = c.contourSfunc(psi=testcontour.psi)
+        f = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
 
         indices = numpy.arange(n, dtype=float)
         assert f(indices) == tight_approx(
             [
-                d - c.get_distance(psi=testcontour.psi)[c.startInd]
-                for d in c.get_distance(psi=testcontour.psi)[c.startInd : c.endInd + 1]
+                d - c.get_distance(psi=testcontour.psi, equilibrium=eq)[c.startInd]
+                for d in c.get_distance(psi=testcontour.psi, equilibrium=eq)[
+                    c.startInd : c.endInd + 1
+                ]
             ]
         )
         assert f(-1.0) == tight_approx(0.0)
         assert f(n + 1.0) == tight_approx(
-            c.get_distance(psi=testcontour.psi)[c.endInd]
-            - c.get_distance(psi=testcontour.psi)[c.startInd]
+            c.get_distance(psi=testcontour.psi, equilibrium=eq)[c.endInd]
+            - c.get_distance(psi=testcontour.psi, equilibrium=eq)[c.startInd]
         )
 
-    def test_contourSfunc_list(self, testcontour):
+    def test_contourSfunc_list(self, testcontour, eq):
         c1 = testcontour.c
         c1.startInd = 2
         c1.endInd = len(c1) - 2
@@ -626,7 +660,7 @@ def test_contourSfunc_list(self, testcontour):
         # the last value from the loop instead of the value when 'sfunc' was added to
         # 'sfunc_list'
         for c in c_list:
-            sfunc_orig = c.contourSfunc(psi=testcontour.psi)
+            sfunc_orig = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
 
             sfunc_list.append(lambda i: sfunc_orig(i) - 3.0)
 
@@ -643,7 +677,7 @@ def test_contourSfunc_list(self, testcontour):
         # 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)
+            sfunc_orig = c.contourSfunc(psi=testcontour.psi, equilibrium=eq)
             return lambda i: sfunc_orig(i) - 3.0
 
         for c in c_list:
@@ -656,24 +690,29 @@ def shift_sfunc(c):
             n2 / (npoints2 - 1.0) * 1.5 * numpy.pi * r - 3.0, abs=4.0e-6
         )
 
-    def test_interpSSperp(self, testcontour):
+    def test_interpSSperp(self, testcontour, eq):
         c = testcontour.c
 
         # Make c.startInd > 0
         c.insert(0, Point2D(c[1].R, 2.0 * c[0].Z - c[1].Z))
 
         # 'vec' argument is in Z-direction, so 's_perp' is displacement in R-direction
-        sfunc, s_perp_total = c.interpSSperp([0.0, 1.0], psi=testcontour.psi)
+        sfunc, s_perp_total = c.interpSSperp(
+            [0.0, 1.0], psi=testcontour.psi, equilibrium=eq
+        )
         assert sfunc(0.0) == tight_approx(0.0)
         assert sfunc(1.0) == pytest.approx(numpy.pi / 2.0, abs=1.0e-6)
         assert sfunc(2.0) == pytest.approx(numpy.pi, abs=2.0e-6)
         assert s_perp_total == tight_approx(2.0)
 
-    def test_FineContour(self, testcontour):
+    def test_FineContour(self, testcontour, eq):
         testcontour.c.refine_width = 1.0e-2
 
         fc = FineContour(
-            testcontour.c, dict(testcontour.c.user_options), psi=testcontour.psi
+            testcontour.c,
+            dict(testcontour.c.user_options),
+            psi=testcontour.psi,
+            equilibrium=eq,
         )
 
         assert fc.totalDistance() == pytest.approx(numpy.pi, abs=1.0e-5)
@@ -689,7 +728,7 @@ def test_FineContour(self, testcontour):
                 testcontour.Z0 + r * numpy.sin(theta), abs=1.0e-4
             )
 
-    def test_finecontour_extent_lower(self, testcontour):
+    def test_finecontour_extent_lower(self, testcontour, eq):
         contour = testcontour.c
 
         contour_initial = contour[0]
@@ -698,15 +737,15 @@ def test_finecontour_extent_lower(self, testcontour):
 
         # This will create contour._fine_contour and check that contour._distance
         # is monotonic at this point
-        contour.get_distance(psi=testcontour.psi)
+        contour.get_distance(psi=testcontour.psi, equilibrium=eq)
 
         contour.prepend(contour_initial)
 
         # Check that after modifying contour, contour._fine_contour still extends far
         # enough to give a monotonic contour._distance
-        contour.get_distance(psi=testcontour.psi)
+        contour.get_distance(psi=testcontour.psi, equilibrium=eq)
 
-    def test_finecontour_extent_upper(self, testcontour):
+    def test_finecontour_extent_upper(self, testcontour, eq):
         contour = testcontour.c
 
         contour_final = contour[-1]
@@ -715,26 +754,12 @@ def test_finecontour_extent_upper(self, testcontour):
 
         # This will create contour._fine_contour and check that contour._distance
         # is monotonic at this point
-        contour.get_distance(psi=testcontour.psi)
+        contour.get_distance(psi=testcontour.psi, equilibrium=eq)
 
         contour.append(contour_final)
         # Check that after modifying contour, contour._fine_contour still extends far
         # enough to give a monotonic contour._distance
-        contour.get_distance(psi=testcontour.psi)
-
-
-class ThisEquilibrium(Equilibrium):
-    def __init__(self, settings=None, wall=None):
-        if settings is None:
-            settings = {}
-        self.user_options = Equilibrium.user_options_factory.add(
-            refine_width=1.0e-5, refine_atol=2.0e-8
-        ).create(settings)
-
-        if wall is not None:
-            self.wall = wall
-
-        super().__init__({})
+        contour.get_distance(psi=testcontour.psi, equilibrium=eq)
 
 
 class TestEquilibrium:
@@ -1281,7 +1306,7 @@ def test_getSqrtPoloidalDistanceFuncBothBoth(self, eqReg):
 
     def test_combineSfuncsPoloidalSpacing(self, eqReg):
         n = len(eqReg)
-        L = eqReg.totalDistance(psi=eqReg.psi)
+        L = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         eqReg.resetNonorthogonalOptions(
             {"nonorthogonal_target_all_poloidal_spacing_length": L}
@@ -1300,7 +1325,7 @@ def sfunc_orthogonal(i):
 
     def test_combineSfuncsPoloidalSpacingRangeLower(self, eqReg):
         n = len(eqReg)
-        L = eqReg.totalDistance(psi=eqReg.psi)
+        L = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         new_settings = {
             "nonorthogonal_target_all_poloidal_spacing_length": L * 40.0 / (n - 1),
@@ -1327,7 +1352,7 @@ def sfunc_orthogonal(i):
 
     def test_combineSfuncsPoloidalSpacingRangeUpper(self, eqReg):
         n = len(eqReg)
-        L = eqReg.totalDistance(psi=eqReg.psi)
+        L = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         new_settings = {
             "nonorthogonal_target_all_poloidal_spacing_length": L * 40.0 / (n - 1),
@@ -1354,7 +1379,7 @@ def sfunc_orthogonal(i):
 
     def test_combineSfuncsPoloidalSpacingRangeBoth(self, eqReg):
         n = len(eqReg)
-        L = eqReg.totalDistance(psi=eqReg.psi)
+        L = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         new_settings = {
             "nonorthogonal_target_all_poloidal_spacing_length": L * 40.0 / (n - 1),
@@ -1397,12 +1422,16 @@ def test_combineSfuncsPerpSpacingIntegrated(self, eqReg):
         eqReg.sin_angle_at_end = 1.0
         eqReg.name = "inner_lower"
 
-        sfunc_orthogonal_original = eqReg.contourSfunc(psi=eqReg.psi)
+        sfunc_orthogonal_original = eqReg.contourSfunc(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )
 
         # as if lower_wall
         intersect_index = 2
         original_start = eqReg.startInd
-        distance_at_original_start = eqReg.get_distance(psi=eqReg.psi)[original_start]
+        distance_at_original_start = eqReg.get_distance(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )[original_start]
 
         d = 1.5 * 3.0 / (n - 1.0)
         eqReg.points.insert(intersect_index, Point2D(d, d))
@@ -1410,7 +1439,9 @@ def test_combineSfuncsPerpSpacingIntegrated(self, eqReg):
         if original_start >= intersect_index:
             original_start += 1
 
-        distance_at_wall = eqReg.get_distance(psi=eqReg.psi)[intersect_index]
+        distance_at_wall = eqReg.get_distance(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )[intersect_index]
 
         sfunc_orthogonal = (
             lambda i: sfunc_orthogonal_original(i)
@@ -1422,12 +1453,14 @@ def test_combineSfuncsPerpSpacingIntegrated(self, eqReg):
         eqReg.extend_lower = intersect_index
         eqReg.extend_upper = 1
 
-        d = eqReg.totalDistance(psi=eqReg.psi)
+        d = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         sfunc = eqReg.combineSfuncs(eqReg, sfunc_orthogonal, [0.0, 1.0], [1.0, 0.0])
 
         assert sfunc(0.0) == tight_approx(0.0)
-        assert sfunc(n - 1.0) == tight_approx(eqReg.totalDistance(psi=eqReg.psi))
+        assert sfunc(n - 1.0) == tight_approx(
+            eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
+        )
 
     def test_combineSfuncsPoloidalSpacingIntegrated(self, eqReg):
         # This test follows roughly the operations in
@@ -1443,7 +1476,9 @@ def test_combineSfuncsPoloidalSpacingIntegrated(self, eqReg):
         eqReg.ny_total = 40
         eqReg.name = "outer_lower"
 
-        sfunc_orthogonal_original = eqReg.contourSfunc(psi=eqReg.psi)
+        sfunc_orthogonal_original = eqReg.contourSfunc(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )
 
         # as if lower_wall
         intersect_index = 2
@@ -1454,9 +1489,13 @@ def test_combineSfuncsPoloidalSpacingIntegrated(self, eqReg):
 
         if original_start >= intersect_index:
             original_start += 1
-        distance_at_original_start = eqReg.get_distance(psi=eqReg.psi)[original_start]
+        distance_at_original_start = eqReg.get_distance(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )[original_start]
 
-        distance_at_wall = eqReg.get_distance(psi=eqReg.psi)[intersect_index]
+        distance_at_wall = eqReg.get_distance(
+            psi=eqReg.psi, equilibrium=eqReg.equilibrium
+        )[intersect_index]
 
         sfunc_orthogonal = (
             lambda i: sfunc_orthogonal_original(i)
@@ -1468,9 +1507,11 @@ def test_combineSfuncsPoloidalSpacingIntegrated(self, eqReg):
         eqReg.extend_lower = intersect_index
         eqReg.extend_upper = 1
 
-        d = eqReg.totalDistance(psi=eqReg.psi)
+        d = eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
 
         sfunc = eqReg.combineSfuncs(eqReg, sfunc_orthogonal)
 
         assert sfunc(0.0) == tight_approx(0.0)
-        assert sfunc(n - 1.0) == tight_approx(eqReg.totalDistance(psi=eqReg.psi))
+        assert sfunc(n - 1.0) == tight_approx(
+            eqReg.totalDistance(psi=eqReg.psi, equilibrium=eqReg.equilibrium)
+        )
diff --git a/integrated_tests/connected_doublenull_nonorthogonal/expected_nonorthogonal.grd.nc b/integrated_tests/connected_doublenull_nonorthogonal/expected_nonorthogonal.grd.nc
index a3db571b..a672f29a 100644
--- a/integrated_tests/connected_doublenull_nonorthogonal/expected_nonorthogonal.grd.nc
+++ b/integrated_tests/connected_doublenull_nonorthogonal/expected_nonorthogonal.grd.nc
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:7dfb841cbbfbc89ab458c1e09e7108572521f40392bbb87164de419eb4353403
-size 638387
+oid sha256:b57f8e20d7c8155dc8d8dcf86d3613b1f247363e0a02b414f78a89a4d2df3b08
+size 671498
diff --git a/integrated_tests/connected_doublenull_orthogonal/expected_orthogonal.grd.nc b/integrated_tests/connected_doublenull_orthogonal/expected_orthogonal.grd.nc
index 8f4e63bb..e862bd74 100644
--- a/integrated_tests/connected_doublenull_orthogonal/expected_orthogonal.grd.nc
+++ b/integrated_tests/connected_doublenull_orthogonal/expected_orthogonal.grd.nc
@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:1155390ac1606a4f89597e4a6699948b9d3aba819ebb4d597f897e71880e0d44
-size 656877
+oid sha256:e1bdde054651788e9977dd51b49a700a9d18e2b0d55a89e5d60a26159f5a501a
+size 687632
diff --git a/integrated_tests/utils.py b/integrated_tests/utils.py
index 7fa2ce95..7c3f3619 100644
--- a/integrated_tests/utils.py
+++ b/integrated_tests/utils.py
@@ -19,9 +19,12 @@
     "module_versions",
     # Variables that have been added. These entries can be removed if/when the
     # expected output is re-generated.
-    "penalty_mask",
-    "closed_wall_R",
-    "closed_wall_Z",
+]
+
+# Variables that are extra sensitive to rounding errors. Tolerances will be increased by
+# 1e2 for these variables.
+expected_large_error_vars = [
+    "ShiftTorsion_ylow",
 ]
 
 
@@ -111,11 +114,21 @@ def run_case(name, inputfile, expectedfile, *, rtol, atol, diagnose, add_noise=N
         .load()
         .drop(expected_different_vars, errors="ignore")
     )
+    if add_noise is not None:
+        # Increase tolerances for some variables that are sensitive to rounding errors.
+        expected_large_errors = expected[expected_large_error_vars]
+        expected = expected.drop(expected_large_error_vars)
+        actual_large_errors = actual[expected_large_error_vars]
+        actual = actual.drop(expected_large_error_vars)
 
     if diagnose:
         check_errors(expected, actual, rtol=rtol, atol=atol)
 
     xrt.assert_allclose(expected, actual, rtol=rtol, atol=atol)
+    if add_noise is not None:
+        xrt.assert_allclose(
+            expected_large_errors, actual_large_errors, rtol=1e2 * rtol, atol=1e2 * atol
+        )
 
     for attrname in expected_different_attrs:
         del expected.attrs[attrname]