Subtle OpenMP problem in 2D/3D spherical coordinates (fixed)

[dullemond]

Recently, Mario Flock and Alicia Trabold reported a strange behavior with OpenMP for spherical coordinates in 2D/3D, where the temperatures near the midplane became slightly lower if you use more than 1 threads in OpenMP. The effect was small, but noticable.

After some serious bug hunting, I found that the problem was due to the fact that the amrray_sint1 etc variables were declared threadprivate, whereas they shouldn't be. The result of this erroneous threadprivate was that for all threads but thread 0 the "self-illumination" of the model did not work. For spherical coordinates with a theta-grid that does not go all the way to the pole, it can (and should) happen that some photons leave the grid but later re-enter again. This does not occur in cartesian coordinates, and it also does not happen in 1D spherical coordinates, or in 2D/3D spherical coordinates where the theta grid spans from 0 to pi. But it can happen if the theta grid spans from, say, 0.3 to pi-0.3, or from 0.3 to pi/2. For simplicity let's call this "self-illumination".

For the serial model, this self-illumination effect works well (as far as I know). However, in OpenMP with setthreads 4 (say), this self-illumination fails for all threads except thread 0. The effect of self-illumination is weak, so the differences are not too big (a few percent at most). But it is still a bug that had to be fixed.

So I removed the threadprivate for the amrray_sint1 etc variables. I also removed the threadprivate for thetagrid_cross_equator. With these fixes, the problem seems to be solved. The results for setthreads 4 are now the same as for setthreads 1 (or serial, which is the same).

So please upgrade to a version of 20 August 2023 or later.

Thanks, Alicia and Mario, for reporting this.