Optimisation of the solid angle computation with aperture obstruction

[obrejank]

This PR is a follow up to #663 providing major performance improvement to all to the 3 functions compute_solid_angle_apertures_light* (which have also all been merged into one) thanks to:

• modification of the computation routines to make them nogil,
• custom algorithm for polygon intersection which (much faster than the python module Polygon and even 3-5x faster than direct calls to GPC, the C library it reslies on)
• OMP parallelisation of the loop on mesh points (the largest loop by far)
• better management of polygon orientation to avoid checking for each point of the plasma.

In my tests, using the Intel compiler, the speedup was about 60x in sequential (OMP_NUM_THREADS=1). The computation seems to be mostly memory-bound so the parallelisation efficiency is not very large (x3 with 8 threads), but it is still better than nothing and I'm not sure much more can be gained.

The function compute_solid_angle_apertures_visibility could probably be merged in the same way but it is not called by the pipeline tests and I did not want to make changes I couldn't validate. Merging it with the rest would allow removing the current dependence to the python module Polygon.

The CI pipeline fails because of the matplotlib._contour error, which is not related to my changes.

[Didou09]

Thank you very much @obrejank , I'll fix the unit tests (matplotlib._contour), and it should work soon !

[Didou09]

OK, the workflow is fixed,

But, when running the full workflow (a 3x3 matrix with python 3.7, 3.8, 3.9 on Ubuntu, Windows and MacOS), it appears that the installation of tofu does not work anymore on MacOS:

https://github.com/ToFuProject/tofu/actions/runs/3620729344/jobs/6103424271

The error message suggests an issue with the compilation of tofu with OpenMP.

I have checked the following :

• I have temporarily removed testing on MacOS to check that the unit tests are all passing on Ubuntu and Windows
• I have checked that the devel branch still works on MacOS => this issue is specific to this branch and is not due to any change on the MacOS platforms themselves.

Hence, this issue is specific to MacOS and specific to this branch.

Any chance you could take a closer look into this MacOS-specifc issue @obrejank ?
Is there anything that you have modified that could have changed the way tofu compiles on MacOS ?

To edit the file that runs the unit tests on the testing matrix, you can play with the following file:
.github/workflows/complete-matrix.yml

[obrejank]

I only had a quick look, but the final error messages says:
error: Command "gcc -Wno-unused-result -Wsign-compare -Wunreachable-code -DNDEBUG -g -fwrapv -O3 -Wall -I/usr/local/opt/sqlite/include -I/usr/local/opt/sqlite/include -I/Users/runner/hostedtoolcache/Python/3.7.15/x64/lib/python3.7/site-packages/numpy/core/include -I/Users/runner/hostedtoolcache/Python/3.7.15/x64/include/python3.7m -I/Users/runner/hostedtoolcache/Python/3.7.15/x64/include/python3.7m -c tofu/geom/_GG.c -o build/temp.macosx-10.15-x86_64-3.7/tofu/geom/_GG.o -O3 -Wall -fno-wrapv" failed with exit status 1
where the last few flags are the "extra flags" defined in setup.py, which should include the -fopenmp flag returned by is_openmp_installed. From what I understand, it can only return nothing if the try except finally clause failed, but I don't have a Mac to test it on and see why the subprocess.call failed. Either the command itself failed or it was the write above.

[obrejank]

If you look at the output in https://github.com/ToFuProject/tofu/actions/runs/3622194828/jobs/6106600138, you will see that MacOS uses clang when you call gcc:
subprocess.call(..) -> CompletedProcess(args=['gcc', '-fopenmp', 'test.c'], returncode=1, stdout=b'', stderr=b"clang: error: unsupported option '-fopenmp'\nclang: error: unsupported option '-fopenmp'\n")
which is absurd and stupid, but confirmed here. Since nothing was changed in recent versions of ToFu, I assume that this comes from an update in github's VMs.

[Didou09]

