[mesa-users] OpenMP data allocation issues

[Kevin Moore]

Hi all,

I’ve got a technical programming question for any OpenMP experts on here.

I’m implementing an alternative semiconvective mixing scheme in MESA, using the other_mlt() hook and my own run_star_extras.f. I basically need to call the Newton solver for each cell to figure out gradT, D, etc., and do this through my own external module that has the init(), do_stuff(), shutdown() functionality similar to most MESA modules. I initialize and allocate my module's array pointers for the Newton solver inside the extras_controls() subroutine in run_star_extras.f, which according to calls to omp_in_parallel() is executed in a non-parallel region. I then call my do_stuff() method inside the subroutine that s% other_mlt points to, and according to omp_in_parallel() this is in a parallel region. I finally deallocate all the Newton solver arrays inside the extras_after_evolve() subroutine in run_star_extras.f

Things work fine with my new prescription in serial evolution (OMP_NUM_THREADS=1), but when I turn on multiple threads things eventually break with "Segmentation fault - invalid memory reference”. All of my array pointers in my external module are set as !$OMP THREADPRIVATE, which is supposed to give each thread its own copy of the data to work with (that fixed terrible problems of multiple threads writing to the same arrays). But, with two threads active things will work fine for about 20 time steps, and then seg fault. With three or more threads, it seg faults the first time I try to call the Newton solver. I can trace the error to the required arrays not being allocated via an associated() test, but they were allocated at the beginning of the program.

My guess as an OpenMP novice is that something may be happening when threads are destroyed and recreated that’s screwing things up, or maybe not all of my array pointers are on the stack and can’t be copied to the threads correctly. I checked the OpenMP 3.1 spec document (http://openmp.org/wp/openmp-specifications/), but couldn’t find anything in the data management section that helped me pinpoint my problem.

My main question for others that have worked on similar issues is, "how should memory be allocated in parallel regions?”. Should I allocate and deallocate arrays around each call to the Newton solver inside the parallel regions? (that seems like it would slow things down a bit) Is what I’m doing correct, but I need to invoke some other OpenMP directives somewhere? Any help or resources would be much appreciated.

-Kevin