ToFuProject / tofu

Project for an open-source python library for synthetic diagnostics and tomography for Fusion devices
https://tofuproject.github.io/tofu/index.html
MIT License
72 stars 11 forks source link

Tutorial for collimators with multiple apertures #952

Closed dvezinet closed 1 week ago

dvezinet commented 1 month ago

And make error msg more informative !

jlball commented 1 month ago

Error message for aperture list only:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jlball/NCAM/tomography/ncam_tomography/builders.py", line 256, in build_NCAM
    coll.add_diagnostic(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class08_Diagnostic.py", line 90, in add_diagnostic
    dref, ddata, dobj = _check._diagnostics(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 602, in _diagnostics
    ) = _diagnostics_check(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 40, in _diagnostics_check
    doptics, lap, lfilt, lcryst, lgrat = _check_doptics_basics(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 208, in _check_doptics_basics
    _err_doptics(key=key, doptics=doptics)
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 387, in _err_doptics
    raise Exception(msg)
Exception: diag 'ncam': arg doptics must be a dict with:
        - keys: key to existing camera
        - values: existing optics (aperture, filter, crystal)
                - as a list of str for regular cameras
                - as a list of (tuples of str) for collimator cameras

Provided:
{'cam0': ('cam0_ap00', 'cam0_ap01', 'cam0_ap10', 'cam0_ap11', 'cam0_ap20', 'cam0_ap21', 'cam0_ap30', 'cam0_ap31', 'cam0_ap40', 'cam0_ap41', 'cam0_ap50', 'cam0_ap51', 'cam0_ap60', 'cam0_ap61', 'cam0_ap70', 'cam0_ap71', 'cam0_ap80', 'cam0_ap81', 'cam0_ap90', 'cam0_ap91', 'cam0_ap100', 'cam0_ap101', 'cam0_ap110', 'cam0_ap111', 'cam0_ap120', 'cam0_ap121', 'cam0_ap130', 'cam0_ap131', 'cam0_ap140', 'cam0_ap141', 'cam0_ap150', 'cam0_ap151', 'cam0_ap160', 'cam0_ap161', 'cam0_ap170', 'cam0_ap171', 'cam0_ap180', 'cam0_ap181')}

Error message for aperture list and path matrix:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/jlball/NCAM/tomography/ncam_tomography/builders.py", line 257, in build_NCAM
    coll.add_diagnostic(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class08_Diagnostic.py", line 90, in add_diagnostic
    dref, ddata, dobj = _check._diagnostics(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 602, in _diagnostics
    ) = _diagnostics_check(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 40, in _diagnostics_check
    doptics, lap, lfilt, lcryst, lgrat = _check_doptics_basics(
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 208, in _check_doptics_basics
    _err_doptics(key=key, doptics=doptics)
  File "/home/jlball/.local/lib/python3.8/site-packages/tofu/data/_class8_check.py", line 387, in _err_doptics
    raise Exception(msg)
Exception: diag 'ncam': arg doptics must be a dict with:
        - keys: key to existing camera
        - values: existing optics (aperture, filter, crystal)
                - as a list of str for regular cameras
                - as a list of (tuples of str) for collimator cameras

Provided:
{'cam0': {'optics': ['cam0_ap00', 'cam0_ap01', 'cam0_ap20', 'cam0_ap21', 'cam0_ap40', 'cam0_ap41', 'cam0_ap60', 'cam0_ap61', 'cam0_ap80', 'cam0_ap81', 'cam0_ap100', 'cam0_ap101', 'cam0_ap120', 'cam0_ap121', 'cam0_ap140', 'cam0_ap141', 'cam0_ap160', 'cam0_ap161', 'cam0_ap180', 'cam0_ap181'], 'paths': array([[ True,  True, False, False, False, False, False, False, False,
        False, False, False, False, False, False, False, False, False,
        False, False],
       [False, False,  True,  True, False, False, False, False, False,
        False, False, False, False, False, False, False, False, False,
        False, False],
       [False, False, False, False,  True,  True, False, False, False,
        False, False, False, False, False, False, False, False, False,
        False, False],
       [False, False, False, False, False, False,  True,  True, False,
        False, False, False, False, False, False, False, False, False,
        False, False],
       [False, False, False, False, False, False, False, False,  True,
         True, False, False, False, False, False, False, False, False,
        False, False],
       [False, False, False, False, False, False, False, False, False,
        False,  True,  True, False, False, False, False, False, False,
        False, False],
       [False, False, False, False, False, False, False, False, False,
        False, False, False,  True,  True, False, False, False, False,
        False, False],
       [False, False, False, False, False, False, False, False, False,
        False, False, False, False, False,  True,  True, False, False,
        False, False],
       [False, False, False, False, False, False, False, False, False,
        False, False, False, False, False, False, False,  True,  True,
        False, False],
       [False, False, False, False, False, False, False, False, False,
        False, False, False, False, False, False, False, False, False,
         True,  True]])}}
jlball commented 1 month ago

Code used to generate the camera and LOS: Most of the activity is wrapped up into the build_NCAM() function:

# make sure to have an interactive backend
import sys
import os
import numpy as np
import matplotlib.pyplot as plt
from transp_profiles_2d import get_transp_profiles_2d
import bsplines2d as bs2d
from scipy.interpolate import RegularGridInterpolator
plt.ion()
#get_ipython().run_line_magic('matplotlib', 'qt5')
import tofu as tf
from pandas import read_csv

#######################################################
# Line of sight configurations
#######################################################

LOS_config_dict = {}

# All 19 lines of sight to be installed on SPARC
LOS_config_dict["all"] = [i for i in range(0, 19)]

# The 10 LOS currently being considered for NCAM
LOS_config_dict["10_los"] = [0, 2, 4, 6, 8, 10, 12, 14, 16, 18]

# 10 LOS NCAM geometry but with up / down asymmetry 
LOS_config_dict["10_los_asymmetric"] = [0, 2, 4, 6, 8, 9, 11, 13, 15, 17]
LOS_config_dict["7_los_asymmetric"] = [4, 6, 8, 9, 11, 13, 15]

LOS_config_dict["midplane_only"] = [7, 8, 9, 10, 11]

def build_NCAM(LOS_config="all", 
               transp_file="12345Q90.CDF", 
               geqdsk_file="SPARC_V1E_transp_3.geq", 
               config="SPARC-V0", 
               collimator_diameters=[1], 
               sensor_collimator_dist=10,
               sensor_side_length=3,
               sensor_shape = "square", 
               angle_tolerance=0,
               endpoint_tolerance=0):
    #######################################################
    # Configure Tofu
    #######################################################
    # simpler config => faster computation and plotting
    conf = tf.load_config(config)

    # instanciate a Diagnostic
    coll = tf.data.Collection()

    #Get emmissivity profile and its position grids
    r_grid = np.linspace(1.85 - 0.6, 1.85 + 0.6, num=100)
    z_grid = np.linspace(-0.6*1.8, 0.6*1.8, num=200)

    z_mesh, r_mesh = np.meshgrid(z_grid, r_grid)

    data_dict = get_transp_profiles_2d(5.973, f'data/{transp_file}', f'data/{geqdsk_file}', r_grid, z_grid)

    # Populate with a rectangular mesh
    coll.add_mesh_2d_rect(
        # naming
        key='m0',
        # parameters
        domain=None,           # if res is provided, assumes a uniform mesh, for rectangular use res = [resR, resZ]
        knots0 = r_grid,         # if res not provided, can provide the R vector of knots directly
        knots1 = z_grid,         # if res not provided, can provide the Z vector of knots directly
        # optional cropping
        crop_poly=conf,     # cropping of the rectangular mesh with the vessel in the Config
        units=['m', 'm'],
    )

    # add a 2d bsplines basis on that mesh
    coll.add_bsplines(key='m0', deg=1)

    # make it time dependent
    nt0 = 2
    t0 = np.linspace(0, 1, nt0)

    # 2d radius
    radius2d = np.cos(0 * t0[:, None, None])*data_dict["sqrt_psiN_oi"]

    coll.add_ref(key='nt0', size=nt0)
    coll.add_data(key='t0', data=t0, ref='nt0', dim='time')
    coll.add_data(key='radius2d', data=radius2d, ref=('nt0','m0_bs1'))

    num_LOS = len(LOS_config_dict[LOS_config])

    # add polar mesh based on radius2d, and associated radial bsplines of degree 2
    coll.add_mesh_1d(
        key='m1',
        knots=np.linspace(0, 1, num_LOS-2),
        subkey='radius2d',
        deg=2,
    )

    # extract time-dependent radius
    t0 = coll.ddata['t0']['data']
    key_radius2d = coll.dobj['mesh']['m1']['subkey'][0]
    radius2d = coll.ddata[key_radius2d]['data']

    # define emissivity
    emiss2d = np.cos(0 * t0[:, None, None])*data_dict["YDTn_oi_n_m3_s"] 

    # add emissivity
    coll.add_data(key='emiss2d', 
                data=emiss2d/(4*np.pi), 
                ref=('nt0', 'm0_bs1'), 
                units='ph/(m3.sr)')

    #######################################################
    # Create lines of sight
    #######################################################

    LOS_data = read_csv("data/CALC-NTRN-NC0-NOMINAL-LOS.csv")
    z_points = LOS_data["Z_DH (m)"]
    angles = -LOS_data["Angle (deg)"]
    x_points = LOS_data["R_DH (m)"]

    # Perturb collimator translation and angle based on normal distribution:
    angle_perturbs = np.random.normal(0, scale=angle_tolerance, size=19)
    z_perturbs = np.random.normal(0, scale=endpoint_tolerance, size=19)

    angles = angles + angle_perturbs
    z_points = z_points + z_perturbs

    # prepare theta for describing the circular outline of apertures
    theta = np.pi * np.linspace(-1, 1, 50)

    collimator_length = 2.863 # m
    sensor_collimator_dist = sensor_collimator_dist/100

    #doptics["cam0"] = {"optics":[], "paths":np.full((num_LOS, 2*num_LOS), False)}
    doptics = {'cam0': []}

    for i, idx in enumerate(LOS_config_dict[LOS_config]):
        # prepare radius of collimator
        if len(collimator_diameters) == 1:
            radius = collimator_diameters[0]/200 # convert from cm to m, diameter to radius
        else:
            if len(collimator_diameters) == len(LOS_config_dict[LOS_config]):
                radius = collimator_diameters[i]/200
            else:
                raise ValueError("number of collimator diameters is different from number of collimators")

        # set keys of camera, aperture 0 and 1 (2 ends of the collimator)
        #kcam = f"ncam_c{idx}"
        kap0 = f'cam0_ap{idx}0'
        kap1 = f'cam0_ap{idx}1'

        # los unit vector
        nin = np.r_[-np.cos(angles[idx]*(np.pi/180)), 0, np.sin(angles[idx]*(np.pi/180))]
        e0 = np.r_[nin[1], nin[0], 0.]
        e0 = e0 / np.linalg.norm(e0)
        e1 = np.cross(nin, e0)

        # center of collimator end
        coll_end_cent = np.r_[x_points[idx], 0, z_points[idx]]
        sensor_cent = coll_end_cent - nin * sensor_collimator_dist

        # geometry dict for single-pixel camera 1d
        if sensor_shape == "square":
            dgeom = {
                'cents_x': np.r_[sensor_cent[0]],
                'cents_y': np.r_[sensor_cent[1]],
                'cents_z': np.r_[sensor_cent[2]],
                'nin_x': nin[0],
                'nin_y': nin[1],
                'nin_z': nin[2],
                'e0_x': e0[0],
                'e0_y': e0[1],
                'e0_z': e0[2],
                'e1_x': e1[0],
                'e1_y': e1[1],
                'e1_z': e1[2],
                'outline_x0': sensor_side_length/200 * np.r_[-1, 1, 1, -1],
                'outline_x1': sensor_side_length/200 * np.r_[-1, -1, 1, 1],
            }

        elif sensor_shape == "circle":
            thetas = np.linspace(0, 2*np.pi, num=50)
            dgeom = {
                'cents_x': np.r_[sensor_cent[0]],
                'cents_y': np.r_[sensor_cent[1]],
                'cents_z': np.r_[sensor_cent[2]],
                'nin_x': nin[0],
                'nin_y': nin[1],
                'nin_z': nin[2],
                'e0_x': e0[0],
                'e0_y': e0[1],
                'e0_z': e0[2],
                'e1_x': e1[0],
                'e1_y': e1[1],
                'e1_z': e1[2],
                'outline_x0': sensor_side_length/200 * np.sin(thetas),
                'outline_x1': sensor_side_length/200 * np.cos(thetas),
            }

        else:
            raise ValueError("Invalid sensor shape specified")

        # add apertures

        dgeom = {
            'cent': coll_end_cent,
            'nin': nin,
            'e0': e0,
            'e1': e1,
            'outline_x0': radius * np.cos(theta),
            'outline_x1': radius * np.sin(theta),
        }

        coll.add_aperture(key=kap0, **dict(dgeom))

        # update center position
        dgeom['cent'] = coll_end_cent + collimator_length * nin
        coll.add_aperture(key=kap1, **dict(dgeom))

        # fill in doptics
        #doptics["cam0"].append((kap1, kap0))
        # doptics["cam0"]["paths"][i, 2*i] = True
        # doptics["cam0"]["paths"][i, 2*i + 1] = True
        doptics["cam0"] += [kap0, kap1]

    # ---------------------------------------------------------------
    # create the diagnostic ncam1 using the info contained in doptics

    # transforming to tuple => collimator camera
    doptics['cam0'] = tuple(doptics['cam0'])
    #print(doptics)

    coll.add_diagnostic(
        key='ncam',
        doptics=doptics,
        compute=True,
        config=conf,
    )   

    return coll
dvezinet commented 1 month ago

OK, I see where the issue is, it's not a bug in tofu:

You're adding apertures and combining them into a diagnostic

But you forgot to add the camera too.

Indeed, the doptics dict should have a key that's refereing to an existing camera, but here 'cam0' refers to nothing as it was never created.

I would suggest that before your loop on lines of sight, you initialize a dgeom dict that will hold the geometry of the camera (the camera being a collection of pixels).

Inside the loop, for each LOS, you append the content of dgeom to declare each pixel In your case all cameras are co-planar, so you only need a unique set of unit vectors (nin, e0, e1)

dgeom = {
    'outline_x0': outline of one pixel (in its own plane), first coordinate,
    'outline_x1': outline of one pixel (in its own plane), second coordinate,
    'cents_x': array of x coordinates of centers of each pixel,
    'cents_y': array of y coordinates of centers of each pixel,
    'cents_z': array of z coordinates of centers of each pixel,
    'nin': unit vector (x, y, z), normalized, normal to the pixel surface, oriented towards the plasma,
    'e0': unit vector (x, y, z), normalized, in which outline_x0 was expressed,
    'e1': unit vector (x, y, z), normalized, such that e2 = np.cross(nin, e1),
}

Then , after the loop, but before you call add_diagnostic(), you call add_camera_1d(key='cam0', dgeom=dgeom) to add the camera so 'cam0' actually corresponds to an existing camera when tofu inspects the doptics dict.

Although that did not correspond to a bug in tofu, i will use this feedback to improve the error message and make it more informative. Let me know how it goes and if you need more help,