From 40863c7a330c91e61f0a8cf7c4816598c5ea4b98 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Mon, 4 Mar 2024 10:58:54 +0000 Subject: [PATCH 1/9] Improve regularisation documentation and add ADMT demo * Add a Regularisation section to the tomography documentation, describing the functionality inside the admt_utils module. * Fixup the docstrings in admt_utils. * Add calculation of "skewed" second derivative operators, which operate along the diagonals of the grid. * Add a bolometer diagonstic system to Generomak. * Add a new example using the Generomak bolometers and the regularisation operators to perform isotropic and ADMT inversions. --- cherab/generomak/diagnostics/__init__.py | 1 + cherab/generomak/diagnostics/bolometers.py | 167 ++++++++++ cherab/tools/inversions/admt_utils.py | 98 ++++-- .../bolometry/admt_tomographic_inversion.py | 285 ++++++++++++++++++ docs/source/tools/tomography.rst | 33 ++ 5 files changed, 556 insertions(+), 28 deletions(-) create mode 100644 cherab/generomak/diagnostics/__init__.py create mode 100644 cherab/generomak/diagnostics/bolometers.py create mode 100644 demos/observers/bolometry/admt_tomographic_inversion.py diff --git a/cherab/generomak/diagnostics/__init__.py b/cherab/generomak/diagnostics/__init__.py new file mode 100644 index 00000000..ee4a0fc8 --- /dev/null +++ b/cherab/generomak/diagnostics/__init__.py @@ -0,0 +1 @@ +from .bolometers import load_bolometers diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py new file mode 100644 index 00000000..e4693712 --- /dev/null +++ b/cherab/generomak/diagnostics/bolometers.py @@ -0,0 +1,167 @@ +""" +Some foil bolometers for measuring total radiated power. +""" +from raysect.core import (Node, Point3D, Vector3D, rotate_basis, + rotate_x, rotate_y, rotate_z, translate) +from raysect.core.math import rotate +from raysect.optical.material import AbsorbingSurface +from raysect.primitive import Box, Subtract + +from cherab.tools.observers import BolometerCamera, BolometerSlit, BolometerFoil + + +# Convenient constants +XAXIS = Vector3D(1, 0, 0) +YAXIS = Vector3D(0, 1, 0) +ZAXIS = Vector3D(0, 0, 1) +ORIGIN = Point3D(0, 0, 0) +# Bolometer geometry, independent of camera. +BOX_WIDTH = 0.05 +BOX_WIDTH = 0.2 +BOX_HEIGHT = 0.07 +BOX_DEPTH = 0.2 +SLIT_WIDTH = 0.004 +SLIT_HEIGHT = 0.005 +FOIL_WIDTH = 0.0013 +FOIL_HEIGHT = 0.0038 +FOIL_CORNER_CURVATURE = 0.0005 +FOIL_SEPARATION = 0.00508 # 0.2 inch between foils + + +def _make_bolometer_camera(slit_sensor_separation, sensor_angles): + """ + Build a single bolometer camera. + + The camera consists of a box with a rectangular slit and 4 sensors, + each of which has 4 foils. + In its local coordinate system, the camera's slit is located at the + origin and the sensors below the z=0 plane, looking up towards the slit. + """ + camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), + upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) + # Hollow out the box + inside_box = Box(lower=camera_box.lower + Vector3D(1e-5, 1e-5, 1e-5), + upper=camera_box.upper - Vector3D(1e-5, 1e-5, 1e-5)) + camera_box = Subtract(camera_box, inside_box) + # The slit is a hole in the box + aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1e-4), + upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 1e-4)) + camera_box = Subtract(camera_box, aperture) + camera_box.material = AbsorbingSurface() + bolometer_camera = BolometerCamera(camera_geometry=camera_box) + # The bolometer slit in this instance just contains targeting information + # for the ray tracing, since we have already given our camera a geometry + # The slit is defined in the local coordinate system of the camera + slit = BolometerSlit(slit_id="Example slit", centre_point=ORIGIN, + basis_x=XAXIS, dx=SLIT_WIDTH, basis_y=YAXIS, dy=SLIT_HEIGHT, + parent=bolometer_camera) + for j, angle in enumerate(sensor_angles): + # 4 bolometer foils, spaced at equal intervals along the local X axis + sensor = Node(name="Bolometer sensor", parent=bolometer_camera, + transform=rotate_y(angle) * translate(0, 0, -slit_sensor_separation)) + for i, shift in enumerate([-1.5, -0.5, 0.5, 1.5]): + # Note that the foils will be parented to the camera rather than the + # sensor, so we need to define their transform relative to the camera. + foil_transform = sensor.transform * translate(shift * FOIL_SEPARATION, 0, 0) + foil = BolometerFoil(detector_id="Foil {} sensor {}".format(i + 1, j + 1), + centre_point=ORIGIN.transform(foil_transform), + basis_x=XAXIS.transform(foil_transform), dx=FOIL_WIDTH, + basis_y=YAXIS.transform(foil_transform), dy=FOIL_HEIGHT, + slit=slit, parent=bolometer_camera, units="Power", + accumulate=False, curvature_radius=FOIL_CORNER_CURVATURE) + bolometer_camera.add_foil_detector(foil) + return bolometer_camera + + +def load_bolometers(parent=None): + """ + Load the Generomak bolometers. + + The Generomak bolometer diagnostic consists of multiple 12-channel + cameras. Each camera has 3 4-channel sensors inside. + + * 2 cameras are located at the midplane with purely-poloidal, + horizontal views. + * 1 camera is located at the top of the machine with purely-poloidal, + vertical views. + * 2 cameras have purely tangential views at the midplane. + * 1 camera has combined poloidal+tangential views, which look like + curved lines of sight in the poloidal plane. It looks at the lower + divertor. + + :param parent: the scenegraph node the bolometers will belong to. + :return: a list of BolometerCamera instances, one for each of the + cameras described above. + """ + poloidal_camera_rotations = [ + 30, -30, # Horizontal poloidal, + -90, # Vertical poloidal, + 0, # Tangential, + 0, # Combined poloidal/tangential, + ] + toroidal_camera_rotations = [ + 0, 0, # Horizontal poloidal + 0, # Vertical poloidal + -40, # Tangential + 40, # Combined poloidal/tangential + ] + radial_camera_rotations = [ + 0, 0, # Horizontal poloidal + 0, # Vertical poloidal + 90, # Tangential + 90, # Combined poloidal/tangential + ] + camera_origins = [ + Point3D(2.5, 0.05, 0), Point3D(2.5, -0.05, 0), # Horizontal poloidal + Point3D(1.3, 0, 1.42), # Vertical poloidal + Point3D(2.5, 0, 0), # Midplane tangential horizontal + Point3D(2.5, 0, -0.2), # Combined poloidal/tangential + ] + slit_sensor_separations = [ + 0.08, 0.08, # Horizontal poloidal + 0.05, # Vertical poloidal + 0.1, # Tangential + 0.1, # Combined poloidal/tangential + ] + all_sensor_angles = [ + [-22.5, -7.5, 7.5, 22.5], [-22.5, -7.5, 7.5, 22.5], # Horizontal poloidal + [-36, -12, 12, 36], # Vertical poloidal + [-18, -6, 6, 18], # Tangential + [-18, -6, 6, 18], # Combined poloidal/tangential + ] + toroidal_angles = [ + 10, 10, # Horizontal poloidal, need to avoid LFS limiters. + 0, # Vertical poloidal, happy to hit LFS limiters. + -15, # Tangential, avoid LFS limiters. + 15, # Combined poloidal/tangential, avoid LFS limiters. + ] + names = [ + 'HozPol1', 'HozPol2', + 'VertPol', + 'TanMid1', + 'TanPol1', + ] + cameras = [] + # FIXME: this for loop definition is really ugly! + for ( + poloidal_rotation, toroidal_rotation, radial_rotation, camera_origin, + slit_sensor_separation, sensor_angles, toroidal_angle, name + ) in zip( + poloidal_camera_rotations, toroidal_camera_rotations, radial_camera_rotations, + camera_origins, + slit_sensor_separations, all_sensor_angles, toroidal_angles, names + ): + camera = _make_bolometer_camera(slit_sensor_separation, sensor_angles) + # FIXME: work out how to combine tangential and poloidal rotations. + camera.transform = ( + rotate_z(toroidal_angle) + * translate(camera_origin.x, camera_origin.y, camera_origin.z) + * rotate_z(toroidal_rotation) + * rotate_x(radial_rotation) + * rotate_y(poloidal_rotation + 90) + * rotate_basis(-ZAXIS, YAXIS) + ) + camera.parent = parent + camera.name = "{} at angle {}".format(name, poloidal_rotation) + cameras.append(camera) + return cameras diff --git a/cherab/tools/inversions/admt_utils.py b/cherab/tools/inversions/admt_utils.py index 549c4d42..74fac859 100644 --- a/cherab/tools/inversions/admt_utils.py +++ b/cherab/tools/inversions/admt_utils.py @@ -36,15 +36,15 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Generate the first and second derivative operators for a regular grid. :param ndarray voxel_vertices: an Nx4x2 array of coordinates of the - vertices of each voxel, (R, Z) + vertices of each voxel, (R, Z) :param dict grid_1d_to_2d_map: a mapping from the 1D array of - voxels in the grid to a 2D array of voxels if they were arranged - spatially. + voxels in the grid to a 2D array of voxels if they were arranged + spatially. :param dict grid_2d_to_1d_map: the inverse mapping from a 2D - spatially-arranged array of voxels to the 1D array. + spatially-arranged array of voxels to the 1D array. - :return dict operators: a dictionary containing the derivative - operators: Dij for i, y ∊ (x, y) and Di for i ∊ (x, y). + :return: a dictionary containing the derivative operators: Dij for + i, j ∊ (x, y) and Di for i ∊ (x, y), Dsp and Dsm. This function assumes that all voxels are rectilinear, with their axes aligned to the coordinate axes. Additionally, all voxels are @@ -62,13 +62,18 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, D_{xx} \equiv \frac{\partial^2}{\partial x^2}\\ D_{xy} \equiv \frac{\partial^2}{\partial x \partial y} - etc. + etc. It also produces two additional operators, Dsp and Dsm, for + second derivatives on the dy/dx = 1 and dy/dx = -1 diagonals + respectively. Note that the standard 2D laplacian (for isotropic regularisation) - can be trivially calculated as L = Dxx * dx + Dyy * dy, where dx and - dy are the voxel width and height respectively. This expression does - not however produce the 2D laplacian derived from the N-dimensional - case. + can be trivially calculated as follows: + + .. math:: + L = (1 - \alpha) (D_{xx} + D_{yy}) + (\alpha / 2) (D_{sp} + D_{sm}) + + α = 2/3 produces the operator used in Carr et. al. RSI 89, 083506 (2018). + α = 1/3 produces the operator with optimal isotropy. """ # Input argument validation: assume rectilinear voxels voxel_vertices = np.asarray(voxel_vertices) @@ -87,6 +92,8 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx = np.zeros((num_cells, num_cells)) Dxy = np.zeros((num_cells, num_cells)) Dyy = np.zeros((num_cells, num_cells)) + Dsp = np.zeros((num_cells, num_cells)) + Dsm = np.zeros((num_cells, num_cells)) # TODO: for now, we assume all voxels have rectangular cross sections # which are approximately identical. As per Ingesson's notation, we # assume voxels are ordered from top left to bottom right, in column-major @@ -99,6 +106,9 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, dx = np.min(abs(dx[dx != 0])).item() dy = np.min(abs(dy[dy != 0])).item() + # Work out how the voxels are ordered: increasing/decreasing in x/y. + xinc, yinc = np.sign(cell_centres[-1] - cell_centres[0]) + # Note that iy increases as y decreases (cells go from top to bottom), # which is the same as Ingesson's notation in equations 37-41 # Use the second version of the second derivative boundary formulae, so @@ -110,8 +120,13 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, # get the 2D mesh coordinates of this cell ix, iy = grid_index_1d_to_2d_map[ith_cell] + iright = ix + xinc + ileft = ix - xinc + iabove = iy + yinc + ibelow = iy - yinc + try: - n_left = grid_index_2d_to_1d_map[ix - 1, iy] # left of n0 + n_left = grid_index_2d_to_1d_map[ileft, iy] # left of n0 except KeyError: at_left = True else: @@ -119,7 +134,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx[ith_cell, n_left] = 1 try: - n_below_left = grid_index_2d_to_1d_map[ix - 1, iy + 1] # below left of n0 + n_below_left = grid_index_2d_to_1d_map[ileft, ibelow] # below left of n0 except KeyError: # KeyError does not necessarily mean bottom AND left pass @@ -127,7 +142,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxy[ith_cell, n_below_left] = 1 / 4 try: - n_below = grid_index_2d_to_1d_map[ix, iy + 1] + n_below = grid_index_2d_to_1d_map[ix, ibelow] except KeyError: at_bottom = True else: @@ -135,14 +150,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dyy[ith_cell, n_below] = 1 try: - n_below_right = grid_index_2d_to_1d_map[ix + 1, iy + 1] + n_below_right = grid_index_2d_to_1d_map[iright, ibelow] except KeyError: pass else: Dxy[ith_cell, n_below_right] = -1 / 4 try: - n_right = grid_index_2d_to_1d_map[ix + 1, iy] + n_right = grid_index_2d_to_1d_map[iright, iy] except KeyError: at_right = True else: @@ -150,14 +165,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxx[ith_cell, n_right] = 1 try: - n_above_right = grid_index_2d_to_1d_map[ix + 1, iy - 1] + n_above_right = grid_index_2d_to_1d_map[iright, iabove] except KeyError: pass else: Dxy[ith_cell, n_above_right] = 1 / 4 try: - n_above = grid_index_2d_to_1d_map[ix, iy - 1] + n_above = grid_index_2d_to_1d_map[ix, iabove] except KeyError: at_top = True else: @@ -165,7 +180,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dyy[ith_cell, n_above] = 1 try: - n_above_left = grid_index_2d_to_1d_map[ix - 1, iy - 1] + n_above_left = grid_index_2d_to_1d_map[ileft, iabove] except KeyError: pass else: @@ -247,14 +262,42 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxy[ith_cell, ith_cell] = -1 Dxy[ith_cell, n_above_left] = -1 + if np.isnan(n_above_left) and not np.isnan(n_below_right): + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_below_right] = 1 + elif np.isnan(n_below_right) and not np.isnan(n_above_left): + Dsm[ith_cell, ith_cell] = -1 + Dsm[ith_cell, n_above_left] = 1 + elif np.isnan(n_above_left) and np.isnan(n_below_right): + Dsm[ith_cell, ith_cell] = 0 + else: + Dsm[ith_cell, ith_cell] = -2 + Dsm[ith_cell, n_above_left] = 1 + Dsm[ith_cell, n_below_right] = 1 + + if np.isnan(n_above_right) and not np.isnan(n_below_left): + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_below_left] = 1 + elif np.isnan(n_below_left) and not np.isnan(n_above_right): + Dsp[ith_cell, ith_cell] = -1 + Dsp[ith_cell, n_above_right] = 1 + elif np.isnan(n_below_left) and np.isnan(n_above_right): + Dsp[ith_cell, ith_cell] = 0 + else: + Dsp[ith_cell, ith_cell] = -2 + Dsp[ith_cell, n_above_right] = 1 + Dsp[ith_cell, n_below_left] = 1 + Dx = Dx / dx Dy = Dy / dy Dxx = Dxx / dx**2 Dyy = Dyy / dy**2 Dxy = Dxy / (dx * dy) + Dsp = Dsp / (dx**2 + dy**2) + Dsm = Dsm / (dx**2 + dy**2) # Package all operators up into a dictionary - operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy) + operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy, Dsp=Dsp, Dsm=Dsm) return operators @@ -263,22 +306,21 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani Calculate the ADMT regularisation operator. :param ndarray voxel_radii: a 1D array of the radius at the centre - of each voxel in the grid + of each voxel in the grid :param tuple derivative_operators: a named tuple with the derivative - operators for the grid, as returned by :func:generate_derivative_operators + operators for the grid, as returned by :func:generate_derivative_operators :param ndarray psi_at_voxels: the magnetic flux at the centre of - each voxel in the grid + each voxel in the grid :param float dx: the width of each voxel. :param float dy: the height of each voxel :param float anisotropy: the ratio of the smoothing in the parallel - and perpendicular directions. - - :return ndarray admt: the ADMT regularisation operator. + and perpendicular directions. + :return: the ADMT regularisation operator. The degree of anisotropy dictates the relative suppression of gradients in the directions parallel and perpendicular to the - magnetic field. For example, `anisotropy=10` implies parallel - gradients in solution are 10 times smaller than perpendicular + magnetic field. For example, ``anisotropy=10`` implies parallel + gradients in the solution are 10 times smaller than perpendicular gradients. This function assumes that all voxels are rectilinear, with their diff --git a/demos/observers/bolometry/admt_tomographic_inversion.py b/demos/observers/bolometry/admt_tomographic_inversion.py new file mode 100644 index 00000000..e27fc518 --- /dev/null +++ b/demos/observers/bolometry/admt_tomographic_inversion.py @@ -0,0 +1,285 @@ +""" +This example demonstrates performing a tomographic reconstruction of a +radiation profile using Cherab's anisotropic diffusion (ADMT) regularisation +utilities. We use the machine geometry, sample bolometers and equilibrium +from Generomak. +""" +import matplotlib.pyplot as plt +import numpy as np + +from raysect.core.math.function.float import Exp2D, Arg2D, Atan4Q2D +from raysect.core.math import translate +from raysect.optical import World +from raysect.optical.material import AbsorbingSurface, VolumeTransform +from raysect.primitive import Cylinder, Subtract + +from cherab.generomak.machine import load_first_wall +from cherab.generomak.equilibrium import load_equilibrium +from cherab.generomak.diagnostics import load_bolometers +from cherab.core.math import sample2d, sample2d_grid, sample2d_points, AxisymmetricMapper +from cherab.tools.emitters import RadiationFunction +from cherab.tools.raytransfer import RayTransferCylinder, RayTransferPipeline0D +from cherab.tools.inversions import admt_utils as admt +from cherab.tools.inversions import invert_regularised_nnls + + +plt.ion() + +################################################################################ +# Define the emissivity profile. +################################################################################ +# The emissivity profile consists of a blob, a ring and part of a ring on the LFS. +# The blob and the ring are Gaussian flux functions. +# The ring is Gaussian in flux and poloidal angle. +# All have equal maximum emissivities, but not necessarily equal total power. +# We use Raysect's function framework to specify an analytic form for the +# emissivity profile, as this is very quick to sample and ray trace. +eq = load_equilibrium() +psin = eq.psi_normalised +axis = eq.magnetic_axis +blob_centre_psin = 0 +blob_width_psin = 0.1 +blob = Exp2D(-0.5 * (psin - blob_centre_psin)**2 / (blob_width_psin**2)) +ring_centre_psin = 0.5 +ring_width_psin = 0.05 +ring = Exp2D(-0.5 * (psin - ring_centre_psin)**2 / (ring_width_psin**2)) +theta = Atan4Q2D(Arg2D('y') - axis.y, Arg2D('x') - axis.x) +lfs_centre_psin = 0.85 +lfs_width_psin = 0.1 +lfs_centre_theta = 0 +lfs_width_theta = 0.5 +lfs = Exp2D(-0.5 * (((psin - lfs_centre_psin) / lfs_width_psin)**2 + + ((theta - lfs_centre_theta) / lfs_width_theta)**2)) +emissivity = blob + ring + lfs +# Assume no emission from these contributors outside the separatrix. +emissivity = emissivity * eq.inside_lcfs + +# Visualise the emissivity profile with the equilibrium overlayed. +plt.figure() +rsamp, zsamp, psisamp = sample2d(psin, (*eq.r_range, 500), (*eq.z_range, 1000)) +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +rsamp, zsamp, emsamp = sample2d(emissivity, (*eq.r_range, 500), (*eq.z_range, 1000)) +im = plt.imshow(emsamp.T, extent=(rsamp[0], rsamp[-1], zsamp[0], zsamp[-1]), cmap='Purples') +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Model emissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.pause(0.5) + + +################################################################################ +# Load the machine wall and diagnostic. +################################################################################ +print("Loading the geometry...") +world = World() +load_first_wall(world, material=AbsorbingSurface()) +bolos = load_bolometers(world) +# Only consider the purely-poloidal cameras for now... +bolos = bolos[:3] + +######################################################################## +# Produce a voxel grid +######################################################################## +print("Producing the voxel grid...") +# Define the centres of each voxel, as an (nx, ny, 2) array. +nx = 40 +ny = 85 +cell_r, cell_dx = np.linspace(0.7, 2.5, nx, retstep=True) +cell_z, cell_dz = np.linspace(-1.8, 1.6, ny, retstep=True) +cell_r_grid, cell_z_grid = np.broadcast_arrays(cell_r[:, None], cell_z[None, :]) +cell_centres = np.stack((cell_r_grid, cell_z_grid), axis=-1) # (nx, ny, 2) array + +# Define the positions of the vertices of the voxels. +cell_vertices_r = np.linspace(cell_r[0] - 0.5 * cell_dx, cell_r[-1] + 0.5 * cell_dx, nx + 1) +cell_vertices_z = np.linspace(cell_z[0] - 0.5 * cell_dz, cell_z[-1] + 0.5 * cell_dz, ny + 1) + +# Build a mask, only including cells within the wall. +mask_2d = sample2d_grid(eq.inside_limiter, cell_r, cell_z) +mask_3d = mask_2d[:, np.newaxis, :] +ncells = mask_3d.sum() + +# We'll use the Ray Transfer frameworks as these voxels are rectangular +# and it's much faster than the Voxel framework for simple cases like this. +ray_transfer_grid = RayTransferCylinder( + radius_outer=cell_vertices_r[-1], + radius_inner=cell_vertices_r[0], + height=cell_vertices_z[-1] - cell_vertices_z[0], + n_radius=nx, n_height=ny, mask=mask_3d, n_polar=1, + transform=translate(0, 0, cell_vertices_z[0]), +) + +######################################################################## +# Calculate the geometry matrix for the grid +######################################################################## +print("Calculating the geometry matrix...") +# The ray transfer object must be in the same world as the bolometers +ray_transfer_grid.parent = world + +sensitivity_matrix = [] +for camera in bolos: + for foil in camera: + # Temporarily override foil pipelines for the sensitivity calculation. + orig_pipelines = foil.pipelines + foil.pipelines = [RayTransferPipeline0D(kind=foil.units)] + # All objects in world have wavelength-independent material properties, + # so it doesn't matter which wavelength range we use (as long as + # max_wavelength - min_wavelength = 1) + foil.min_wavelength = 1 + foil.max_wavelength = 2 + foil.spectral_bins = ray_transfer_grid.bins + foil.observe() + sensitivity_matrix.append(foil.pipelines[0].matrix) + # Restore original pipelines for subsequent observe calls. + foil.pipelines = orig_pipelines +sensitivity_matrix = np.asarray(sensitivity_matrix) + + +################################################################################ +# Generate the regularisation operators. +################################################################################ +print("Generating regularisation operators...") +# generate_derivative_operators requires two mappings, one from a flat list of +# voxels to the original 2D grid, and one for the 2D grid coordinates to the +# flat list of voxels. We could build these by hand, but the RayTransferCylinder +# object helpfully provides the data already so we just need to convert from +# arrays to dictionaries. +grid_index_1d_to_2d_map = {} +for k, idx2d in enumerate(ray_transfer_grid.invert_voxel_map()): + # We want the x and z elements, as the Ray Transfer grid is 3D and this + # inversion is going to be in 2D. + grid_index_1d_to_2d_map[k] = (idx2d[0].item(), idx2d[2].item()) +grid_index_2d_to_1d_map = {} +nx, _, ny = ray_transfer_grid.voxel_map.shape +for i in range(nx): + for j in range(ny): + voxel_index = ray_transfer_grid.voxel_map[i, 0, j] + if voxel_index != -1: + grid_index_2d_to_1d_map[(i, j)] = voxel_index +# We now need an (Nx4x2) array of voxel vertices, which can be easily calculated. +voxel_centres = np.array([cell_centres[grid_index_1d_to_2d_map[i]] + for i in range(ray_transfer_grid.bins)]) +vertex_displacements = np.array([[-cell_dx/2, -cell_dz/2], + [-cell_dx/2, cell_dz/2], + [cell_dx/2, cell_dz/2], + [cell_dx/2, -cell_dz/2]]) +# Combine the (N,2) and (4,2) arrays to get an (N,4,2) array. +voxel_vertices = voxel_centres[:, None, :] + vertex_displacements[None, :, :] +derivative_operators = admt.generate_derivative_operators( + voxel_vertices, grid_index_1d_to_2d_map, grid_index_2d_to_1d_map +) + +# As described in the docstring for generate_derivative_operators, we can +# calculate a 2D laplacian operator for "isotropic" smoothing easily: +alpha = 1/3 # Optimal isotropy +aligned = derivative_operators['Dxx'] * cell_dx**2 + derivative_operators['Dyy'] * cell_dz**2 +skewed = (derivative_operators['Dsp'] + derivative_operators['Dsm']) * (cell_dx**2 + cell_dz**2) +laplacian = (1 - alpha) * aligned + (alpha / 2) * skewed +# We could also use alpha = 2/3, which would produce an operator akin to the one +# used in Carr et. al. RSI 89, 083506 (2018). + +# We can also derive an anistoropic regularisation operator, which calculates the +# amount of un-smoothness parallel and perpendicular to the magnetic field lines. +# For this we need the radii of the voxels and the magnetic flux at each voxel, +# along with a few other inputs. +voxel_radii = voxel_centres[:, 0] +psi_at_voxels = sample2d_points(eq.psi_normalised, voxel_centres) +# We also need to decide on the degree of anisotropy we expect, i.e. how much more +# smooth the radiation is along the field lines vs perpendicular to them. +# The optimal value will depend on the problem at hand. +anisotropy = 50 +admt_operator = admt.calculate_admt( + voxel_radii, derivative_operators, psi_at_voxels, cell_dx, cell_dz, anisotropy +) + +################################################################################ +# Forward model the measurements. +################################################################################ +print("Modelling the measurement values...") +# Create an emitting object whose emission is defined by the analytic form we +# produced earlier. As the emission depends on the equilibrium, this object +# should have an extent no larger than the equilibrium reconstruction extent. +# We actually make the emitter slightly smaller than the equilibrium region to +# avoid numerical precision issues creating attempts to calculate the emissivity +# outside of the equlibrium domain. +CYLINDER_RADIUS = eq.r_range[-1] - 1e-6 +CYLINDER_HEIGHT = eq.z_range[-1] - eq.z_range[0] - 2e-6 +CYLINDER_SHIFT = eq.z_range[0] + 1e-6 +emitter = Cylinder(radius=CYLINDER_RADIUS, height=CYLINDER_HEIGHT, + transform=translate(0, 0, CYLINDER_SHIFT)) +# Cut out middle of cylinder as well: equilibrium not defined here. +emitter = Subtract(emitter, Cylinder(radius=eq.r_range[0] + 1e-6, height=10, + transform=translate(0, 0, -5))) +emission_function_3d = AxisymmetricMapper(emissivity) +emitting_material = VolumeTransform(RadiationFunction(emission_function_3d), + transform=emitter.transform.inverse()) +emitter.material = emitting_material +emitter.parent = world + +# Calculate the line-integral bolometer measurements by observing the emitter +# with all bolometers. The measurements should have the same channel order as +# the sensitivity matrix. +all_measurements = [] +for camera in bolos: + all_measurements.extend(camera.observe()) + + +################################################################################ +# Perform the inversions. +################################################################################ +print("Performing inversions...") +# We'll use NNLS with regularisation. The hyperparameters have been chosen by +# hand but techniques such as the discrepancy principle or L curve optimisation +# could also be used to determine them. That is out of the scope of this demo. +isotropic_alpha = 1e-10 +isotropic_inversion, _ = invert_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=isotropic_alpha, + tikhonov_matrix=laplacian, +) + +admt_alpha = 1e-10 +admt_inversion, _ = invert_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=admt_alpha, + tikhonov_matrix=admt_operator, +) + + +################################################################################ +# Plot the inversion results. +################################################################################ +emiss2d = np.zeros((nx, ny)) + +# Isotropic +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = isotropic_inversion[index1d] +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Isotropic regularisation") + +# Anisotropic. +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = admt_inversion[index1d] +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Anisotropic regularisation") + +plt.ioff() +plt.show() diff --git a/docs/source/tools/tomography.rst b/docs/source/tools/tomography.rst index 3f3984f6..523f41a1 100644 --- a/docs/source/tools/tomography.rst +++ b/docs/source/tools/tomography.rst @@ -119,3 +119,36 @@ Use spectral pipelines from Raysect if you need these features. .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline1D .. autoclass:: cherab.tools.raytransfer.pipelines.RayTransferPipeline2D + + +Regularisation +-------------- + +Some of the inversion methods take a regularisation operator, which provides +additional constraints to help achieved unique solutions to ill-posed +tomography problems. Many regularisation schemes impose constraints on the smoothness +of the resulting solution, with this smoothness quantified by the second derivative +of the solution. Two such regularisation schemes are common in fusion applications: + +#. Isotropic smoothing, where the solution has the same smoothness in all directions. +#. Anisotropic smoothing, so-called "anisotropic diffusion model tomography" (ADMT), + where the solution is smoother parallel to the magnetic field and less smooth + perpendicular to the magnetic field. + +Cherab provides some utility functions to assist in calculating appropriate +operators using these (and other) derivative-based regularisation schemes. These can be used +on inversion grids defined using both the Voxel and Ray Transfer frameworks, and passed +directly to the inversion methods in Cherab which take regularisation operators, such as +cherab.tools.inversions.invert_constrained_sart and cherab.tools.inversions.invert_regularised_nnls. + +The routines to calculate derivative operators for inversion grids, and further to calculate +the ADMT operator for a given set of derivative operators and magnetic field, are taken from +work published by L. C. Ingesson in `JET-R(99)08`_. + + +.. autofunction:: cherab.tools.inversions.admt_utils.generate_derivative_operators + +.. autofunction:: cherab.tools.inversions.admt_utils.calculate_admt + + +.. _JET-R(99)08: http://www.euro-fusionscipub.org/wp-content/uploads/2014/11/JETR99008.pdf From e7cfa3dad7176222495fab3de5d0a55b3dde4da6 Mon Sep 17 00:00:00 2001 From: "Lovell, Jack J" Date: Tue, 19 May 2026 15:56:35 +0200 Subject: [PATCH 2/9] Tidy up the for loop creating generomak bolometer sensors --- cherab/generomak/diagnostics/bolometers.py | 135 ++++++++++----------- 1 file changed, 65 insertions(+), 70 deletions(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e4693712..e64d9dec 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -20,6 +20,7 @@ BOX_WIDTH = 0.2 BOX_HEIGHT = 0.07 BOX_DEPTH = 0.2 +THICKNESS = 1e-3 SLIT_WIDTH = 0.004 SLIT_HEIGHT = 0.005 FOIL_WIDTH = 0.0013 @@ -39,13 +40,13 @@ def _make_bolometer_camera(slit_sensor_separation, sensor_angles): """ camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) - # Hollow out the box - inside_box = Box(lower=camera_box.lower + Vector3D(1e-5, 1e-5, 1e-5), - upper=camera_box.upper - Vector3D(1e-5, 1e-5, 1e-5)) + # Hollow out the box: it has 1 mm thick walls. + inside_box = Box(lower=camera_box.lower + Vector3D(THICKNESS, THICKNESS, THICKNESS), + upper=camera_box.upper - Vector3D(THICKNESS, THICKNESS, THICKNESS)) camera_box = Subtract(camera_box, inside_box) - # The slit is a hole in the box - aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1e-4), - upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 1e-4)) + # The slit is a hole in the box. Make it thicker than the wall. + aperture = Box(lower=Point3D(-SLIT_WIDTH / 2, -SLIT_HEIGHT / 2, -1.1 * THICKNESS), + upper=Point3D(SLIT_WIDTH / 2, SLIT_HEIGHT / 2, 0.1 * THICKNESS)) camera_box = Subtract(camera_box, aperture) camera_box.material = AbsorbingSurface() bolometer_camera = BolometerCamera(camera_geometry=camera_box) @@ -93,75 +94,69 @@ def load_bolometers(parent=None): :return: a list of BolometerCamera instances, one for each of the cameras described above. """ - poloidal_camera_rotations = [ - 30, -30, # Horizontal poloidal, - -90, # Vertical poloidal, - 0, # Tangential, - 0, # Combined poloidal/tangential, - ] - toroidal_camera_rotations = [ - 0, 0, # Horizontal poloidal - 0, # Vertical poloidal - -40, # Tangential - 40, # Combined poloidal/tangential - ] - radial_camera_rotations = [ - 0, 0, # Horizontal poloidal - 0, # Vertical poloidal - 90, # Tangential - 90, # Combined poloidal/tangential - ] - camera_origins = [ - Point3D(2.5, 0.05, 0), Point3D(2.5, -0.05, 0), # Horizontal poloidal - Point3D(1.3, 0, 1.42), # Vertical poloidal - Point3D(2.5, 0, 0), # Midplane tangential horizontal - Point3D(2.5, 0, -0.2), # Combined poloidal/tangential - ] - slit_sensor_separations = [ - 0.08, 0.08, # Horizontal poloidal - 0.05, # Vertical poloidal - 0.1, # Tangential - 0.1, # Combined poloidal/tangential - ] - all_sensor_angles = [ - [-22.5, -7.5, 7.5, 22.5], [-22.5, -7.5, 7.5, 22.5], # Horizontal poloidal - [-36, -12, 12, 36], # Vertical poloidal - [-18, -6, 6, 18], # Tangential - [-18, -6, 6, 18], # Combined poloidal/tangential - ] - toroidal_angles = [ - 10, 10, # Horizontal poloidal, need to avoid LFS limiters. - 0, # Vertical poloidal, happy to hit LFS limiters. - -15, # Tangential, avoid LFS limiters. - 15, # Combined poloidal/tangential, avoid LFS limiters. - ] - names = [ - 'HozPol1', 'HozPol2', - 'VertPol', - 'TanMid1', - 'TanPol1', - ] + camera_properties = { + 'HozPol1': {}, # Horizontal poloidal + 'HozPol2': {}, # Horizontal poloidal, + 'VertPol': {}, # Vertical poloidal + 'TanMid1': {}, # Tangential + 'TanPol1': {} # Combined poloidal/tangential + } + # poloidal rotations + camera_properties['HozPol1']['rotation_poloidal'] = 30 + camera_properties['HozPol2']['rotation_poloidal'] = -30 + camera_properties['VertPol']['rotation_poloidal'] = -90 + camera_properties['TanMid1']['rotation_poloidal'] = 0 + camera_properties['TanPol1']['rotation_poloidal'] = 0 + # toroidal rotation + camera_properties['HozPol1']['rotation_toroidal'] = 0 + camera_properties['HozPol2']['rotation_toroidal'] = 0 + camera_properties['VertPol']['rotation_toroidal'] = 0 + camera_properties['TanMid1']['rotation_toroidal'] = -40 + camera_properties['TanPol1']['rotation_toroidal'] = 40 + # radial rotation + camera_properties['HozPol1']['rotation_radial'] = 0 + camera_properties['HozPol2']['rotation_radial'] = 0 + camera_properties['VertPol']['rotation_radial'] = 0 + camera_properties['TanMid1']['rotation_radial'] = 90 + camera_properties['TanPol1']['rotation_radial'] = 90 + # origins + camera_properties['HozPol1']['origin'] = Point3D(2.5, 0.05, 0) + camera_properties['HozPol2']['origin'] = Point3D(2.5, -0.05, 0) + camera_properties['VertPol']['origin'] = Point3D(1.3, 0, 1.42) + camera_properties['TanMid1']['origin'] = Point3D(2.5, 0, 0) + camera_properties['TanPol1']['origin'] = Point3D(2.5, 0, -0.2) + # slit sensor separations + camera_properties['HozPol1']['slit_sensor_separation'] = 0.08 + camera_properties['HozPol2']['slit_sensor_separation'] = 0.08 + camera_properties['VertPol']['slit_sensor_separation'] = 0.05 + camera_properties['TanMid1']['slit_sensor_separation'] = 0.1 + camera_properties['TanPol1']['slit_sensor_separation'] = 0.1 + # sensor angles + camera_properties['HozPol1']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] + camera_properties['HozPol2']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] + camera_properties['VertPol']['sensor_angles'] = [-36, -12, 12, 36] + camera_properties['TanMid1']['sensor_angles'] = [-18, -6, 6, 18] + camera_properties['TanPol1']['sensor_angles'] = [-18, -6, 6, 18] + # toroidal angles + camera_properties['HozPol1']['toroidal_angle'] = 10 # need to avoid LFS limiters + camera_properties['HozPol2']['toroidal_angle'] = 10 # need to avoid LFS limiters + camera_properties['VertPol']['toroidal_angle'] = 0 # happy to hit LFS limiters + camera_properties['TanMid1']['toroidal_angle'] = -15 # avoid LFS limiters + camera_properties['TanPol1']['toroidal_angle'] = 15 # avoid LFS limiters + cameras = [] - # FIXME: this for loop definition is really ugly! - for ( - poloidal_rotation, toroidal_rotation, radial_rotation, camera_origin, - slit_sensor_separation, sensor_angles, toroidal_angle, name - ) in zip( - poloidal_camera_rotations, toroidal_camera_rotations, radial_camera_rotations, - camera_origins, - slit_sensor_separations, all_sensor_angles, toroidal_angles, names - ): - camera = _make_bolometer_camera(slit_sensor_separation, sensor_angles) + for name, prop in camera_properties.items(): + camera = _make_bolometer_camera(prop['slit_sensor_separation'], prop['sensor_angles']) # FIXME: work out how to combine tangential and poloidal rotations. camera.transform = ( - rotate_z(toroidal_angle) - * translate(camera_origin.x, camera_origin.y, camera_origin.z) - * rotate_z(toroidal_rotation) - * rotate_x(radial_rotation) - * rotate_y(poloidal_rotation + 90) + rotate_z(prop['toroidal_angle']) + * translate(prop['origin'].x, prop['origin'].y, prop['origin'].z) + * rotate_z(prop['rotation_toroidal']) + * rotate_x(prop['rotation_radial']) + * rotate_y(prop['rotation_poloidal'] + 90) * rotate_basis(-ZAXIS, YAXIS) ) camera.parent = parent - camera.name = "{} at angle {}".format(name, poloidal_rotation) + camera.name = "{} at angle {}".format(name, prop['rotation_poloidal']) cameras.append(camera) return cameras From e6d233d09be497c8372310de6f2221844926a4cc Mon Sep 17 00:00:00 2001 From: "Lovell, Jack J" Date: Wed, 27 May 2026 12:04:58 +0100 Subject: [PATCH 3/9] Remove unused import --- cherab/generomak/diagnostics/bolometers.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e64d9dec..a2412348 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -3,7 +3,6 @@ """ from raysect.core import (Node, Point3D, Vector3D, rotate_basis, rotate_x, rotate_y, rotate_z, translate) -from raysect.core.math import rotate from raysect.optical.material import AbsorbingSurface from raysect.primitive import Box, Subtract From c7b504f60bf64bb0de768ee11100ce195b52390d Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Wed, 17 Jun 2026 10:26:43 +0100 Subject: [PATCH 4/9] Improve Generomak bolometer geometry * Fix the series of transforms to position and orient the bolometer cameras. * Update geometry parameters such that the lines of sight all follow a consistant convention and specify the channel ordering in the docstring. This required adding the ability to rotate the sensors about the viewing axis to reverse the channel order of the horizontal-poloidal camera. * Correct the docstring of load_bolometers to specify the correct number of 4-channel sensors in each camera. * Expand docstrings and comments to better describe the intricacies of setting up the bolometer geometry. --- cherab/generomak/diagnostics/bolometers.py | 123 +++++++++++++++------ 1 file changed, 91 insertions(+), 32 deletions(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index a2412348..e522fb36 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -16,7 +16,7 @@ ORIGIN = Point3D(0, 0, 0) # Bolometer geometry, independent of camera. BOX_WIDTH = 0.05 -BOX_WIDTH = 0.2 +BOX_WIDTH = 0.1 BOX_HEIGHT = 0.07 BOX_DEPTH = 0.2 THICKNESS = 1e-3 @@ -28,14 +28,23 @@ FOIL_SEPARATION = 0.00508 # 0.2 inch between foils -def _make_bolometer_camera(slit_sensor_separation, sensor_angles): +def _make_bolometer_camera(slit_sensor_separation, sensor_angles, sensor_rotations): """ Build a single bolometer camera. The camera consists of a box with a rectangular slit and 4 sensors, each of which has 4 foils. + In its local coordinate system, the camera's slit is located at the - origin and the sensors below the z=0 plane, looking up towards the slit. + origin with its width along the X axis and its height along the y + axis, and the sensors are below the z=0 plane looking up towards the + slit. + + The sensors are rotated by sensor_angles about the y axis to form a + fan, and by sensor_rotations about the axis defined by the line + between the slit and the sensor. A rotation of 180 degrees flips + the sensor upside down and therefore reverses the spatial ordering + of lines of sight relative to a rotation of 0 degrees. """ camera_box = Box(lower=Point3D(-BOX_WIDTH / 2, -BOX_HEIGHT / 2, -BOX_DEPTH), upper=Point3D(BOX_WIDTH / 2, BOX_HEIGHT / 2, 0)) @@ -55,10 +64,14 @@ def _make_bolometer_camera(slit_sensor_separation, sensor_angles): slit = BolometerSlit(slit_id="Example slit", centre_point=ORIGIN, basis_x=XAXIS, dx=SLIT_WIDTH, basis_y=YAXIS, dy=SLIT_HEIGHT, parent=bolometer_camera) - for j, angle in enumerate(sensor_angles): + for j, (angle, rotation) in enumerate(zip(sensor_angles, sensor_rotations)): # 4 bolometer foils, spaced at equal intervals along the local X axis - sensor = Node(name="Bolometer sensor", parent=bolometer_camera, - transform=rotate_y(angle) * translate(0, 0, -slit_sensor_separation)) + sensor = Node(name="Bolometer sensor", parent=bolometer_camera) + sensor.transform = ( + rotate_y(angle) + * rotate_z(rotation) + * translate(0, 0, -slit_sensor_separation) + ) for i, shift in enumerate([-1.5, -0.5, 0.5, 1.5]): # Note that the foils will be parented to the camera rather than the # sensor, so we need to define their transform relative to the camera. @@ -77,8 +90,8 @@ def load_bolometers(parent=None): """ Load the Generomak bolometers. - The Generomak bolometer diagnostic consists of multiple 12-channel - cameras. Each camera has 3 4-channel sensors inside. + The Generomak bolometer diagnostic consists of multiple 16-channel + cameras. Each camera has 4 4-channel sensors inside. * 2 cameras are located at the midplane with purely-poloidal, horizontal views. @@ -89,10 +102,40 @@ def load_bolometers(parent=None): curved lines of sight in the poloidal plane. It looks at the lower divertor. + Channel ordering is as follows: + * Poloidal channels are ordered anti-clockwise by line-of-sight: + channel 1 of HozPol1 views the top of the machine and channel 16 + HozPol2 views the bottom of the machine. Similarly, channel 1 of + VertPol views the high field side and channel 16 views the low + field side. + * Tangential channels are ordered by increasing tangency radius: + channel 1 of TanMid1 has its tangency radius on the high field + side and channel 16 has its tangency radius on the low field side. + * The combined tangential/poloidal channels follow both conventions: + channel 1 views the high field side and channel 16 views the low + field side. + :param parent: the scenegraph node the bolometers will belong to. :return: a list of BolometerCamera instances, one for each of the cameras described above. """ + # The coordinate system conventions are as follows. All angles are in + # degrees and increase clockwise when viewing along the relevant axes: + # y axis for poloidal rotation, z axis for toroidal rotation and x axis + # for radial rotation. + # - rotation_poloidal: viewing angle of the slit in the poloidal plane, + # with 0 being horizontally inwards. + # - rotation_toroidal: viewing angle of the slit in the toroidal plane, + # with 0 being purely radial. + # - rotation_radial: rotation about the radial axis, 0 being vertically upwards. + # - origin: position of the slit relative to the (x, z) poloidal plane i.e. y=0. + # - slit_sensor_separation: distance between slit and each 4-channel sensor. + # - sensor_angles: angle between slit normal and sensor normal. + # - sensor_rotations: rotation angle about the slit-sensor vector, enables + # reversing the order of lines of sight spatially within + # each sensor. + # - toroidal_angle: the angle of the poloidal plane in which the origin is + # definied, with 0 being the (x, z) plane. camera_properties = { 'HozPol1': {}, # Horizontal poloidal 'HozPol2': {}, # Horizontal poloidal, @@ -105,7 +148,7 @@ def load_bolometers(parent=None): camera_properties['HozPol2']['rotation_poloidal'] = -30 camera_properties['VertPol']['rotation_poloidal'] = -90 camera_properties['TanMid1']['rotation_poloidal'] = 0 - camera_properties['TanPol1']['rotation_poloidal'] = 0 + camera_properties['TanPol1']['rotation_poloidal'] = -25 # toroidal rotation camera_properties['HozPol1']['rotation_toroidal'] = 0 camera_properties['HozPol2']['rotation_toroidal'] = 0 @@ -113,30 +156,36 @@ def load_bolometers(parent=None): camera_properties['TanMid1']['rotation_toroidal'] = -40 camera_properties['TanPol1']['rotation_toroidal'] = 40 # radial rotation - camera_properties['HozPol1']['rotation_radial'] = 0 - camera_properties['HozPol2']['rotation_radial'] = 0 - camera_properties['VertPol']['rotation_radial'] = 0 - camera_properties['TanMid1']['rotation_radial'] = 90 - camera_properties['TanPol1']['rotation_radial'] = 90 - # origins - camera_properties['HozPol1']['origin'] = Point3D(2.5, 0.05, 0) - camera_properties['HozPol2']['origin'] = Point3D(2.5, -0.05, 0) + camera_properties['HozPol1']['rotation_radial'] = -90 + camera_properties['HozPol2']['rotation_radial'] = -90 + camera_properties['VertPol']['rotation_radial'] = -90 + camera_properties['TanMid1']['rotation_radial'] = 0 + camera_properties['TanPol1']['rotation_radial'] = 0 + # origins relative to the poloidal (x, z) plane + camera_properties['HozPol1']['origin'] = Point3D(2.45, 0.05, 0) + camera_properties['HozPol2']['origin'] = Point3D(2.45, -0.05, 0) camera_properties['VertPol']['origin'] = Point3D(1.3, 0, 1.42) camera_properties['TanMid1']['origin'] = Point3D(2.5, 0, 0) - camera_properties['TanPol1']['origin'] = Point3D(2.5, 0, -0.2) - # slit sensor separations + camera_properties['TanPol1']['origin'] = Point3D(2.2, 0, -0.8) + # slit-sensor separations camera_properties['HozPol1']['slit_sensor_separation'] = 0.08 camera_properties['HozPol2']['slit_sensor_separation'] = 0.08 camera_properties['VertPol']['slit_sensor_separation'] = 0.05 camera_properties['TanMid1']['slit_sensor_separation'] = 0.1 - camera_properties['TanPol1']['slit_sensor_separation'] = 0.1 - # sensor angles - camera_properties['HozPol1']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] - camera_properties['HozPol2']['sensor_angles'] = [-22.5, -7.5, 7.5, 22.5] - camera_properties['VertPol']['sensor_angles'] = [-36, -12, 12, 36] - camera_properties['TanMid1']['sensor_angles'] = [-18, -6, 6, 18] - camera_properties['TanPol1']['sensor_angles'] = [-18, -6, 6, 18] - # toroidal angles + camera_properties['TanPol1']['slit_sensor_separation'] = 0.15 + # sensor angles relative to the slit + camera_properties['HozPol1']['sensor_angles'] = [22.5, 7.5, -7.5, -22.5] + camera_properties['HozPol2']['sensor_angles'] = [22.5, 7.5, -7.5, -22.5] + camera_properties['VertPol']['sensor_angles'] = [36, 12, -12, -36] + camera_properties['TanMid1']['sensor_angles'] = [18, 6, -6, -18] + camera_properties['TanPol1']['sensor_angles'] = [-12, -4, 4, 12] + # sensor rotation relative to the slit + camera_properties['HozPol1']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['HozPol2']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['VertPol']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['TanMid1']['sensor_rotations'] = [0, 0, 0, 0] + camera_properties['TanPol1']['sensor_rotations'] = [180, 180, 180, 180] + # toroidal angles about which to rotate the poloidal plane camera_properties['HozPol1']['toroidal_angle'] = 10 # need to avoid LFS limiters camera_properties['HozPol2']['toroidal_angle'] = 10 # need to avoid LFS limiters camera_properties['VertPol']['toroidal_angle'] = 0 # happy to hit LFS limiters @@ -145,17 +194,27 @@ def load_bolometers(parent=None): cameras = [] for name, prop in camera_properties.items(): - camera = _make_bolometer_camera(prop['slit_sensor_separation'], prop['sensor_angles']) - # FIXME: work out how to combine tangential and poloidal rotations. + camera = _make_bolometer_camera( + prop['slit_sensor_separation'], + prop['sensor_angles'], + prop['sensor_rotations'], + ) + # The transform is applied as follows: + # 1. Point the camera along the inward radial direction in the (x, z) plane. + # 2. Make the poloidal, radial and toroidal rotations while the camera is at + # the origin. + # 3. Move the camera to its position relative to the (x, z) plane. + # 4. Rotate the (x, z) plane to the correct toroidal angle. + # Transforms are applied right-to-left (or bottom-to-top with one per line): camera.transform = ( rotate_z(prop['toroidal_angle']) * translate(prop['origin'].x, prop['origin'].y, prop['origin'].z) * rotate_z(prop['rotation_toroidal']) + * rotate_y(prop['rotation_poloidal']) * rotate_x(prop['rotation_radial']) - * rotate_y(prop['rotation_poloidal'] + 90) - * rotate_basis(-ZAXIS, YAXIS) + * rotate_basis(-XAXIS, ZAXIS) ) camera.parent = parent - camera.name = "{} at angle {}".format(name, prop['rotation_poloidal']) + camera.name = name cameras.append(camera) return cameras From 2d838a6f13a1b31b88c7b8d42c86683fa8ba8ff4 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Wed, 17 Jun 2026 10:36:00 +0100 Subject: [PATCH 5/9] Update transform comment with more logical order --- cherab/generomak/diagnostics/bolometers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cherab/generomak/diagnostics/bolometers.py b/cherab/generomak/diagnostics/bolometers.py index e522fb36..05e2115e 100644 --- a/cherab/generomak/diagnostics/bolometers.py +++ b/cherab/generomak/diagnostics/bolometers.py @@ -201,7 +201,7 @@ def load_bolometers(parent=None): ) # The transform is applied as follows: # 1. Point the camera along the inward radial direction in the (x, z) plane. - # 2. Make the poloidal, radial and toroidal rotations while the camera is at + # 2. Make the radial, poloidal and toroidal rotations while the camera is at # the origin. # 3. Move the camera to its position relative to the (x, z) plane. # 4. Rotate the (x, z) plane to the correct toroidal angle. From 6113e861a846633e3759b3f93c77c91ef22e0b81 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Wed, 17 Jun 2026 17:14:49 +0100 Subject: [PATCH 6/9] Improve generate_derivative_operators function * More thorough commenting of the function to make it easier to follow what is going on and why. * Exploit the sparsity of the derivative opreators but using a dictionary-of-keys representation instead of dense numpy arrays. A dictionary is faster to insert single elements into inside the loop than a numpy array too. * Remove the `np.isnan` checkc which are slow for single elements, replace with a cheap check for `None`. * Support returning the operators as Scipy sparse arrays for use downstream. Default to returning as Numpy arrays for backwards compatibility. These changes result in an order of magnitude speedup in this function. --- cherab/tools/inversions/admt_utils.py | 79 ++++++++++++++++++++------- 1 file changed, 60 insertions(+), 19 deletions(-) diff --git a/cherab/tools/inversions/admt_utils.py b/cherab/tools/inversions/admt_utils.py index 74fac859..9ad9be49 100644 --- a/cherab/tools/inversions/admt_utils.py +++ b/cherab/tools/inversions/admt_utils.py @@ -28,10 +28,14 @@ from collections.abc import Mapping import numpy as np +try: + from scipy.sparse import coo_array as coo +except ImportError: # Scipy < 1.8, deprecated from 1.18 + from scipy.sparse import coo_matrix as coo def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, - grid_index_2d_to_1d_map): + grid_index_2d_to_1d_map, sparse=False): r""" Generate the first and second derivative operators for a regular grid. @@ -42,6 +46,8 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, spatially. :param dict grid_2d_to_1d_map: the inverse mapping from a 2D spatially-arranged array of voxels to the 1D array. + :param sparse: return the operators as sparse matrices if True, or + as dense matrices if False. :return: a dictionary containing the derivative operators: Dij for i, j ∊ (x, y) and Di for i ∊ (x, y), Dsp and Dsm. @@ -87,13 +93,14 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, num_cells = voxel_vertices.shape[0] cell_centres = np.mean(voxel_vertices, axis=1) # Individual derivative operators - Dx = np.zeros((num_cells, num_cells)) - Dy = np.zeros((num_cells, num_cells)) - Dxx = np.zeros((num_cells, num_cells)) - Dxy = np.zeros((num_cells, num_cells)) - Dyy = np.zeros((num_cells, num_cells)) - Dsp = np.zeros((num_cells, num_cells)) - Dsm = np.zeros((num_cells, num_cells)) + # Store derivative operators in dictionary-of-keys sparse array format. + Dx = {} + Dy = {} + Dxx = {} + Dxy = {} + Dyy = {} + Dsp = {} + Dsm = {} # TODO: for now, we assume all voxels have rectangular cross sections # which are approximately identical. As per Ingesson's notation, we # assume voxels are ordered from top left to bottom right, in column-major @@ -115,8 +122,9 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, # that we only need to consider nearest neighbours for ith_cell in range(num_cells): at_top, at_bottom, at_left, at_right = False, False, False, False - n_left, n_right, n_below, n_above = np.nan, np.nan, np.nan, np.nan - n_above_left, n_above_right, n_below_left, n_below_right = np.nan, np.nan, np.nan, np.nan + n_left, n_right, n_below, n_above = None, None, None, None + n_above_left, n_above_right, n_below_left, n_below_right = None, None, None, None + # get the 2D mesh coordinates of this cell ix, iy = grid_index_1d_to_2d_map[ith_cell] @@ -125,6 +133,7 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, iabove = iy + yinc ibelow = iy - yinc + # Handle voxels not at the edges/corners of the grid. try: n_left = grid_index_2d_to_1d_map[ileft, iy] # left of n0 except KeyError: @@ -186,14 +195,18 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, else: Dxy[ith_cell, n_above_left] = -1 / 4 + + # Cases which are the same throughout the matrix. + Dxx[ith_cell, ith_cell] = -2 + Dyy[ith_cell, ith_cell] = -2 + + + # Handle cases at the edges/corners top_left = at_top and at_left top_right = at_top and at_right bottom_left = at_bottom and at_left bottom_right = at_bottom and at_right - Dxx[ith_cell, ith_cell] = -2 - Dyy[ith_cell, ith_cell] = -2 - if at_left: Dx[ith_cell, ith_cell] = -1 Dx[ith_cell, n_right] = 1 @@ -262,32 +275,51 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dxy[ith_cell, ith_cell] = -1 Dxy[ith_cell, n_above_left] = -1 - if np.isnan(n_above_left) and not np.isnan(n_below_right): + + # Handle the "skewed" operators. + if n_above_left is None and n_below_right is not None: Dsm[ith_cell, ith_cell] = -1 Dsm[ith_cell, n_below_right] = 1 - elif np.isnan(n_below_right) and not np.isnan(n_above_left): + elif n_below_right is None and n_above_left is not None: Dsm[ith_cell, ith_cell] = -1 Dsm[ith_cell, n_above_left] = 1 - elif np.isnan(n_above_left) and np.isnan(n_below_right): + elif n_above_left is None and n_below_right is None: Dsm[ith_cell, ith_cell] = 0 else: Dsm[ith_cell, ith_cell] = -2 Dsm[ith_cell, n_above_left] = 1 Dsm[ith_cell, n_below_right] = 1 - if np.isnan(n_above_right) and not np.isnan(n_below_left): + if n_above_right is None and n_below_left is not None: Dsp[ith_cell, ith_cell] = -1 Dsp[ith_cell, n_below_left] = 1 - elif np.isnan(n_below_left) and not np.isnan(n_above_right): + elif n_below_left is None and n_above_right is not None: Dsp[ith_cell, ith_cell] = -1 Dsp[ith_cell, n_above_right] = 1 - elif np.isnan(n_below_left) and np.isnan(n_above_right): + elif n_below_left is None and n_above_right is None: Dsp[ith_cell, ith_cell] = 0 else: Dsp[ith_cell, ith_cell] = -2 Dsp[ith_cell, n_above_right] = 1 Dsp[ith_cell, n_below_left] = 1 + + # Although we've stored the operators as dictionaries of keys, it turns out to be + # more convenient to construct a COOrdinate sparse matrix rather than a DOK one + # in Scipy. We then convert that to CSR representation for efficient numerical + # operations later. + def dok_to_sparse(D): + row, col = zip(*D.keys()) + vals = list(D.values()) + return coo((vals, (row, col)), shape=(num_cells, num_cells)).tocsr() + + Dx = dok_to_sparse(Dx) + Dy = dok_to_sparse(Dy) + Dxx = dok_to_sparse(Dxx) + Dyy = dok_to_sparse(Dyy) + Dxy = dok_to_sparse(Dxy) + Dsp = dok_to_sparse(Dsp) + Dsm = dok_to_sparse(Dsm) Dx = Dx / dx Dy = Dy / dy Dxx = Dxx / dx**2 @@ -296,6 +328,15 @@ def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, Dsp = Dsp / (dx**2 + dy**2) Dsm = Dsm / (dx**2 + dy**2) + # If the user requests dense matrices, convert them after performing all the scaling. + if not sparse: + Dx = Dx.toarray() + Dy = Dy.toarray() + Dxx = Dxx.toarray() + Dxy = Dxy.toarray() + Dsp = Dsp.toarray() + Dsm = Dsm.toarray() + # Package all operators up into a dictionary operators = dict(Dx=Dx, Dy=Dy, Dxx=Dxx, Dyy=Dyy, Dxy=Dxy, Dsp=Dsp, Dsm=Dsm) return operators From 388c1fcedd231f1e9f40bf7ba44345aa0d993974 Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Tue, 7 Jul 2026 16:09:10 +0100 Subject: [PATCH 7/9] Support generating sparse ADMT operators This is done automatically if the derivative operators are sparse. A dense operator is returned otherwise. --- cherab/tools/inversions/admt_utils.py | 27 +++++++++++++++++++-------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/cherab/tools/inversions/admt_utils.py b/cherab/tools/inversions/admt_utils.py index 9ad9be49..abd1a9ad 100644 --- a/cherab/tools/inversions/admt_utils.py +++ b/cherab/tools/inversions/admt_utils.py @@ -28,10 +28,11 @@ from collections.abc import Mapping import numpy as np +from scipy.sparse import issparse try: - from scipy.sparse import coo_array as coo + from scipy.sparse import coo_array as coo, diags_array as diags except ImportError: # Scipy < 1.8, deprecated from 1.18 - from scipy.sparse import coo_matrix as coo + from scipy.sparse import coo_matrix as coo, diags def generate_derivative_operators(voxel_vertices, grid_index_1d_to_2d_map, @@ -348,7 +349,7 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani :param ndarray voxel_radii: a 1D array of the radius at the centre of each voxel in the grid - :param tuple derivative_operators: a named tuple with the derivative + :param dict derivative_operators: a dictionary with the derivative operators for the grid, as returned by :func:generate_derivative_operators :param ndarray psi_at_voxels: the magnetic flux at the centre of each voxel in the grid @@ -377,6 +378,10 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani This means it is suitable for use in Cherab's inversion methods, such as NNLS and SART. + + If the derivative operators are sparse matrices, the returned admt + operator is also a sparse matrix. Otherwise a dense matrix is + returned. """ Dpar = np.full(psi_at_voxels.shape, 1) Dperp = Dpar / anisotropy @@ -428,11 +433,17 @@ def calculate_admt(voxel_radii, derivative_operators, psi_at_voxels, dx, dy, ani + (Dperp - Dpar) * (dpsidxdy * dpsidx + dpsidxx * dpsidy) + ddiff_term_cy + dnorm_term_cy + toroidal_term_cy ) / normalisation - cx = np.diag(cx) - cy = np.diag(cy) - cxx = np.diag(cxx) - cyy = np.diag(cyy) - cxy = np.diag(cxy) + if all(issparse(d) for d in derivative_operators.values()): + # Make sparse versions of the diagonal matrices. + diag = diags + else: + # Dense versions using Numpy. + diag = np.diag + cx = diag(cx) + cy = diag(cy) + cxx = diag(cxx) + cyy = diag(cyy) + cxy = diag(cxy) admt_operator = cx @ Dx + cy @ Dy + cxx @ Dxx + 2 * cxy @ Dxy + cyy @ Dyy admt_operator *= np.sqrt(dx * dy) return admt_operator From 50048d8f69db8b816aeb5063f48464723495176a Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Tue, 7 Jul 2026 16:10:27 +0100 Subject: [PATCH 8/9] Add a regularised NNLS inversion using sparse matrices An alternative to `invert_regularised_nnls` using sparse weight and penalty matrices. The numerical algorithm used is slightly different to `scipy.optimize.nnls`, as it uses the TRF variant of a bounded `scipy.optimize.lsq_linear` with the lower bound set to 0 to enforce positivity. Since the results differ due to floating point precision, the sparse variant is implemented as a separate function to maintain backwards compatibility. --- cherab/tools/inversions/__init__.py | 2 +- cherab/tools/inversions/nnls.py | 61 +++++++++++++++++++++++++++++ 2 files changed, 62 insertions(+), 1 deletion(-) diff --git a/cherab/tools/inversions/__init__.py b/cherab/tools/inversions/__init__.py index 00b67a9b..128a0a0d 100644 --- a/cherab/tools/inversions/__init__.py +++ b/cherab/tools/inversions/__init__.py @@ -19,7 +19,7 @@ from .sart import invert_sart, invert_constrained_sart from .opencl import SartOpencl -from .nnls import invert_regularised_nnls +from .nnls import invert_regularised_nnls, invert_sparse_regularised_nnls from .lstsq import invert_regularised_lstsq from .svd import invert_svd from .voxels import Voxel, AxisymmetricVoxel, VoxelCollection, ToroidalVoxelGrid, UnityVoxelEmitter diff --git a/cherab/tools/inversions/nnls.py b/cherab/tools/inversions/nnls.py index 34779f71..4daf9cf9 100644 --- a/cherab/tools/inversions/nnls.py +++ b/cherab/tools/inversions/nnls.py @@ -19,6 +19,10 @@ import numpy as np import scipy +try: + from scipy.sparse import lil_array as lil, eye_array as eye +except ImportError: # Scipy < 1.8, deprecated from 1.18 + from scipy.sparse import lil_matrix as lil, eye def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs): @@ -70,3 +74,60 @@ def invert_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None x_vector, rnorm = scipy.optimize.nnls(c_matrix / vmax, d_vector / vmax, **kwargs) return x_vector, rnorm * vmax + + +def invert_sparse_regularised_nnls(w_matrix, b_vector, alpha=0.01, tikhonov_matrix=None, **kwargs): + r""" + Solves :math:`\mathbf{b} = \mathbf{W} \mathbf{x}` for the vector :math:`\mathbf{x}`, + using Tikhonov regulariastion. + + This is a thin wrapper around scipy.optimize.lsq_linear which modifies + the arguments to include the supplied Tikhonov regularisation matrix and + enforces bounds to avoid negativity. + + The values of w_matrix, b_vector and alpha * tikhonov_matrix are notmalised + by max(b_vector) before passing them to scipy.optimize.nnls(). + + :param w_matrix: The sensitivity matrix describing the coupling between the + detectors and the voxels. Must be an array with shape :math:`(N_d, N_s)`. May be either + a dense array or a sparse matrix or array. + :param np.ndarray b_vector: The measured power/radiance vector with shape :math:`(N_d)`. + :param float alpha: The regularisation hyperparameter :math:`\alpha` which determines + the regularisation strength of the tikhonov matrix. + :param np.ndarray tikhonov_matrix: The tikhonov regularisation matrix operator, an array + with shape :math:`(N_s, N_s)`. If None, the identity matrix is used. + :param \**kwargs: Keyword arguments passed to scipy.optimize.nnls. + :return: (x, norm), the solution vector and the residual norm. + + .. code-block:: pycon + + >>> from cherab.tools.inversions import invert_regularised_nnls + >>> x, norm = invert_regularised_nnls(w_matrix, b_vector, tikhonov_matrix=tikhonov_matrix) + """ + + m, n = w_matrix.shape + + if tikhonov_matrix is None: + tikhonov_matrix = eye(n) + + tikhonov_matrix = alpha * tikhonov_matrix + + # Extend W to have form ... + c_matrix = lil((m+n, n)) + c_matrix[0:m, :] = w_matrix[:, :] + c_matrix[m:, :] = tikhonov_matrix[:, :] + c_matrix = c_matrix.tocsr() + + # Extend b to have form ... + d_vector = np.zeros(m+n) + d_vector[0:m] = b_vector[:] + + # Normalise c_matrix and d_vector to avoid possible issues with the nnls termination criteria. + vmax = d_vector.max() + + res = scipy.optimize.lsq_linear(c_matrix / vmax, d_vector / vmax, bounds=(0, np.inf), **kwargs) + + x_vector = res.x + rnorm = np.linalg.norm(res.fun) + + return x_vector, rnorm * vmax From 212484621a6bb2ac96026c72ca6638351280918b Mon Sep 17 00:00:00 2001 From: Jack Lovell Date: Tue, 7 Jul 2026 16:14:16 +0100 Subject: [PATCH 9/9] More improvements to the ADMT demo * Produce additional inversions which incorporate Generomak's tangentially-viewing channels, for comparison with the poloidal-only inversions. * Correctly scale the emissivity. * Fix removing the ray transfer grid from the world before forward-modelling the bolometer measurements. * Better explain why we're building the 1D <-> 2D maps by hand instead of just using the ray transfer matrix capabilities directly: this is for generality. * Use the sparse variant of the NNLS inversion for speed and memory efficiency. The grid is large enough for this to be worthwhile. --- .../bolometry/admt_tomographic_inversion.py | 130 ++++++++++++++++-- 1 file changed, 116 insertions(+), 14 deletions(-) diff --git a/demos/observers/bolometry/admt_tomographic_inversion.py b/demos/observers/bolometry/admt_tomographic_inversion.py index e27fc518..7ef18e47 100644 --- a/demos/observers/bolometry/admt_tomographic_inversion.py +++ b/demos/observers/bolometry/admt_tomographic_inversion.py @@ -20,7 +20,7 @@ from cherab.tools.emitters import RadiationFunction from cherab.tools.raytransfer import RayTransferCylinder, RayTransferPipeline0D from cherab.tools.inversions import admt_utils as admt -from cherab.tools.inversions import invert_regularised_nnls +from cherab.tools.inversions import invert_sparse_regularised_nnls plt.ion() @@ -78,7 +78,8 @@ load_first_wall(world, material=AbsorbingSurface()) bolos = load_bolometers(world) # Only consider the purely-poloidal cameras for now... -bolos = bolos[:3] +poloidal_bolos = bolos[:3] +tangential_bolos = bolos[3:] # Includes midplane and divertor tangential ######################################################################## # Produce a voxel grid @@ -99,7 +100,7 @@ # Build a mask, only including cells within the wall. mask_2d = sample2d_grid(eq.inside_limiter, cell_r, cell_z) mask_3d = mask_2d[:, np.newaxis, :] -ncells = mask_3d.sum() +ncells = int(mask_3d.sum()) # We'll use the Ray Transfer frameworks as these voxels are rectangular # and it's much faster than the Voxel framework for simple cases like this. @@ -119,7 +120,7 @@ ray_transfer_grid.parent = world sensitivity_matrix = [] -for camera in bolos: +for camera in poloidal_bolos: for foil in camera: # Temporarily override foil pipelines for the sensitivity calculation. orig_pipelines = foil.pipelines @@ -136,6 +137,10 @@ foil.pipelines = orig_pipelines sensitivity_matrix = np.asarray(sensitivity_matrix) +# Remove the ray transfer object from the world so it doesn't interfere with +# later observations. +ray_transfer_grid.parent = None + ################################################################################ # Generate the regularisation operators. @@ -143,9 +148,10 @@ print("Generating regularisation operators...") # generate_derivative_operators requires two mappings, one from a flat list of # voxels to the original 2D grid, and one for the 2D grid coordinates to the -# flat list of voxels. We could build these by hand, but the RayTransferCylinder -# object helpfully provides the data already so we just need to convert from -# arrays to dictionaries. +# flat list of voxels. We could build these by hand - and in the general case +# they must be built by hand - but the RayTransferCylinder object we're using +# helpfully provides the data already so we just need to convert from arrays +# to dictionaries. grid_index_1d_to_2d_map = {} for k, idx2d in enumerate(ray_transfer_grid.invert_voxel_map()): # We want the x and z elements, as the Ray Transfer grid is 3D and this @@ -167,8 +173,13 @@ [cell_dx/2, -cell_dz/2]]) # Combine the (N,2) and (4,2) arrays to get an (N,4,2) array. voxel_vertices = voxel_centres[:, None, :] + vertex_displacements[None, :, :] +# The derivative operators are (ncells x ncells) matrices which are sparse. We +# have quite a lot of cells (around 2100), so it's more efficient to generate +# and use sparse matrices here, though dense ones will be returned by default +# for backwards compatibility. +sparse = True derivative_operators = admt.generate_derivative_operators( - voxel_vertices, grid_index_1d_to_2d_map, grid_index_2d_to_1d_map + voxel_vertices, grid_index_1d_to_2d_map, grid_index_2d_to_1d_map, sparse ) # As described in the docstring for generate_derivative_operators, we can @@ -222,7 +233,7 @@ # with all bolometers. The measurements should have the same channel order as # the sensitivity matrix. all_measurements = [] -for camera in bolos: +for camera in poloidal_bolos: all_measurements.extend(camera.observe()) @@ -230,17 +241,19 @@ # Perform the inversions. ################################################################################ print("Performing inversions...") -# We'll use NNLS with regularisation. The hyperparameters have been chosen by +# We'll use NNLS with regularisation. Since the number of voxels is reasonably +# large (around 2100), we'll use the sparse variant of the NNLS inversion for +# memory and computational efficiency. The hyperparameters have been chosen by # hand but techniques such as the discrepancy principle or L curve optimisation # could also be used to determine them. That is out of the scope of this demo. isotropic_alpha = 1e-10 -isotropic_inversion, _ = invert_regularised_nnls( +isotropic_inversion, _ = invert_sparse_regularised_nnls( sensitivity_matrix, all_measurements, alpha=isotropic_alpha, tikhonov_matrix=laplacian, ) admt_alpha = 1e-10 -admt_inversion, _ = invert_regularised_nnls( +admt_inversion, _ = invert_sparse_regularised_nnls( sensitivity_matrix, all_measurements, alpha=admt_alpha, tikhonov_matrix=admt_operator, ) @@ -254,6 +267,94 @@ # Isotropic for index1d, indices2d in grid_index_1d_to_2d_map.items(): emiss2d[indices2d] = isotropic_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Isotropic regularisation,\npoloidal channels") + +# Anisotropic. +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = admt_inversion[index1d] +emiss2d *= 4 * np.pi +plt.figure() +im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') +plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, + levels=np.linspace(0, 1, 10), colors=['k']*9 + ['red']) +plt.xlabel("R[m]") +plt.ylabel("Z[m]") +plt.colorbar(im, label="Inverted\nEmissivity [W/m3]") +plt.xlim([rsamp[0], rsamp[-1]]) +plt.ylim([zsamp[0], zsamp[-1]]) +plt.gca().set_aspect('equal') +plt.title("Anisotropic regularisation\npoloidal channels") + +plt.pause(0.5) + + +######################################################################## +# Can we get a better inversion by including tangential information? +######################################################################## +print("Augmenting the geometry matrix with tangential bolos...") +# The ray transfer object must be in the same world as the bolometers, +# and the plasma emitter must be absent. +ray_transfer_grid.parent = world +emitter.parent = None + + +# sensitivity_matrix = [] +sensitivity_matrix = sensitivity_matrix.tolist() +for camera in tangential_bolos: + for foil in camera: + # Temporarily override foil pipelines for the sensitivity calculation. + orig_pipelines = foil.pipelines + foil.pipelines = [RayTransferPipeline0D(kind=foil.units)] + # All objects in world have wavelength-independent material properties, + # so it doesn't matter which wavelength range we use (as long as + # max_wavelength - min_wavelength = 1) + foil.min_wavelength = 1 + foil.max_wavelength = 2 + foil.spectral_bins = ray_transfer_grid.bins + foil.observe() + sensitivity_matrix.append(foil.pipelines[0].matrix) + # Restore original pipelines for subsequent observe calls. + foil.pipelines = orig_pipelines +sensitivity_matrix = np.asarray(sensitivity_matrix) + +ray_transfer_grid.parent = None + + +print("Adding tangential bolometer measurements...") +emitter.parent = world +for camera in tangential_bolos: + all_measurements.extend(camera.observe()) + + +print("Performing new inversions...") +isotropic_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=isotropic_alpha, + tikhonov_matrix=laplacian, +) + +admt_inversion, _ = invert_sparse_regularised_nnls( + sensitivity_matrix, all_measurements, alpha=admt_alpha, + tikhonov_matrix=admt_operator, +) + +print("Plotting results...") +emiss2d = np.zeros((nx, ny)) + +# Isotropic +for index1d, indices2d in grid_index_1d_to_2d_map.items(): + emiss2d[indices2d] = isotropic_inversion[index1d] +emiss2d *= 4 * np.pi plt.figure() im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, @@ -264,11 +365,12 @@ plt.xlim([rsamp[0], rsamp[-1]]) plt.ylim([zsamp[0], zsamp[-1]]) plt.gca().set_aspect('equal') -plt.title("Isotropic regularisation") +plt.title("Isotropic regularisation,\nall channels") # Anisotropic. for index1d, indices2d in grid_index_1d_to_2d_map.items(): emiss2d[indices2d] = admt_inversion[index1d] +emiss2d *= 4 * np.pi plt.figure() im = plt.imshow(emiss2d.T, extent=(cell_r[0], cell_r[-1], cell_z[0], cell_z[-1]), cmap='Purples') plt.contour(rsamp, zsamp, psisamp.T, linewidths=0.5, alpha=0.3, @@ -279,7 +381,7 @@ plt.xlim([rsamp[0], rsamp[-1]]) plt.ylim([zsamp[0], zsamp[-1]]) plt.gca().set_aspect('equal') -plt.title("Anisotropic regularisation") +plt.title("Anisotropic regularisation\nall channels") plt.ioff() plt.show()