diff --git a/xbout/__init__.py b/xbout/__init__.py index 9511b93d..b3b4029c 100644 --- a/xbout/__init__.py +++ b/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 diff --git a/xbout/cherab/triangulate.py b/xbout/cherab/triangulate.py index a26f73a1..5da9666d 100644 --- a/xbout/cherab/triangulate.py +++ b/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: """ diff --git a/xbout/wall.py b/xbout/wall.py new file mode 100644 index 00000000..f7e9fe25 --- /dev/null +++ b/xbout/wall.py @@ -0,0 +1,116 @@ +""" +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 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 (hasattr(data, "rlim") and hasattr(data, "zlim")): + 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) +