Thread initialization preferences

[pshriwise]

Context: I'm looking for some advice on how to handle initialization for threads with OpenMP in libMesh. I'm currently working on a library, XDG, that provides a DAGMC-like boundary to boundary ray tracing capability and tracklength tallies for libMesh meshes (and MFEM eventually). It is thread-safe and the downstream application (OpenMC) will access it from many threads using OpenMP threaded paralellism.

Trouble: OpenMC sets the number of OpenMP threads when it starts up, but inside the libMeshInit constructor the number of OpenMP threads will be reset to the number of threads passed to the constructor (default of1). libMeshInit is called inside the intermediate library XDG and as a result it now needs to communicate the number of threads to libMesh. I can currently do this by either a) exposing a threads argument in XDG's config/constructor and passing it through or b) collecting the number of threads set for OpenMP with omc_get_num_threads and passing it on to the libMeshInit constructor.

a) is overhead that isn't needed used XDG and b) requires me to add OpenMP as a dependency solely to query the number of threads set. Right now I'm doing b) because, despite the extra dependency, it's easier to reverse.

For clarity, XDG operations are the same regardless of the number of threads, and in this context libMesh doesn't either. This is entirely to ensure that number of threads isn't changed for OpenMC. And for some history, this previously wasn't an issue because OpenMC was initializing libMesh directly and could pass the number of threads OpenMC as needed.

An ideal situation for me would be to set the OMP_NUM_THREADS variable and have that be used, but I'm open to any suggestions for alternative approaches!

[roystgnr]

We have a --n-threads command line option that overrides the constructor argument. I'd call that more of a "workaround" than a "solution", though.

IMHO (b) is a perfectly fine solution.

I would say it also wouldn't be crazy to give libMesh a set_n_threads(n) setter to go along with the n_threads() getter; the setter would recreate the task_scheduler, call omp_set_num_threads if we're using OpenMP, etc. The biggest problem with that would be that, if we rely on libMesh (via TIMPI) to do the MPI_Init call, then we rely on our initial n_threads value to determine whether we should auto-request threading support from MPI, and I'm not sure we have any way to tell MPI we changed our minds later. The best a setter could do if we needed increased threading support and couldn't get it would be to check for the inconsistency and then scream and die with a nice error immediately rather than risk dying with a complicated error much later.

Checking $OMP_NUM_THREADS would have been a fine way to do things if we'd been doing that from the get-go, but if we suddently started checking it at this late date I bet there would be a lot of users who just happen to have it set to n_procs for their environment and who suddenly started trying to spawn n_procs*n_procs threads in what they'd intended to be an n_procs*1 threaded MPI-parallelism-only run.