cherab interface: Calculate heat fluxes to walls

docs/cherab.rst

@@ -0,0 +1,59 @@
+Cherab interface
+================
+
+Cherab (https://www.cherab.info/) is a python library for forward
+modelling diagnostics based on spectroscopic plasma emission.  It is
+based on the Raysect (http://www.raysect.org/) scientific ray-tracing
+framework.
+
+Triangulation
+-------------
+
+Before Cherab can be used, a triangulated mesh must be generated. This
+mesh is stored in an attribute ``cherab_grid``. This can be attached to
+a Dataset or DataArray by calling the ``bout.with_cherab_grid()`` accessor::
+
+  ds = ds.bout.with_cherab_grid()
+
+Following operations will generate a grid if needed, but if this
+attribute is already present then the grid is not
+recalculated. Calling this method on a dataset before performing
+Cherab operations improves efficiency.
+
+
+Wall heat fluxes
+----------------
+
+To calculate radiation heat fluxes to walls, first create an ``xbout.AxisymmetricWall``
+object that represents a 2D (R, Z) axisymmetric wall. This can be done by
+reading the wall coordinates from a GEQDSK equilibrium file::
+
+   wall = xbout.wall.read_geqdsk("geqdsk")
+
+Triangulate a grid and attach to the dataset::
+
+  ds = ds.bout.with_cherab_grid()
+
+Extract a data array. This must be 2D (x, y) so select a single time
+slice and toroidal angle e.g. excitation radiation::
+
+  da = -bd['Rd+_ex'].isel(t=1, zeta=0)
+
+Remove potentially invalid data in the guard cells::
+
+  da[:2,:] = 0.0
+  da[-2:,:] = 0.0
+
+Attach data to the Cherab triangulated mesh, returning an
+``xbout.cherab.TriangularData`` object that Cherab can work with::
+
+  data = da.bout.as_cherab_data()
+
+Calculate the wall heat fluxes::
+
+  result = data.wall_flux(cat_wall, pixel_samples=1000)
+
+This result is a list of dicts, one for each wall element. Those dicts
+contain coordinates "rz1" and "rz2", results "power_density" and
+"total_power".
+

docs/index.rst

@@ -21,6 +21,7 @@ suggestions.
 
    loading_data
    accessor_methods
+   cherab
    extending_xbout
 
 .. toctree::

xbout/__init__.py

@@ -1,6 +1,6 @@
 from .load import open_boutdataset, collect
 
-from . import geometries
+from . import geometries, wall
 from .geometries import register_geometry, REGISTERED_GEOMETRIES
 
 from .boutdataset import BoutDatasetAccessor

xbout/boutdataarray.py

@@ -1165,7 +1165,7 @@ def as_cherab_emitter(
             parent=parent,
             cylinder_zmin=cylinder_zmin,
             cylinder_zmax=cylinder_zmax,
-            cylinder_rmin=cyliner_rmin,
-            cylinder_rmax=cyliner_rmax,
+            cylinder_rmin=cylinder_rmin,
+            cylinder_rmax=cylinder_rmax,
             step=step,
         )

xbout/cherab/triangulate.py

@@ -7,7 +7,6 @@
 """
 
 import numpy as np
-import xarray as xr
 
 
 class TriangularData:
@@ -74,7 +73,7 @@ def to_emitter(
             material=emitting_material,
         )
 
-    def plot_2d(self, ax=None, nr: int = 150, nz: int = 150):
+    def plot_2d(self, ax=None, nr: int = 150, nz: int = 150, colorbar: bool = True):
         """
         Make a 2D plot of the data
 
@@ -95,12 +94,136 @@ def plot_2d(self, ax=None, nr: int = 150, nz: int = 150):
         image = ax.imshow(
             emiss_sampled.T, origin="lower", extent=(r.min(), r.max(), z.min(), z.max())
         )
-        fig.colorbar(image)
+        if colorbar:
+            fig.colorbar(image)
         ax.set_xlabel("r")
         ax.set_ylabel("z")
 
         return ax
 
+    def wall_flux(
+        self,
+        wall,
+        pixel_samples: int = 2000,
+        wall_detector_offset: float = 0.001,
+        toroidal_width: float = 0.01,
+        step: float = 0.01,
+    ):
+        """
+        Calculate power onto segments of an axisymmetric wall
+
+        Based on the Cherab manual here:
+        https://www.cherab.info/demonstrations/radiation_loads/symmetric_power_load.html#symmetric-power-load
+
+        Parameters
+        ----------
+
+        wall : AxisymmetricWall
+            Defines a closed loop in R-Z, that defines an axisymmetric wall.
+
+        pixel_samples : The number of samples per wall segment
+
+        wall_detector_offset
+            Distance of detector pixels from wall [meters].
+            Slightly displaced to avoid numerical overlap (ray trapping)
+
+        toroidal_width
+            toroidal width of the detectors [meters]
+
+        step
+            Ray path step length [meters]
+
+        """
+        from raysect.core import Point3D, Vector3D, translate, rotate_basis
+        from raysect.optical import World
+        from raysect.optical.observer import PowerPipeline0D
+        from raysect.optical.material import AbsorbingSurface
+        from raysect.optical.observer.nonimaging.pixel import Pixel
+        from cherab.tools.primitives import axisymmetric_mesh_from_polygon
+
+        world = World()
+
+        # Create a plasma emitter
+        emitter = self.to_emitter(parent=world, step=step)
+
+        # Create the walls around the plasma
+        wall_polygon = wall.to_polygon()  # [npoints, 2] array
+
+        # create a 3D mesh from the 2D polygon outline using symmetry
+        wall_mesh = axisymmetric_mesh_from_polygon(wall_polygon)
+        wall_mesh.parent = world
+        wall_mesh.material = AbsorbingSurface()
+
+        result = []
+        for rz1, rz2 in wall:
+            p1 = Point3D(rz1[0], 0, rz1[1])
+            p2 = Point3D(rz2[0], 0, rz2[1])
+
+            # evaluate y_vector
+            y_vector_full = p1.vector_to(p2)
+
+            if y_vector_full.length < 1e-5:
+                continue  # Skip; points p1, p2 are the same
+
+            y_vector = y_vector_full.normalise()
+            y_width = y_vector_full.length
+
+            # evaluate normal_vector (inward pointing)
+            normal_vector = y_vector.cross(Vector3D(0, 1, 0)).normalise()
+
+            # evaluate the central point of the detector
+            # Note: Displaced in the normal direction by wall_detector_offset
+            detector_center = (
+                p1 + y_vector_full * 0.5 + wall_detector_offset * normal_vector
+            )
+
+            # Use the power pipeline to record total power arriving at the surface
+            power_data = PowerPipeline0D()
+
+            # Note: Displace the pixel location
+            pixel_transform = translate(
+                detector_center.x, detector_center.y, detector_center.z
+            ) * rotate_basis(normal_vector, y_vector)
+
+            # Use pixel_samples argument to increase amount of sampling and reduce noise
+            pixel = Pixel(
+                [power_data],
+                x_width=toroidal_width,
+                y_width=y_width,
+                name=f"wall-{rz1}-{rz2}",
+                spectral_bins=1,
+                transform=pixel_transform,
+                parent=world,
+                pixel_samples=pixel_samples,
+            )
+
+            pixel_area = toroidal_width * y_width
+
+            # Start collecting samples
+            pixel.observe()
+
+            # Estimate power density
+            power_density = power_data.value.mean / pixel_area
+
+            # For checking energy conservation.
+            # Revolve this tile around the CYLINDRICAL z-axis to get total power collected by these tiles.
+            # Add up all the tile contributions to get total power collected.
+            detector_radius = np.sqrt(detector_center.x**2 + detector_center.y**2)
+            observed_total_power = power_density * (
+                y_width * 2 * np.pi * detector_radius
+            )
+
+            result.append(
+                {
+                    "rz1": rz1,
+                    "rz2": rz2,
+                    "power_density": power_density,
+                    "total_power": observed_total_power,
+                }
+            )
+
+        return result
+
 
 class Triangulate:
     """

xbout/wall.py

@@ -0,0 +1,135 @@
+"""
+Routines to read and represent wall geometries
+"""
+
+import numpy as np
+
+
+class AxisymmetricWall:
+    def __init__(self, Rs, Zs):
+        """
+        Defines a 2D (R,Z) axisymmetric wall
+
+        Parameters
+        ----------
+
+        Rs : list or 1D array
+            Major radius coordinates [meters]
+        Zs : list or 1D array
+            Vertical coordinates [meters]
+
+        """
+        if len(Rs) != len(Zs):
+            raise ValueError("Rs and Zs arrays have different lengths")
+
+        # Ensure that the members are numpy arrays
+        self.Rs = np.array(Rs)
+        self.Zs = np.array(Zs)
+
+    def __iter__(self):
+        """
+        Iterate over wall elements
+
+        Each iteration returns a pair of (R,Z) pairs:
+            ((R1, Z1), (R2, Z2))
+
+        These pairs define wall segment.
+        """
+        return iter(
+            zip(zip(self.Rs, self.Zs), zip(np.roll(self.Rs, -1), np.roll(self.Zs, -1)))
+        )
+
+    def refine(self, factor: int):
+        """
+        Split wall elements into factor elements
+
+        Returns a new AxisymmetricWall
+        """
+
+        # Wall is a closed loop
+        xold = np.linspace(0, 1, len(self.Rs), endpoint=False)
+        xnew = np.linspace(0, 1, len(self.Rs) * factor, endpoint=False)
+
+        Rs = np.interp(xnew, xold, self.Rs, period=1.0)
+        Zs = np.interp(xnew, xold, self.Zs, period=1.0)
+
+        return AxisymmetricWall(Rs, Zs)
+
+    def to_polygon(self):
+        """
+        Returns a 2D Numpy array [npoints, 2]
+        Index 0 is major radius (R) in meters
+        Index 1 is height (Z) in meters
+        """
+        return np.stack((self.Rs, self.Zs), axis=-1)
+
+    def plot(self, linestyle="k-", ax=None):
+        """
+        Plot the wall on given axis. If no axis
+        is given then a new figure is created.
+
+        Returns
+        -------
+
+        The matplotlib axis containing the plot
+        """
+
+        import matplotlib.pyplot as plt
+
+        if ax is None:
+            fig, ax = plt.subplots()
+
+        ax.plot(self.Rs, self.Zs, linestyle)
+        return ax
+
+
+def read_geqdsk(filehandle):
+    """
+    Read wall geometry from a GEQDSK file.
+
+    Note: Requires the freeqdsk package
+    """
+
+    if isinstance(filehandle, str):
+        with open(filehandle, "r") as f:
+            return read_geqdsk(f)
+
+    from freeqdsk import geqdsk
+
+    data = geqdsk.read(filehandle)
+    # rlim and zlim should be 1D arrays of wall coordinates
+
+    if not (("rlim" in data) and ("zlim" in data)):
+        raise ValueError(f"Wall coordinates not found in GEQDSK file")
+
+    return AxisymmetricWall(data["rlim"], data["zlim"])
+
+
+def read_csv(filehandle, delimiter=","):
+    """
+    Parameters
+    ----------
+
+    filehandle: File handle
+        Must contain two columns, for R and Z coordinates [meters]
+
+    delimier : character
+        A single character that separates fields
+
+    Notes:
+    - Uses the python `csv` module
+    """
+    import csv
+
+    reader = csv.reader(filehandle, delimiter=delimiter)
+    Rs = []
+    Zs = []
+    for row in reader:
+        if len(row) == 0:
+            continue  # Skip empty rows
+        if len(row) != 2:
+            raise ValueError(f"CSV row should contain two columns: {row}")
+        Rs.append(float(row[0]))
+        Zs.append(float(row[1]))
+
+    return AxisymmetricWall(Rs, Zs)