Feature CUDA compatability

.github/workflows/self-hosted-ci.yml

@@ -39,6 +39,13 @@ jobs:
           compiler-modules: "PrgEnv-cray craype-accel-nvidia60 cdt-cuda/22.05 cudatoolkit/11.2.0_3.39-2.1__gf93aa1c"
           # OpenMP flags from Nichols Romero (Argonne)
           fcflags: "-hnoacc -homp -O0"
+        - config-name: cuda-kernels
+          # Fall back to OpenACC 
+          rte-kernels: extern 
+          compiler-modules: "PrgEnv-nvidia nvidia craype-accel-nvidia60 cdt-cuda/21.09 cudatoolkit/11.2.0_3.39-2.1__gf93aa1c !cray-libsci_acc"
+          fcflags: "-g -O3 -acc -gopt -Mallocatable=03 -Mpreprocess -Minfo"
+          experimental: true
+
     env:
       # Core variables:
       FC: ftn
@@ -53,7 +60,8 @@ jobs:
     #
     # Checks-out repository under $GITHUB_WORKSPACE
     #
-    - uses: actions/checkout@v3
+    - name: Check out Fortran code
+      uses: actions/checkout@v3
     #
     # Check out data
     #
@@ -63,6 +71,15 @@ jobs:
         repository: earth-system-radiation/rrtmgp-data
         path: rrtmgp-data
     #
+    # Check out CUDA kernels if needed 
+    #
+    - name: Check out CUDA kernels 
+      if: matrix.config-name == 'cuda-kernels'
+      uses: actions/checkout@v3
+      with:
+        repository: earth-system-radiation/rte-rrtmgp-cuda-kernels
+        path: cuda-kernels
+    #
     # Finalize build environment
     #
     - name: Finalize build environment
@@ -90,14 +107,22 @@ jobs:
         # SLURM jobs, user home directories and HDF5 file locking are
         # incompatible on Daint:
         echo 'HDF5_USE_FILE_LOCKING=FALSE' >> "${GITHUB_ENV}"
-    #
     # Build libraries, examples and tests
     #
     - name: Build libraries
       run: |
         $FC --version
         make -j8 libs
     #
+    # Build library of CUDA kernels; copy to build directory and overwrite defaults
+    #
+    - name: Build CUDA kernels
+      if: matrix.config-name == 'cuda-kernels'
+      run: |
+        make -C cuda-kernels/build 
+        cp cuda-kernels/build/librtecudakernels.a build/librtekernels.a
+        cp cuda-kernels/build/librrtmgpcudakernels.a build/librrtmgpkernels.a
+    #
     # Run examples and tests (expect success)
     #
     - name: Build and run examples and tests (expect success)

examples/all-sky/Makefile

@@ -7,7 +7,7 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build
 # RRTMGP library, module files
 #
 LDFLAGS   += -L$(RRTMGP_BUILD)
-LIBS      += -lrrtmgp -lrte
+LIBS      += -lrrtmgpf -lrtef -lrrtmgpkernels -lrtekernels 
 FCINCLUDE += -I$(RRTMGP_BUILD)
 
 # netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. 

examples/rfmip-clear-sky/Makefile

@@ -7,7 +7,7 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build
 # RRTMGP library, module files
 #
 LDFLAGS   += -L$(RRTMGP_BUILD)
-LIBS      += -lrrtmgp -lrte
+LIBS      += -lrrtmgpf -lrtef -lrrtmgpkernels -lrtekernels 
 FCINCLUDE += -I$(RRTMGP_BUILD)
 
 # netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. 

rrtmgp-frontend/mo_gas_optics_rrtmgp.F90