But the devel branch seems to run fin through all tests on MacOS... (I re-ran all them 1 hour ago), so this issue seems to be specific to this branch for some obscure reason...

[Didou09]

A possible fix could be to drop temporarily the support for MacOS, but several users are using MacOS.... It has to be fixed at some point

[Didou09]

So, I tried providing the correct argument:

But it fails because the default CLANG install on MacOS does not support OpenMP:

https://github.com/ToFuProject/tofu/actions/runs/3622194828/jobs/6106600138

subprocess.call(..) -> CompletedProcess(args=['gcc', '-fopenmp=libomp', 'test.c'], returncode=1, stdout=b'', stderr=b"clang: error: unsupported argument 'libomp' to option 'fopenmp='\nclang: error: unsupported argument

https://www.positioniseverything.net/clang-error-unsupported-option-fopenmp

And apparently we have to make sure it's llvm that's installed and not just the default clang

What I still don't understand is why it works on other branches...

[Didou09]

Identified source of compiling issue:

• In the prevous versions, openmp was an optional dependency, indeed the only module that was explicitly doing a cimport openmp was ./tofu/geom/_openmp_tools.pyx, but it was a conditional import

The rest of the module is coded with a condition on the existence of openmp:

• On this branch, _GG.pyx now has an unconditional cimport openmp

So basically, the compiling issue we have here was probably existing already in the previous version, but it was not apparent because for MacOS TOFU_OPENMP_ENABLED was False and MacOS tests were run without parallelization.

We have 2 options:

1. Implement the test on TOFU_OPENMP_ENABLED in _GG.pyx (=> lots of extra work)
2. Make sure openmp works on MacOS => better but difficult

[Didou09]

I was testing this a bit more extensively and found a bug that can be reproduced with the MWE script attached.

It creates a set of sensors, with 2 apertures each.
Each sensors is bigger than the previous.
The resulting solid angle map should increase accordingly (on devel).

On this branch, the solid angle map becomes degenerate when the sensor becomes big, and even decreases instead of increasing, for some reason.
It may look like a rounding error?
NCAM_tutorial_tofu.txt

import numpy as np
import matplotlib.pyplot as plt
plt.switch_backend('Qt5Agg')
plt.ion()

import tofu as tf


conf = tf.load_config('SPARC-V0')

# --------------------------------
# Instanciate an empty Collection

coll = tf.data.Collection()

# add mesh
coll.add_mesh_2d_rect(
    key='m0',
    res=0.05,
    domain=[(1.3, 2.4), (-1, 1)],
    crop_poly=conf,
)



# ----------------------
# use the built-in method

coll.add_camera_pinhole(
    #keys / identifiers
    key='cam0',
    key_pinhole='pinhole0',
    key_diag='ncam0',
    cam_type='1d',
    # coordinates of the pinhole (x, y, z) or (R, Z, phi)
    R=5,
    z=0,
    phi=0,
    # 3 angles to orientate the camera
    # theta is the angle vs horizontal plane (0 => pointing towards R=inf)  in a poloidal plane
    theta=np.pi,
    # dphi is an angle vs the poloidal plane
    dphi=0.,
    # tilt is an angular tilt of the camera around its own axis
    tilt=np.pi/2.,
    # length between pinhole and sensor plane
    focal=10,
    # nb / size of pixels
    pix_nb=5,
    pix_size=0.1,
    pix_spacing=1.,
    # pinhole: circular or rectangular
    pinhole_radius=0.1,
    pinhole_size=None,
    # tokamak (Config) to be associated to it
    config=conf,
    # shall LOS computation start right away?
    compute=True,
)

# --------------------------
# extract los pts from ncam0

# pts (2, ncam)
ptsx, ptsy, ptsz = coll.get_rays_pts('ncam0_cam0_los')

# units vectors of LOS (1, ncam)
vectx, vecty, vectz = coll.get_rays_vect('ncam0_cam0_los')

# prepare sensor size
size = 0.001

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

# prepare radius of collimator
radius = 0.02

# prepare distance between sensor and apertures
d0 = 0.2
d1 = 3.2

# ------------------
# Build the geometry

# get starting points of
doptics = {}
ncam = ptsx.shape[1]
size = size * (np.arange(ncam)*10 + 1)
for ii in range(ncam):
# for ii in range(1):

    # set keys of camera, aperture 0 and 1 (2 ends of the collimator)
    kcam = f'ncam1_c{ii}'
    kap0 = f'cam1_ap{ii}0'
    kap1 = f'cam1_ap{ii}1'

    # center of sensor
    cent = np.r_[ptsx[0, ii], ptsy[0, ii], ptsz[0, ii]]

    # los unit vector
    nin = np.r_[vectx[0, ii], vecty[0, ii], vectz[0, ii]]
    e0 = np.r_[-nin[1], nin[0], 0.]
    e0 = e0 / np.linalg.norm(e0)
    e1 = np.cross(nin, e0)

    # geometry dict for single-pixel camera 1d
    dgeom = {
        'cents_x': np.r_[cent[0]],
        'cents_y': np.r_[cent[1]],
        'cents_z': np.r_[cent[2]],
        'nin_x': np.r_[-1],
        'nin_y': np.r_[0],
        'nin_z': np.r_[0],
        'e0_x': np.r_[0],
        'e0_y': np.r_[-1],
        'e0_z': np.r_[0],
        'e1_x': np.r_[0],
        'e1_y': np.r_[0],
        'e1_z': np.r_[1],
        'outline_x0': size[ii] * np.r_[-1, 1, 1, -1],
        'outline_x1': size[ii] * np.r_[-1, -1, 1, 1],
    }

    # add single-pixel 1d camera
    coll.add_camera_1d(key=kcam, dgeom=dgeom)

    # add apertures
    dgeom = {
        'cent': cent + d0 * nin,
        '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'] = cent + d1 * nin
    dgeom['outline_x0'] = radius * np.cos(theta)
    dgeom['outline_x1'] = radius * np.sin(theta)
    coll.add_aperture(key=kap1, **dict(dgeom))

    # fill in doptics
    doptics[kcam] = [kap0, kap1]


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

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


# --------------------------
# Compute vos

dvos, dref = coll.compute_diagnostic_vos(
    # which diag and mesh to use
    key_diag='ncam1',
    key_cam=None,
    key_mesh='m0',
    # sampling resolution (to be tuned by multiple tests until you find reasonable convergence)
    res_RZ=0.03,
    res_phi=0.01,
    # extra flags
    visibility=False,
    verb=None,
    store=False,
)


# --------------------
# prepare plot

fig = plt.figure(figsize=(15, 6))
ax = fig.add_axes([0.03, 0.02, 0.95, 0.95], aspect='equal')
ax.set_xlabel('R (m)', size=12, fontweight='bold')
ax.set_ylabel('Z (m)', size=12, fontweight='bold')

dax = {'cross': ax}
for ii in range(ncam):

    kcam = f'ncam1_c{ii}'

    # cam0
    dax = coll.plot_diagnostic_vos(
        key='ncam1',
        key_cam=kcam,
        indch=None,
        elements='o',
        dvos=dvos,
        # specify colormaps
        proj='cross',
        dax=dax,
        plot_config=conf,
        plot_colorbar=False,
    )

    # add detector size
    ax.text(
        np.hypot(ptsx[0, ii], ptsy[0, ii]) + 0.05,
        ptsz[0, ii],
        f'{size[ii]}',
        c='k',
        horizontalalignment='left',
        verticalalignment='center',
    )

    # add max solid angle
    sang = dvos[kcam]['sang']['data']
    ax.text(
        np.hypot(ptsx[1, ii], ptsy[1, ii]) - 0.05,
        ptsz[1, ii],
        f"{np.max(sang):3.2e}",
        c='k',
        horizontalalignment='right',
        verticalalignment='center',
    )

On devel:

On this branch: