Get started on AUG

[Didou09]

Hi @fsciortino, hopefully this brief tutorial will help you get started fast:

In tofu, a tokamak geometry is called a Configuration (Config class), to account for the fact that a tokamak with the same name can exist in various configurations (with a limiter, then with a divertor, then with a new antenna...).

To create arbitrary cameras on AUG, start by loading the AUG geometry using the tf.load_config() routine.

You can get a list of the available tokamak geometries and their different versions by typing using the routine without giving arguments, which will raise an error printing the list of choices to pick from:

In [1]: import tofu as tf

In [2]: tf.load_config()
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
<ipython-input-2-b68e02c00c2c> in <module>
----> 1 tf.load_config()

~/ToFu_All/tofu_git/tofu/tofu/geom/utils.py in create_config(case, Exp, Type, Lim, Bump_posextent, R, r, elong, Dshape, divlow, divup, nP, returnas, strict, SavePath, path)
    987     elif not any(lc):
    988         msg = get_available_config(verb=False, returnas=str)
--> 989         raise Exception(msg)
    990 
    991     # Get config, either from known case or geometrical parameterization

Exception: A config is the geometry of a tokamak
You can define your own, see online tutorial at:
    https://tofuproject.github.io/tofu/auto_examples/tutorials/tuto_plot_create_geometry.html
tofu also also provides some pre-defined config ready to load
They are available via their name or via shortcuts

    unique names    shortcuts
    ------------    -------------------------
    - AUG-V1        ['AUG-V1', 'AUG']
    - COMPASS-V0    ['COMPASS-V0', 'COMPASS']
    - COMPASS-V1    ['COMPASS-V1']
    - COMPASS2-V0   ['COMPASS2-V0']
    - DEMO-2019     ['DEMO-2019', 'DEMO']
    - ITER-V0       ['ITER-V0']
    - ITER-V1       ['ITER-V1', 'A2']
    - ITER-V2       ['ITER-V2', 'B4', 'ITER']
    - JET-V0        ['JET-V0', 'JET']
    - NSTX-V0       ['NSTX-V0', 'NSTX']
    - SPARC-V0      ['SPARC-V0']
    - SPARC-V1      ['SPARC-V1', 'SPARC']
    - SPARC-V2      ['SPARC-V2']
    - TCV-V0        ['TCV-V0', 'TCV']
    - TOMAS-V0      ['TOMAS-V0', 'TOMAS']
    - WEST-Sep      ['WEST-Sep', 'A3']
    - WEST-V0       ['WEST-V0']
    - WEST-V1       ['WEST-V1', 'A1']
    - WEST-V2       ['WEST-V2', 'B1']
    - WEST-V3       ['WEST-V3', 'B2']
    - WEST-V4       ['WEST-V4', 'B3', 'WEST']

  => to get a pre-defined config, call for example:
    config = tf.geom.utils.create_config('ITER')

 => to also load coils add '-coils' to the config name, e.g.:
    config = tf.geom.utils.create_config('ITER-coils')

