Merge pull request #191 from boutproject/fix-xbout-grid-loading

Set the dimension for R_closed_wall and Z_closed_wall to 'closed_wall'

Author: johnomotani

.github/workflows/pythonpackage.yml

@@ -4,7 +4,10 @@
 name: Python package
 
 on:
-  [push, pull_request]
+  push:
+  pull_request:
+  schedule:
+    - cron: '25 1 3 * *'
 
 concurrency:
   group: ${{ github.workflow}}-${{ github.ref }}
@@ -18,7 +21,7 @@ jobs:
     if: always()
     strategy:
       matrix:
-        python-version: ['3.10', '3.11', '3.12']
+        python-version: ['3.10', '3.11', '3.12', '3.13']
       fail-fast: false
 
     steps:
@@ -42,7 +45,7 @@ jobs:
     if: always()
     strategy:
       matrix:
-        python-version: ['3.10', '3.11', '3.12']
+        python-version: ['3.10', '3.11', '3.12', '3.13']
       fail-fast: false
 
     steps:

doc/whats-new.md

@@ -4,19 +4,27 @@ Release history
 ### Bug fixes
 - When using the GUI, if there is an error in `TokamakEquilibrium` object
   creation, still plot the equilibrium data (#186).
+  By [John Omotani](https://github.com/johnomotani)
+- 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)
 
 
 ### New features
 
 - Radial grid line construction can recover from failure and generate
   a rough grid to help visual inspection when option
   `follow_perpendicular_recover` is set to True (#175)
+  By [Ben Dudson](https://github.com/bendudson)
 - A `View` menu enables the grid plot to be customised, with cell edges, corners,
   grid lines and other components (#176).
+  By [Ben Dudson](https://github.com/bendudson)
 - `penalty_mask` is calculated and written to the grid file, based on intersection
   of the grid with the wall. This enables immersed boundary conditions.
+  By [Ben Dudson](https://github.com/bendudson)
 - Wall coordinates are written to output grid as `closed_wall_R` and `closed_wall_Z`
   (#176)
+  By [Ben Dudson](https://github.com/bendudson)
 
 0.5.2 (13th March 2023)
 -------------------------

hypnotoad/core/equilibrium.py

@@ -1854,11 +1854,18 @@ def checkFineContourExtend(self, *, psi):
         minind = numpy.argmin(distances)
         # if minind > 0, or the distance to point 1 is less than the distance between
         # point 0 and point 1 of the fine_contour, then fine_contour extends past p so
-        # does not need to be extended
-        if minind == 0 and distances[1] > numpy.sqrt(
-            numpy.sum(
-                (fine_contour.positions[1, :] - fine_contour.positions[0, :]) ** 2
+        # does not need to be extended.
+        # Include some tolerance to allow for rounding errors when the first point on
+        # the FineContour and the first point on the PsiContour are in 'the same place'.
+        if (
+            minind == 0
+            and distances[1]
+            - numpy.sqrt(
+                numpy.sum(
+                    (fine_contour.positions[1, :] - fine_contour.positions[0, :]) ** 2
+                )
             )
+            > 1.0e-13
         ):
             ds = fine_contour.distance[1] - fine_contour.distance[0]
             n_extend_lower = max(int(numpy.ceil(distances[0] / ds)), 1)
@@ -1874,10 +1881,17 @@ def checkFineContourExtend(self, *, psi):
         # if minind < len(distances)-1, or the distance to the last point is less than
         # the distance between the last and second-last of the fine_contour, then
         # fine_contour extends past p so does not need to be extended
-        if minind == len(distances) - 1 and distances[-2] > numpy.sqrt(
-            numpy.sum(
-                (fine_contour.positions[-1, :] - fine_contour.positions[-2, :]) ** 2
+        # Include some tolerance to allow for rounding errors when the last point on
+        # the FineContour and the last point on the PsiContour are in 'the same place'.
+        if (
+            minind == len(distances) - 1
+            and distances[-2]
+            - numpy.sqrt(
+                numpy.sum(
+                    (fine_contour.positions[-1, :] - fine_contour.positions[-2, :]) ** 2
+                )
             )
+            > 1.0e-13
         ):
             ds = fine_contour.distance[-1] - fine_contour.distance[-2]
             n_extend_upper = max(int(numpy.ceil(distances[-1] / ds)), 1)

hypnotoad/core/mesh.py

@@ -3728,8 +3728,16 @@ def writeGridfile(self, filename):
                 f.write("psi_bdry_gfile", self.equilibrium.psi_bdry_gfile)
 
             if hasattr(self.equilibrium, "closed_wallarray"):
-                f.write("closed_wall_R", self.equilibrium.closed_wallarray[:, 0])
-                f.write("closed_wall_Z", self.equilibrium.closed_wallarray[:, 1])
+                f.write(
+                    "closed_wall_R",
+                    self.equilibrium.closed_wallarray[:, 0],
+                    dims=("closed_wall",),
+                )
+                f.write(
+                    "closed_wall_Z",
+                    self.equilibrium.closed_wallarray[:, 1],
+                    dims=("closed_wall",),
+                )
 
             # write the 2d fields
             for name in self.fields_to_output:

hypnotoad/scripts/hypnotoad_plot_grid_cells.py

@@ -137,25 +137,25 @@ def main():
         # kwarg for something else.
         ds_region.isel(
             x=slice(xin, xout_centre), theta=slice(ylow, yup_centre)
-        ).plot.scatter("R", "Z", marker=".", color="k", s=1)
+        ).plot.scatter(x="R", y="Z", marker=".", color="k", s=1)
 
         # plot radial grid lines
         plt.plot(
-            ds_region["Rxy_corners"].isel(
+            ds_region["Rxy_lower_left_corners"].isel(
                 x=slice(xin, xout_corner_rad), theta=slice(ylow, yup_corner)
             ),
-            ds_region["Zxy_corners"].isel(
+            ds_region["Zxy_lower_left_corners"].isel(
                 x=slice(xin, xout_corner_rad), theta=slice(ylow, yup_corner)
             ),
             color="k",
         )
 
         # plot poloidal grid lines
         plt.plot(
-            ds_region["Rxy_corners"]
+            ds_region["Rxy_lower_left_corners"]
             .isel(x=slice(xin_pol, xout_corner_pol), theta=slice(ylow, yup_corner))
             .T,
-            ds_region["Zxy_corners"]
+            ds_region["Zxy_lower_left_corners"]
             .isel(x=slice(xin_pol, xout_corner_pol), theta=slice(ylow, yup_corner))
             .T,
             color="k",
@@ -175,8 +175,12 @@ def main():
             # neighbouring in the global grid, because if they are the boundary between
             # regions is not (or 'not really') a branch cut.
             plt.plot(
-                ds_region["Rxy_corners"].isel(x=slice(xin, xout_corner_rad), theta=-1),
-                ds_region["Zxy_corners"].isel(x=slice(xin, xout_corner_rad), theta=-1),
+                ds_region["Rxy_lower_left_corners"].isel(
+                    x=slice(xin, xout_corner_rad), theta=-1
+                ),
+                ds_region["Zxy_lower_left_corners"].isel(
+                    x=slice(xin, xout_corner_rad), theta=-1
+                ),
                 color="r",
                 linewidth=3,
                 zorder=1000,
@@ -187,8 +191,12 @@ def main():
             # By arbitrary choice, plot from the region(s) inside the separatrix, so
             # highlight the outer edge.
             plt.plot(
-                ds_region["Rxy_corners"].isel(x=-1, theta=slice(ylow, yup_corner)),
-                ds_region["Zxy_corners"].isel(x=-1, theta=slice(ylow, yup_corner)),
+                ds_region["Rxy_lower_left_corners"].isel(
+                    x=-1, theta=slice(ylow, yup_corner)
+                ),
+                ds_region["Zxy_lower_left_corners"].isel(
+                    x=-1, theta=slice(ylow, yup_corner)
+                ),
                 color="b",
                 linewidth=3,
                 zorder=999,
@@ -197,10 +205,10 @@ def main():
         if targets:
             if ds_region.regions[r].connection_lower_y is None:
                 plt.plot(
-                    ds_region["Rxy_corners"].isel(
+                    ds_region["Rxy_lower_left_corners"].isel(
                         x=slice(xin, xout_corner_rad), theta=ylow
                     ),
-                    ds_region["Zxy_corners"].isel(
+                    ds_region["Zxy_lower_left_corners"].isel(
                         x=slice(xin, xout_corner_rad), theta=ylow
                     ),
                     color="k",
@@ -210,10 +218,10 @@ def main():
             if ds_region.regions[r].connection_upper_y is None:
                 yval = -1 if yup_corner is None else yup_corner
                 plt.plot(
-                    ds_region["Rxy_corners"].isel(
+                    ds_region["Rxy_lower_left_corners"].isel(
                         x=slice(xin, xout_corner_rad), theta=yval
                     ),
-                    ds_region["Zxy_corners"].isel(
+                    ds_region["Zxy_lower_left_corners"].isel(
                         x=slice(xin, xout_corner_rad), theta=yval
                     ),
                     color="k",

integrated_tests/connected_doublenull_with_rounding_error_nonorthogonal/runtest.py

@@ -10,8 +10,8 @@
 
 diagnose = False
 
-rtol = 1.0e-9
-atol = 5.0e-10
+rtol = 2.0e-9
+atol = 2.0e-9
 
 # make sure we are in the test directory
 os.chdir(Path(__file__).parent)

pyproject.toml

@@ -26,7 +26,7 @@ classifiers = [
 
 requires-python = ">=3.10"
 dependencies = [
-    "boututils~=0.1.7",
+    "boutdata~=0.3.0",
     "dill~=0.3,!=0.3.5,!=0.3.5.1",
     "func_timeout~=4.3",
     "matplotlib~=3.7",