Make LSODE integrator thread-safe in OpenMP parallel loops by adding THREADPRIVATE statements

[yantosca]

Overview

This is the companion PR to #118. We have now made the LSODE integrator thread-safe by adding $OMP THREADPRIVATE statements to the following variables:

1. RWORK in routine KppLsode (in int/lsode.f90);

   • This was also declared with SAVE to make it a global variable

2. IWORK in routine KppLsode (in int/lsode.f90)

   • This was also declared with SAVE to make it a global variable

3. All instances of common block /DLS001/ throughout int/lsode.f90

Validation

I used our fork of the KPP-Standalone, and modified the kpp_standalone.F90 driver program to include an OpenMP loop:

    !$OMP PARALLEL DO             &
    !$OMP DEFAULT( SHARED       ) &
    !$OMP PRIVATE( N, TIN, TOUT )
    DO N = 1, 10

       !$OMP CRITICAL
       print*, '### instance ', N
       !$OMP END CRITICAL

       ! Set ENV
       TIN          = T
       TOUT         = T + OperatorTimestep

       ! Set initial concentrations (C) and reacton rates (RCONST)
       ! to values read from the input file
       C            = Cinit
       RCONST       = R

       ! Integrate the mechanism for an operator timestep
       CALL Integrate( TIN, TOUT, ICNTRL, RCNTRL, ISTATUS, RSTATE, IERR )

    ENDDO
    !$OMP END PARALLEL DO

This will perform an integration on the same initial conditions 10 times. The C array is overwritten by each integration. This is sufficient to test if there is a parallelization error.

When compiled with OpenMP parallelization turned off, the integrations are done in order from 1..10 and we get this screen output:

$ ./kpp_standalone.exe Beijing_L1_20190701_0040.txt ref.log
 Processing sample: Beijing_L1_20190701_0040.txt
 Output file: ref.log
 ### instance            1
 ### instance            2
 ### instance            3
 ### instance            4
 ### instance            5
 ### instance            6
 ### instance            7
 ### instance            8
 ### instance            9
 ### instance           10
 Number of internal timesteps (from 3D run):    12
 Number of internal timesteps ( standalone):    92
 Hexit (from 3D run):     497.80
 Hexit ( standalone):      80.40
Warning: Number of internal steps do not match 3D grid cell
Warning: final timestep does not match 3D grid cell within 0.1%

final concentrations are written to the ref.log file. The difference in Hexit is because the initial conditions came from a run using the Rosenbrock integrator.

Then I recompiled the executable with -fopenmp , and we get this screen output:

$ ./kpp_standalone.exe Beijing_L1_20190701_0040.txt devmp.log
 Processing sample: Beijing_L1_20190701_0040.txt
 Output file: devmp.log
 ### instance           10
 ### instance            5
 ### instance            6
 ### instance            2
 ### instance            4
 ### instance            3
 ### instance            9
 ### instance            7
 ### instance            1
 ### instance            8
 Number of internal timesteps (from 3D run):    12
 Number of internal timesteps ( standalone):    92
 Hexit (from 3D run):     497.80
 Hexit ( standalone):      80.40
Warning: Number of internal steps do not match 3D grid cell
Warning: final timestep does not match 3D grid cell within 0.1%

Here you can see the instances are being done in non-sequential order, as the parallel loop is active. We get identical Hexit as the run with parallelization turned off. Also, the output in devmp.log was identical to ref.log as determined by a diff command.

[RolfSander]

Thanks for asking me to review this. However, I don't have much
experience with OpenMP and "thread-safe" coding, so I will not be able
to contribute much here.

[yantosca]

I spoke too soon. Compiling with GNU 12.2.0 resulted in this error:

[ 68%] Building Fortran object src/GEOS-Chem/KPP/fullchem/CMakeFiles/KPP.dir/gckpp_Integrator.F90.o
/n/netscratch/jacob_lab/Lab/ryantosca/tests/kpptests/lsode/test_lsode/CodeDir/src/GEOS-Chem/KPP/fullchem/gckpp_Integrator.F90:254:27:

  254 |       INTEGER :: IWORK(LIW), IPAR(1), ITOL, ITASK,         &
      |                           ^
Error: ‘iwork’ causes a section type conflict with ‘rwork’
/n/netscratch/jacob_lab/Lab/ryantosca/tests/kpptests/lsode/test_lsode/CodeDir/src/GEOS-Chem/KPP/fullchem/gckpp_Integrator.F90:253:33:

  253 |       REAL(kind=dp) :: RWORK(LRW), RPAR(1)
      |                                 ^
note: ‘rwork’ was declared here

I will look into this further.

[jimmielin]

Thanks @yantosca for these changes (and the helpful comments)! I am not sure about the section type conflict error or what could be causing it, so unfortunately I cannot help on that.

[yantosca]

@jimmielin: I wasn't able to make much headway with the section-type error. I decided to replace the COMMON /DLS001/ block with a derived-type object that is held !$OMP THREADPRIVATE. Still shaking down some more issues but at least this approach doesn't incur the section-type error. Will request a re-review once I get it working.