Add a View menu, improve grid plotting

doc/grid-file.rst

@@ -138,6 +138,27 @@ BoutMesh writes three poloidal coordinates to the grid file:
        full poloidal turn around a closed flux surface. Not calculated on open
        flux surfaces.
 
+Wall location
++++++++++++++
+
+If a wall is defined in the equilibrium then the coordinates of the
+closed wall is saved:
+
+.. list-table::
+   :widths: 30 70
+
+   * - ``closed_wall_R``
+
+     - Major radius locations [in meters] of points on the wall. This
+       array forms a closed loop so the last element is the same as
+       the first.
+
+   * - ``closed_wall_Z``
+
+     - Height locations [in meters] of points on the wall, the same
+       number of points as ``closed_wall_R``. This array forms a
+       closed loop so the last element is the same as the first.
+
 2D arrays
 ---------
 
@@ -211,6 +232,21 @@ Magnetic field quantities
 
      - Total magnetic field.
 
+Boundary quantities
++++++++++++++++++++
+
+.. list-table::
+   :widths: 30 70
+
+   * - ``penalty_mask``
+
+     - A 2D mask indicating whether a cell is inside or outside the
+       wall. It's value is 1 for cells entirely outside the wall; 0
+       for cells entirely inside the wall. Cells that cross the wall
+       are given a penalty proportional to the fraction of the cell
+       poloidal length that is inside the wall.
+
+
 Integral quantities
 +++++++++++++++++++
 

doc/whats-new.md

@@ -1,6 +1,21 @@
 Release history
 ===============
 
+### Bug fixes
+
+
+### 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)
+- A `View` menu enables the grid plot to be customised, with cell edges, corners,
+  grid lines and other components (#176).
+- `penalty_mask` is calculated and written to the grid file, based on intersection
+  of the grid with the wall. This enables immersed boundary conditions.
+- Wall coordinates are written to output grid as `closed_wall_R` and `closed_wall_Z`
+  (#176)
+
 0.5.2 (13th March 2023)
 -------------------------
 

hypnotoad/cases/circular.py

@@ -293,9 +293,7 @@ def d2psidr2_r(self, r):
             def func(x):
                 return (
                     B0 / (np.sqrt(1.0 - x**2 / R0**2) * self.q(x))
-                    + B0
-                    * x**2
-                    / (R0**2 * (1.0 - x**2 / R0**2) ** 1.5 * self.q(x))
+                    + B0 * x**2 / (R0**2 * (1.0 - x**2 / R0**2) ** 1.5 * self.q(x))
                     - B0
                     * x
                     * self.dqdr(x)
@@ -326,10 +324,7 @@ def psi_r(self, r):
 
                 def func(x):
                     return (
-                        B0
-                        * R0**2
-                        / coef_array[0]
-                        * (1.0 - np.sqrt(1.0 - x**2 / R0**2))
+                        B0 * R0**2 / coef_array[0] * (1.0 - np.sqrt(1.0 - x**2 / R0**2))
                     )
 
                 self._psi_r = func
@@ -346,8 +341,7 @@ def func(x):
                                 * (
                                     -1.0
                                     + np.sqrt(
-                                        (a1 * (-(x**2) + R0**2))
-                                        / (a0 + a1 * R0**2)
+                                        (a1 * (-(x**2) + R0**2)) / (a0 + a1 * R0**2)
                                     )
                                 )
                             )
@@ -356,8 +350,7 @@ def func(x):
                                 * (
                                     1.0
                                     + np.sqrt(
-                                        (a1 * (-(x**2) + R0**2))
-                                        / (a0 + a1 * R0**2)
+                                        (a1 * (-(x**2) + R0**2)) / (a0 + a1 * R0**2)
                                     )
                                 )
                             )

hypnotoad/core/equilibrium.py

@@ -4792,10 +4792,7 @@ def psi(i):
                             * numpy.pi
                             * CosInt_m_j1
                             * numpy.sin(b * numpy.pi / n)
-                            + b**2
-                            * numpy.pi
-                            * numpy.cos(b * numpy.pi / n)
-                            * SinInt_b_n
+                            + b**2 * numpy.pi * numpy.cos(b * numpy.pi / n) * SinInt_b_n
                             - b
                             * j1(i)
                             * numpy.pi
@@ -4812,14 +4809,8 @@ def psi(i):
                             * numpy.pi
                             * numpy.cos(b * numpy.pi / n)
                             * SinInt_b_n
-                            + n**2
-                            * numpy.pi
-                            * numpy.cos(b * numpy.pi / n)
-                            * SinInt_b_n
-                            - b**2
-                            * numpy.pi
-                            * numpy.cos(b * numpy.pi / n)
-                            * SinInt_j1
+                            + n**2 * numpy.pi * numpy.cos(b * numpy.pi / n) * SinInt_b_n
+                            - b**2 * numpy.pi * numpy.cos(b * numpy.pi / n) * SinInt_j1
                             + b
                             * j1(i)
                             * numpy.pi
@@ -4836,10 +4827,7 @@ def psi(i):
                             * numpy.pi
                             * numpy.cos(b * numpy.pi / n)
                             * SinInt_j1
-                            - n**2
-                            * numpy.pi
-                            * numpy.cos(b * numpy.pi / n)
-                            * SinInt_j1
+                            - n**2 * numpy.pi * numpy.cos(b * numpy.pi / n) * SinInt_j1
                         )
                         + grad_upper / 2 * j2(i) * (b + n) / (b + j2(i))
                         - grad_upper