Pick a low-version (simpler) if you want faster computations, a high-version if you want details. The routine returns a Config instance, just type it and press enter to get a list of its content (in terms of structures, nb of points in each (R, Z) polygon, number of toroidal occurences (0 => axisymmetric)...

In [3]: conf = tf.load_config('AUG')

In [4]: conf
Out[4]: 
<bound method Config.plot of tot. Struct  tot. occur  tot. points
-----------  ----------  -----------
33           33          595        

class  Name     SaveName                                                nP  noccur  move  color                     visible
-----  -------  ------------------------------------------------------  --  ------  ----  ------------------------  -------
Ves    VESiR    TFG_Ves_ExpAUG_VESiR_sh00000_Vers1.5.0-142-g1ccea1b8    36  0       None  ( 0.0,  0.0,  0.0,  1.0)  True   
PFC    D2cdome  TFG_PFC_ExpAUG_D2cdome_sh00000_Vers1.5.0-142-g1ccea1b8  14  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2cdomL  TFG_PFC_ExpAUG_D2cdomL_sh00000_Vers1.5.0-142-g1ccea1b8  14  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2cdomR  TFG_PFC_ExpAUG_D2cdomR_sh00000_Vers1.5.0-142-g1ccea1b8  14  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2ci1    TFG_PFC_ExpAUG_D2ci1_sh00000_Vers1.5.0-142-g1ccea1b8    12  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2ci2    TFG_PFC_ExpAUG_D2ci2_sh00000_Vers1.5.0-142-g1ccea1b8    10  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2cTPib  TFG_PFC_ExpAUG_D2cTPib_sh00000_Vers1.5.0-142-g1ccea1b8  23  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2cTPic  TFG_PFC_ExpAUG_D2cTPic_sh00000_Vers1.5.0-142-g1ccea1b8  12  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2cTPi   TFG_PFC_ExpAUG_D2cTPi_sh00000_Vers1.5.0-142-g1ccea1b8   18  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBG2   TFG_PFC_ExpAUG_D2dBG2_sh00000_Vers1.5.0-142-g1ccea1b8   28  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBl1   TFG_PFC_ExpAUG_D2dBl1_sh00000_Vers1.5.0-142-g1ccea1b8   13  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBl2   TFG_PFC_ExpAUG_D2dBl2_sh00000_Vers1.5.0-142-g1ccea1b8   14  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBl3   TFG_PFC_ExpAUG_D2dBl3_sh00000_Vers1.5.0-142-g1ccea1b8   13  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBu1   TFG_PFC_ExpAUG_D2dBu1_sh00000_Vers1.5.0-142-g1ccea1b8   22  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBu2   TFG_PFC_ExpAUG_D2dBu2_sh00000_Vers1.5.0-142-g1ccea1b8   11  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBu3   TFG_PFC_ExpAUG_D2dBu3_sh00000_Vers1.5.0-142-g1ccea1b8   11  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D2dBu4   TFG_PFC_ExpAUG_D2dBu4_sh00000_Vers1.5.0-142-g1ccea1b8   8   0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D3BG10   TFG_PFC_ExpAUG_D3BG10_sh00000_Vers1.5.0-142-g1ccea1b8   19  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    D3BG1    TFG_PFC_ExpAUG_D3BG1_sh00000_Vers1.5.0-142-g1ccea1b8    20  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    ICRHa    TFG_PFC_ExpAUG_ICRHa_sh00000_Vers1.5.0-142-g1ccea1b8    61  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    LIM09    TFG_PFC_ExpAUG_LIM09_sh00000_Vers1.5.0-142-g1ccea1b8    8   0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    PClow    TFG_PFC_ExpAUG_PClow_sh00000_Vers1.5.0-142-g1ccea1b8    6   0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    PCup     TFG_PFC_ExpAUG_PCup_sh00000_Vers1.5.0-142-g1ccea1b8     6   0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    SBi      TFG_PFC_ExpAUG_SBi_sh00000_Vers1.5.0-142-g1ccea1b8      42  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPLT1    TFG_PFC_ExpAUG_TPLT1_sh00000_Vers1.5.0-142-g1ccea1b8    17  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPLT2    TFG_PFC_ExpAUG_TPLT2_sh00000_Vers1.5.0-142-g1ccea1b8    12  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPLT3    TFG_PFC_ExpAUG_TPLT3_sh00000_Vers1.5.0-142-g1ccea1b8    15  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPLT4    TFG_PFC_ExpAUG_TPLT4_sh00000_Vers1.5.0-142-g1ccea1b8    23  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPLT5    TFG_PFC_ExpAUG_TPLT5_sh00000_Vers1.5.0-142-g1ccea1b8    22  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPRT2    TFG_PFC_ExpAUG_TPRT2_sh00000_Vers1.5.0-142-g1ccea1b8    29  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPRT3    TFG_PFC_ExpAUG_TPRT3_sh00000_Vers1.5.0-142-g1ccea1b8    14  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPRT4    TFG_PFC_ExpAUG_TPRT4_sh00000_Vers1.5.0-142-g1ccea1b8    12  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   
PFC    TPRT5    TFG_PFC_ExpAUG_TPRT5_sh00000_Vers1.5.0-142-g1ccea1b8    16  0       None  ( 0.8,  0.8,  0.8,  0.8)  True   >

In [5]: lax = conf.plot()   # will return a list of 2 matplotlib axes handles

You can add random colors to cheer it up:

In [6]: conf.set_colors_random()

In [7]: lax = conf.plot()

You then easily create a Mesh2D instance based on this geometry, first create the instance (a generic data holder class that can handle several meshes), then add a rectangular mesh using the configuration as a cropping limit, with resolution 10 cm, and also add bsplines of degree 1 to it:

 In [8]: mesh = tf.data.Mesh2D()

In [9]: mesh.add_mesh(crop_poly=conf, res=0.1, deg=1, key='mesh1')  # give it a unique name (key)

Again, type the variable and pressing enter will tell what's inside, and there is a plotting method (that will return a dict of axes, this is currently being refactored for harmonization)

In [10]: mesh
Out[10]: 
group  nb. ref  nb. data
-----  -------  --------
R      2        3       
Z      2        3       

ref key        group  size  nb. data
-------------  -----  ----  --------
mesh1-R-knots  R      17    1       
mesh1-Z-knots  Z      28    1       
mesh1-R-cents  R      16    2       
mesh1-Z-cents  Z      27    2       

data key       shape     ref                                 group       units  dim       quant  name     monot           source 
-------------  --------  ----------------------------------  ----------  -----  --------  -----  -------  --------------  -------
mesh1-R-knots  (17,)     ('mesh1-R-knots',)                  ('R',)      m      distance  R      R        (True,)         unknown
mesh1-Z-knots  (28,)     ('mesh1-Z-knots',)                  ('Z',)      m      distance  Z      Z        (True,)         unknown
mesh1-R-cents  (16,)     ('mesh1-R-cents',)                  ('R',)      m      distance  R      R        (True,)         unknown
mesh1-Z-cents  (27,)     ('mesh1-Z-cents',)                  ('Z',)      m      distance  Z      Z        (True,)         unknown
mesh1-crop     (16, 27)  ('mesh1-R-cents', 'mesh1-Z-cents')  ('R', 'Z')  a.u.   bool      bool   unknown  (False, False)  unknown

mesh   type  knots                               cents                               ref                                 shape     variable  crop        crop-thresh
-----  ----  ----------------------------------  ----------------------------------  ----------------------------------  --------  --------  ----------  -----------
mesh1  rect  ('mesh1-R-knots', 'mesh1-Z-knots')  ('mesh1-R-cents', 'mesh1-Z-cents')  ('mesh1-R-cents', 'mesh1-Z-cents')  (16, 27)  False     mesh1-crop  3          

In [11]: dax =  mesh.plot_mesh()

You can also plot any individual 2D bspline providing its index in the grid:

In [12]: dax = mesh.plot_bsplines(ind=(10, 10))

Now you can add two 1d pinhole cameras as two set of lines of sight using the geometry utility:

In [13]: cam0 = tf.geom.utils.create_CamLOS1D(pinhole=[2.4, 0, 0], orientation=[np.pi, 0, 0], focal=0.1, sensor_nb=100, sensor_size=0.1, config=conf, Name='cam1', Diag='SXR')

In [14]: cam1 = tf.geom.utils.create_CamLOS1D(pinhole=[1.15, 0, 0], orientation=[-np.pi/4, 0, 0], focal=0.1, sensor_nb=100, sensor_size=0.1, config=conf, Name='cam2', Diag='SXR')

Feel free to plot in astatic or interactive figure each camera, or to concatenate them into a single total camera: The plot_touch() method provides a left-clickable interactive figure to let you select and LOS and get it highlighted in the cross-section and horizontal plot. By default, the plot shows the length of the LOS, and the color of each LOS is inherited from the structural element on which it ends. You can also plot the incidence angle of the LOS on sais structural element if you want, using quant='angles'

In [15]: lax = cam0.plot()

In [16]: cam1.plot_touch()

In [17]: camtot = cam0 + cam1   # concatenates cameras together

In [18]: camtot.plot_touch(quant='angles')

You can then decide to compute the geometry matrix of the mesh you created using the total camera you created too (add the geometry matrix to the generic data handler mesh), with a spatial resolution of mm for the integrals. Then plot to explore the content of the geometry matrix

In [19]: mesh.add_geometry_matrix(cam=camtot, res=0.005)
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/data/_mesh_comp.py:193: UserWarning: ind is not 2d anymore!
  warnings.warn("ind is not 2d anymore!")
Geom. matrix, chan 200 / 200   

In [20]: dax = mesh.plot_geometry_matrix()
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/data/_mesh_comp.py:193: UserWarning: ind is not 2d anymore!
  warnings.warn("ind is not 2d anymore!")

[fsciortino]

thanks @Didou09 for this nice overview! While I play with these methods, two first questions come to mind:

1. How can I define a 1D LOS simply based on start and end points? It seems that create_CamLOS1D allows for LOS definitions using a pinhole and 3 angles. Is that what you would recommend one to use in this case?

2. Obvious follow-up: once I define LOS geometries and pre-computed the geometry matrix, how can I provide some measurements and run the tomography?

It isn't so clear to me what is being plotted in the last figure above (obtained by running mesh.plot_geometry_matrix()). I'm looking more into it now, but it would be helpful if you could clarify here too, please.

[Didou09]

@fsciortino , can we discuss now ?

[Didou09]

In a nutshell:

1. Defining rays from start and end points

# start_pts = (3, N) array of floats containing the (X, Y, Z) coordinates of points
# end_pts = (3, N) array of floats containing the (X, Y, Z) coordinates of points

# Derive the unit vectors for each LOS
vect = end_pts - start_pts
vect = vect / sqrt(np.sum(vect**2, axis=0))

# Provide start_pts and vect as a tuple under argument dgeom (dictionary of geometric paramaters) to CamLOS1D
# Give it a Name, Diag, Exp
# Don't forget to give the configuration too

conf = tf.load_config('AUG')
cam = tf.geom.CamLOS1D(dgeom=(start_pts, vect), config=conf, Name='Cam1', Diag='SXR', Exp='AUG')

[fsciortino]

@Didou09 I'm running into the problem that my LoS are defined as starting from a point below the divertor structure and therefore the mesh.add_mesh(crop_poly=conf, res=0.1, deg=1, key='mesh1') crops the LoS in the wrong way -- see the image below.

I need to correct the cropping both for the incoming and outgoing part of each ray. I could try to crop this myself, but there's probably already something in ToFu to deal with this properly...?

[Didou09]

There are at least 2 ways to do this:

1/ Once the camera is created, it holds the tokamak configuration as an attribute, you can modify that attribute to tell the camera that it should ignore one particular structure element when it computes the ray-tracing, and force the ray-tracing to be re-computed

Load the configuration and create a camera (either a pinhole camera using tf.geom.utils.create_CamLOS1D() or directly by specifying start_point and vect using tf.geom.CamLOS1D(), both will be united in a single routine in future versions):

In [1]: conf = tf.load_config('AUG')

In [2]: config.set_colors_random()   # assign random colors to each element to spot easily the name of the one you're interested in

n [3]: config.plot()

Create the camera, observe that it is not as desired, then deactivate the PFC named D2cdomR in the compuitation, and re-run the ray-tracing computation:

In [4]: cam = tf.geom.utils.create_CamLOS1D(pinhole=[1.48, 0, -1.15], orientation=[0,0,0], sensor_nb=10, sensor_size=0.05, focal=0.1, config=conf, Name='test', Diag='test', Exp='AUG')

In [5]: cam.plot()

In [6]: cam.config    # config is now an attribute of cam, see what's inside, we want to set compute=False for D2cdomR

In [7]: cam.config.PFC.D2cdomR.set_compute(False)

In [7]: cam.config.PFC.D2cdomR.set_compute(False)

In [8]: cam.compute_dgeom()  # re-compute dict of geometry parameters

In [9]: cam.plot()      # double check

[Didou09]

Then to compute a synthetic signal from an arbitrary emissivity function: 1/ define an emissivity function, it should take at least pts (3, N) X, Y, Z point coordinates array as argument, and optionally t (time vector).

• it t is not provided, the function shall return a (N,) array (one value per pts)
• If t is provided, it shall return a (nt. N) array 2/ Feed it to cam.calc_signal(), along with the desired spatial resolution res (m) for LOS discretization, and optionally a time vector 3/ It will give the synthetic signal as a sig array (nt, nchan), a units str, and optionally plot the result in an interactive figure

In [57]: def emiss(pts, t=None):
    ...:     R = np.hypot(pts[0, :], pts[1, :])
    ...:     Z = pts[2, :]
    ...:     e0 = np.exp(-(R-1.575)**2/0.1**2 - (Z+1.15)**2/0.01**2)
    ...:     if t is not None:
    ...:         e0 = e0[None, :] * (1. + 0.1*np.cos(t[:, None]))
    ...:     return e0
    ...: 

n [58]: sig, unit = cam.calc_signal(emiss, res=0.001, t=np.linspace(0,10,100))

(clic on axes to select time / channel)

[fsciortino]

@Didou09 thanks for these extra steps. This all works great! So, I suppose that the final steps to do a tomography are:

1. add measured signals for each of my 1D LOS cameras
2. add uncertainties for each signal above...?
3. run a minimization that finds the bspline coefficients which, when multiplying each spline, reproduce a local emissivity map that reproduces all the line-integrated signals that I provided.

How can I do these final steps? Also, any advice on how to assess the quality of the tomographic reconstruction? (e.g. standard plots to look at a \chi^2, or plots that compare residua, etc.) Thanks!

[Didou09]

Hi @fsciortino,

For your specific problem, I've made a few minor upgrade to tofu to make your life a bit easier. They are already available on the devel branch and will be officially release in the next version (1.5.1).

2 options:

• If you are in a hurry and can't wait for official releases, just clone the devel branch to always get the last upgrades (read-only)
• If you prefer a clean pip / conda install and are in no hurry, wait for the release and upgrade when available (probably next week)

In both cases, you can now:

• load a simpler version of the AUG geometry ('AUG-V0'), boiling down to a single closed polygon representing the first wall. This saves you the need to de-activate some structural elements in the ray-tracing computation
• Create a custom mesh by specifying both limits in (R, Z) via domain and a polygon for cropping (you had to choose between these 2 options before PR #606)

To load a simple version of AUG and create a camera:

In [1]: import tofu as tf

In [2]: conf = tf.load_config('AUG-V0')

In [3]: cam0 = tf.geom.utils.create_CamLOS1D(pinhole=[1.48, 0, -1.15], orientation=[0,0,0], sensor_nb=10, sensor_size=0.01, focal=0.02, config=conf, Name='test', Diag='test', Exp='AUG')

In [4]: cam2 = tf.geom.utils.create_CamLOS1D(pinhole=[1.46, 0, -1.10], orientation=[0,0,0], sensor_nb=10, sensor_size=0.01, focal=0.02, config=conf, Name='test', Diag='test', Exp='AUG')

In [5]: cam = cam0 + cam1 + cam2

In [6]: cam.plot()

Then create a mesh with the relevant limits and desired bsplines degree, then compte the geometry matrix using the concatenated camera and plot it, choosing a basis function and channel indew such that you can better understand the plots (and add the AUG geometry on top):

In [7]: mesh = tf.data.Mesh2D()

In [8]: mesh.add_mesh(domain=[[1.45, 1.65], [-1.22, -1.05]], crop_poly=conf, res=0.008, deg=0, key='mesh1')

In [9]: mesh.add_geometry_matrix(cam=cam, res=0.001)
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/data/_mesh_comp.py:193: UserWarning: ind is not 2d anymore!
  warnings.warn("ind is not 2d anymore!")
Geom. matrix, chan 30 / 30   

In [10]: dax = mesh.plot_geometry_matrix(indbf=50, indchan=15, cam=cam)
/Home/DV226270/ToFu_All/tofu_git/tofu/tofu/data/_mesh_comp.py:193: UserWarning: ind is not 2d anymore!
  warnings.warn("ind is not 2d anymore!")

In [11]: conf.plot(proj='cross', lax=dax['crosstot']['ax'])
Out[11]: <AxesSubplot:xlabel='Z (m)', ylabel='R (m)'>

Then define a custom emissivity function, use it to compute synthetic signal:

In [12]: def emiss(pts, t=None):
    ...:     R = np.hypot(pts[0, :], pts[1, :])
    ...:     Z = pts[2, :]
    ...:     e0 = np.exp(-(R-1.575)**2/0.1**2 - (Z+1.15)**2/0.01**2)
    ...:     if t is not None:
    ...:         e0 = e0[None, :] * (1. + 0.1*np.cos(t[:, None]))
    ...:     return e0

In [13]: time = np.linspace(0, 10, 11)

In [14]: sig, unit = cam.calc_signal(emiss, res=0.001, t=time)

Now you can compute an inversion.

An inversion requires several intermediates (like the regularity operator...). As a preliminary to computing the inversion, we need to store the time vector we used and the synthetic data, it takes 2 extra lines, using a tofu-specific syntax:

In [15]: mesh.add_ref(key='time', data=time, group='time')

In [16]: mesh.add_data(key='synthetic data', data=sig.data, ref=('time', 'chan0'))

You can then check that mesh indeed contains the time vector and data array, which means we will now be able to access them by key:

You can now add an inversion, using this data, to the mesh instance:

In [16]: mesh.add_inversion(key_data='synthetic data', sigma=None, operator='D1N2', verb=1)
    time step 1 / 11    chi2n = 5.417e-01    niter = 7
    time step 2 / 11    chi2n = 5.031e-01    niter = 3
    time step 3 / 11    chi2n = 4.270e-01    niter = 4
    time step 4 / 11    chi2n = 3.842e-01    niter = 4
    time step 5 / 11    chi2n = 4.090e-01    niter = 4
    time step 6 / 11    chi2n = 4.821e-01    niter = 4
    time step 7 / 11    chi2n = 5.384e-01    niter = 3
    time step 8 / 11    chi2n = 5.209e-01    niter = 3
    time step 9 / 11    chi2n = 4.480e-01    niter = 4
    time step 10 / 11    chi2n = 3.899e-01    niter = 4
    time step 11 / 11    chi2n = 3.952e-01    niter = 3
0.17695121094584465 s
Post-formatting results...

And plot the result (select the time step to plot using indt=...):

In [79]: mesh.plot_inversion(indt=0)

On the right hand side of the plot:

• original data (dashed line) and reconstructed data (solid line), versus channel index
• error, versus channel index
• value of the potential the algorithm minimizes, vs time
• nb of iterations for convergence (vs tine)

Questions you should keep in mind when performing an inversion:

• Do you have a good geometrical coverage (i.e.: your cameras, are they wisely placed ?)
• Do you have a relevant mesh resolution ? Should you use a finer mesh ? a uniform or a variable mesh size ?
• Are you using the most relevant basis functions for your problem ? (here: the degree of the bsplines)
• Did you choose a relevant regularity operator ? (here D1N2 = minimize the squared norm of the gradient this table for a larger choice)
• Is the criterion for convergence relevant ? (by default conv_crit=1.e-4, but you can change it)

They are other parameters you can tweak, but these are probably the most important points.

[Didou09]

Tips:

• When in doubt because your inversion looks very ugly, always check the geometry matrix (ex: look for pixels not detected by any line of sight)
• If you're trying to reconstruct an emissivity field that is strong close to the mesh edge, prefer pixels (bsplines of degree 0), because higher-order bsplines tend to smooth the solution close to the edges (an effect that can be circumvented at the cost of a finer mesh).

[fsciortino]

@Didou09 I've made some progress using some synthetic data from SOLPS simulations, but the result seems wrong. In the figure below, I'm plotting the reconstruction on the left and the original signal on the right. I've overplotted the 2D simulation grid on top of the tomography result just for ease of comparison.

I've used bsplines of order 0 to start with. When trying to use order 1, the inversion gave me an assertion error that I couldn't decipher...

I set up the problem as time-independent and checked that the lines of sights as given to tofu are correctly representing the way I created synthetic line-integrated data.

Interpreting the figure that comes up when doing

dax = mesh.plot_geometry_matrix(indbf=50, indchan=15, cam=cam)

doesn't seem trivial. When running this command (with this choice of indbf and indchan), I get the following:

What exactly am I looking at here? Do you see anything wrong? Thanks for your help!

[Didou09]

Hi @fsciortino , can you also include a plot of all your lines of sight ? (just do cam.plot() and copy-paste the figure, si I get an idea of the global geometry)

Also, can you copy-paste the whole figure produced by plot_inversion() ?

[Didou09]

plot_geometry_matrix() shows that your grid probably too thin with respect to your geometrical coverage (i.e.: you have a few pixels which are blank in the middle of the mesh. This is bad for the resolution.

• solution 1: try a slightly coarser grid
• solution 2: try bsplines of degree 1 (and send me the exception message you get if any)

[fsciortino]

Here is the result of cam.plot(), zoomed into the lower divertor:

(FYI, we don't actually have all these views; this is just an idealized test)

[fsciortino]

... and the result of mesh.plot_inversion:

[Didou09]

Thanks, so apparently the inversion gives measurements that are a pretty good fit to the synthetic data (apart from a few isolated channels). If course you can guess that it is not satisfactory, this is an illustration of the infinity of solutions to inverse problems.

My guess to improve the results:

• slightly coarser grid + bsplines of degree 1 (send the error message if any)

Also, which operator did you choose ?

[fsciortino]

OK, changing the mesh resolution to 5cm now gives this (with bsplines order of 1, which now worked fine -- I may have done something silly previously):

Does the order of magnitude of the data that I give to mesh.add_data matter for convergence?

I'm using the D1N2 operator.

[Didou09]

How did you compute the synthetic signal that you use as input for inversion ?

Just to clarify: do you compute the synthetic signal using the SOLPS simulation of the whole cross-section and then try to run an inversion on a small fraction only of that cross-section ? If this is what you are trying to do, it cannot work fine, because the solution you're looking for cannot account for the whole signal (i.e.: for the part of the signal that comes from other parts of the cross-section)

[fsciortino]

Essentially, yes. I've been computing the synthetic signals from SOLPS on the whole cross-section. I'm trying to reconstruct the lower divertor signal under the assumption that the signal measured in the upper SOL is negligible... but I will try again dropping this assumption and only using LOS that only go through the lower divertor. I'll get back to you with those results!

[Didou09]

Yes, try that, it should give better results

[Didou09]

Be conservative for a first test: only keep the LOS that do not get anything from outside of the region you're meshing. Then, if the result is satisfactory, You'll try to gradually add carefully chosen LOS

[fsciortino]

With this set of views

I'm now getting this

[fsciortino]

I'll check again my line-integration routine for synthetic data -- I wouldn't expect it to be a problem, but since I don't seem to get good results, it seems worth another check.

[Didou09]

On the second figure, there are 3 LOS that show a clear error around channel indices ~17, ~55, ~101. explore the geometry of these LOS to understand why: are they edge channels ?

More generally, explore the channels with the largest errors to understand how to better match your grid and basis functions with your geometry, and to identify issues.

We can arrange a discussion on Monday if you want.