@@ -260,8 +260,10 @@ function gas_optics_int(this,                             &
     !
     ! Gas optics
     !
-    !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
-    !$omp target enter data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
+    !$acc         enter data copyin(play, plev, tlay) 
+    !$omp target  enter data map(to:play, plev, tlay)
+    !$acc         enter data create(   jtemp, jpress, tropo, fmajor, jeta)
+    !$omp target  enter data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
     error_msg = compute_gas_taus(this,                       &
                                  ncol, nlay, ngpt, nband,    &
                                  play, plev, tlay, gas_desc, &
@@ -274,7 +276,7 @@ function gas_optics_int(this,                             &
     ! External source -- check arrays sizes and values
     ! input data sizes and values
     !
-    !$acc enter data copyin(tsfc, tlev) ! Should be fine even if tlev is not supplied
+    !$acc enter data        copyin(tsfc, tlev) ! Should be fine even if tlev is not supplied
     !$omp target enter data map(to:tsfc, tlev)
 
     if(check_extents) then
@@ -328,6 +330,8 @@ function gas_optics_int(this,                             &
     !$omp target exit data map(release:tsfc)
     !$acc        exit data      delete(jtemp, jpress, tropo, fmajor, jeta)
     !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
+    !$acc        exit data delete(play, plev, tlay)
+    !$omp target exit data map(release:play, plev, tlay)
   end function gas_optics_int
   !------------------------------------------------------------------------------------------
   !
@@ -373,23 +377,27 @@ function gas_optics_ext(this,                         &
     !
     ! Gas optics
     !
-    !$acc enter data create(jtemp, jpress, tropo, fmajor, jeta)
-    !$omp target enter data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
+    !$acc  data       copyin(play, plev, tlay) 
+    !$omp target data map(to:play, plev, tlay)
+    !$acc  data       create(   jtemp, jpress, tropo, fmajor, jeta)
+    !$omp target data map(alloc:jtemp, jpress, tropo, fmajor, jeta)
     error_msg = compute_gas_taus(this,                       &
                                  ncol, nlay, ngpt, nband,    &
                                  play, plev, tlay, gas_desc, &
                                  optical_props,              &
                                  jtemp, jpress, jeta, tropo, fmajor, &
                                  col_dry)
-    !$acc exit data delete(jtemp, jpress, tropo, fmajor, jeta)
-    !$omp target exit data map(release:jtemp, jpress, tropo, fmajor, jeta)
+    !$acc end         data
+    !$omp end target  data 
+    !$acc end         data
+    !$omp end target  data 
     if(error_msg  /= '') return
 
     ! ----------------------------------------------------------
     !
     ! External source function is constant
     !
-    !$acc enter data create(toa_src)
+    !$acc        enter data create(   toa_src)
     !$omp target enter data map(alloc:toa_src)
     if(check_extents) then
       if(.not. extents_are(toa_src, ncol, ngpt)) &
@@ -404,7 +412,7 @@ function gas_optics_ext(this,                         &
           toa_src(icol,igpt) = this%solar_source(igpt)
        end do
     end do
-    !$acc exit data copyout(toa_src)
+    !$acc        exit data copyout( toa_src)
     !$omp target exit data map(from:toa_src)
   end function gas_optics_ext
   !------------------------------------------------------------------------------------------
@@ -710,6 +718,7 @@ function compute_gas_taus(this,                       &
                           ncol, nlay, ngpt, optical_props%p)
         end select
       end if
+      ! Interpolation coefficients copyout 
       !$acc end        data
       !$omp end target data
       if (present(col_dry)) then

rrtmgp-kernels/api/mo_gas_optics_rrtmgp_kernels.F90

@@ -1,3 +1,14 @@
+! This code is part of
+! RRTM for GCM Applications - Parallel (RRTMGP)
+!
+! Contact: Eli Mlawer and Robert Pincus
+! email:  rrtmgp@aer.com
+!
+! Copyright 2015-,  Atmospheric and Environmental Research and
+! Regents of the University of Colorado.  All right reserved.
+!
+! Use and duplication is permitted under the terms of the
+!    BSD 3-clause license, see http://opensource.org/licenses/BSD-3-Clause
 module mo_gas_optics_rrtmgp_kernels
   use mo_rte_kind,      only : wp, wl
   use mo_rte_util_array,only : zero_array

rte-frontend/mo_optical_props.F90

@@ -275,8 +275,7 @@ function init_base(this, band_lims_wvn, band_lims_gpt, name) result(err_message)
     !
     ! Assignment
     !
-    if(allocated(this%band2gpt     )) deallocate(this%band2gpt)
-    if(allocated(this%band_lims_wvn)) deallocate(this%band_lims_wvn)
+    call this%finalize_base()
     allocate(this%band2gpt     (2,size(band_lims_wvn,2)), &
              this%band_lims_wvn(2,size(band_lims_wvn,2)))
     this%band2gpt      = band_lims_gpt_lcl
@@ -287,11 +286,13 @@ function init_base(this, band_lims_wvn, band_lims_gpt, name) result(err_message)
     ! Make a map between g-points and bands
     !   Efficient only when g-point indexes start at 1 and are contiguous.
     !
-    if(allocated(this%gpt2band)) deallocate(this%gpt2band)
     allocate(this%gpt2band(maxval(band_lims_gpt_lcl)))
     do iband=1,size(band_lims_gpt_lcl,dim=2)
       this%gpt2band(band_lims_gpt_lcl(1,iband):band_lims_gpt_lcl(2,iband)) = iband
     end do
+    !$acc        enter data copyin(this%band2gpt, this%gpt2band, this%band_lims_wvn)
+    !$omp target enter data map(to:this%band2gpt, this%gpt2band, this%band_lims_wvn)
+
   end function init_base
   !-------------------------------------------------------------------------------------------------
   function init_base_from_copy(this, spectral_desc) result(err_message)
@@ -326,6 +327,8 @@ end function is_initialized_base
   subroutine finalize_base(this)
     class(ty_optical_props),    intent(inout) :: this
 
+    !$acc        exit data     delete( this%band2gpt, this%gpt2band, this%band_lims_wvn)
+    !$omp target exit data map(release:this%band2gpt, this%gpt2band, this%band_lims_wvn)
     if(allocated(this%band2gpt)) deallocate(this%band2gpt)
     if(allocated(this%gpt2band)) deallocate(this%gpt2band)
     if(allocated(this%band_lims_wvn)) &
@@ -1094,12 +1097,15 @@ end function get_ngpt
   !> The first and last g-point of all bands at once
   !> dimension (2, nbands)
   !>
-  pure function get_band_lims_gpoint(this)
-    class(ty_optical_props), intent(in) :: this
-    integer, dimension(size(this%band2gpt,dim=1), size(this%band2gpt,dim=2)) &
-                                        :: get_band_lims_gpoint
+  function get_band_lims_gpoint(this)
+    class(ty_optical_props), target, intent(in) :: this
+    integer, dimension(:,:), pointer :: get_band_lims_gpoint
 
-    get_band_lims_gpoint = this%band2gpt
+    if(this%is_initialized()) then
+      get_band_lims_gpoint => this%band2gpt
+    else
+      get_band_lims_gpoint => NULL()
+    end if 
   end function get_band_lims_gpoint
   !>--------------------------------------------------------------------------------------------------------------------
   !>
@@ -1121,15 +1127,14 @@ end function convert_band2gpt
   !> Lower and upper wavenumber of all bands
   !> (upper and lower wavenumber by band) = band_lims_wvn(2,band)
   !>
-  pure function get_band_lims_wavenumber(this)
-    class(ty_optical_props), intent(in) :: this
-    real(wp), dimension(size(this%band_lims_wvn,1), size(this%band_lims_wvn,2)) &
-                                        :: get_band_lims_wavenumber
+  function get_band_lims_wavenumber(this)
+    class(ty_optical_props), target, intent(in) :: this
+    real(wp), dimension(:,:), pointer           :: get_band_lims_wavenumber
 
     if(this%is_initialized()) then
-      get_band_lims_wavenumber(:,:) = this%band_lims_wvn(:,:)
+      get_band_lims_wavenumber => this%band_lims_wvn
     else
-      get_band_lims_wavenumber(:,:) = 0._wp
+      get_band_lims_wavenumber => NULL()
     end if
   end function get_band_lims_wavenumber
   !>--------------------------------------------------------------------------------------------------------------------
@@ -1151,15 +1156,14 @@ end function get_band_lims_wavelength
   !> Bands for all the g-points at once
   !> dimension (ngpt)
   !>
-  pure function get_gpoint_bands(this)
-    class(ty_optical_props), intent(in) :: this
-    integer, dimension(size(this%gpt2band,dim=1)) &
-                                        :: get_gpoint_bands
+  function get_gpoint_bands(this)
+    class(ty_optical_props), target, intent(in) :: this
+    integer, dimension(:), pointer :: get_gpoint_bands
 
     if(this%is_initialized()) then
-      get_gpoint_bands(:) = this%gpt2band(:)
+      get_gpoint_bands => this%gpt2band
     else
-      get_gpoint_bands(:) = 0
+      get_gpoint_bands => NULL()
     end if
   end function get_gpoint_bands
   !>--------------------------------------------------------------------------------------------------------------------
@@ -1196,7 +1200,7 @@ end function expand
   !>
   !> Are the bands of two objects the same? (same number, same wavelength limits)
   !>
-  pure function bands_are_equal(this, that)
+  function bands_are_equal(this, that)
     class(ty_optical_props), intent(in) :: this, that
     logical                             :: bands_are_equal
 
@@ -1212,7 +1216,7 @@ end function bands_are_equal
   !> Is the g-point structure of two objects the same?
   !>   (same bands, same number of g-points, same mapping between bands and g-points)
   !>
-  pure function gpoints_are_equal(this, that)
+  function gpoints_are_equal(this, that)
     class(ty_optical_props), intent(in) :: this, that
     logical                             :: gpoints_are_equal
 

rte-frontend/mo_rte_lw.F90

@@ -142,6 +142,7 @@ function rte_lw(optical_props, top_at_1, &
                              0.2009319137_wp, 0.2292411064_wp, 0.0698269799_wp, 0._wp, &
                              0.1355069134_wp, 0.2034645680_wp, 0.1298475476_wp, 0.0311809710_wp], &
                              [max_gauss_pts, max_gauss_pts])
+    !$acc declare create(gauss_Ds, gauss_wts)
     ! ------------------------------------------------------------------------------------
     ncol  = optical_props%get_ncol()
     nlay  = optical_props%get_nlay()
@@ -310,7 +311,7 @@ function rte_lw(optical_props, top_at_1, &
     ! Compute the radiative transfer...
     !
     !$acc        data create(   sfc_emis_gpt, flux_up_loc, flux_dn_loc, gpt_flux_up, gpt_flux_dn)
-    !$omp target data map(alloc:sfc_emis_gpt, flux_up_loc, flux_dn_loc, gpt_flux_up, gpt_flux_dn)
+    !$omp target data map(alloc:sfc_emis_gpt, flux_up_loc, flux_dn_loc, gpt_flux_up, gpt_flux_dn) map(to:gauss_wts, gauss_Ds)
     call expand_and_transpose(optical_props, sfc_emis, sfc_emis_gpt)
     if(check_values) error_msg =  optical_props%validate()
     if(len_trim(error_msg) == 0) then ! Can't do an early return within OpenACC/MP data regions

tests/Makefile

@@ -7,7 +7,7 @@ RRTMGP_BUILD = $(RRTMGP_ROOT)/build
 # RRTMGP library, module files
 #
 LDFLAGS   += -L$(RRTMGP_BUILD)
-LIBS      += -lrrtmgp -lrte
+LIBS      += -lrrtmgpf -lrtef -lrrtmgpkernels -lrtekernels 
 FCINCLUDE += -I$(RRTMGP_BUILD)
 
 # netcdf Fortran module files has to be in the search path or added via environment variable FCINCLUDE e.g. 

tests/check_variants.F90

@@ -154,6 +154,7 @@ program rte_clear_sky_regression
   !
   call stop_on_err(gas_conc_array(1)%get_subset(1, ncol, gas_concs))
   call gas_conc_array(1)%reset()
+  print *, "Reset gas concs" ! Without this line there's an OpenACC error
   deallocate(gas_conc_array)
   ! ----------------------------------------------------------------------------
   ! load data into classes

tests/intel-codecov.sh

@@ -12,7 +12,6 @@ export PROF_DIR=$PWD
 # Environment variables for netCDF Fortran and C installations
 #
 export NFHOME=${HOME}/Applications/${FC}
-export NCHOME=/opt/local
 
 #
 # An Anaconda environent with modules needed for other python scripts
@@ -32,7 +31,7 @@ make -C ${RRTMGP_BUILD} -j 4 || exit 1
 #
 cd ${RRTMGP_ROOT}/examples/rfmip-clear-sky || exit 1
 export FCFLAGS+=" -I${RRTMGP_BUILD} -I${NFHOME}/include"
-export LDFLAGS+=" -L${RRTMGP_BUILD} -L${NFHOME}/lib -L${NCHOME}/lib -lrte -lrrtmgp -lnetcdff -lnetcdf"
+export LDFLAGS+=" -L${RRTMGP_BUILD} -L${NFHOME}/lib -lrte -lrrtmgp -lnetcdff"
 make clean || exit 1
 make -j 4  || exit 1
 python ./stage_files.py