Feature 1926 gridstat openmp

met/configure.ac

@@ -8,6 +8,13 @@ AC_CONFIG_HEADERS([config.h])
 
 AM_INIT_AUTOMAKE([1.9 foreign])
 
+# OpenMP
+
+AC_OPENMP()
+
+CPPFLAGS="${CPPFLAGS} ${OPENMP_CFLAGS}"
+LDFLAGS="${LDFLAGS} ${OPENMP_CFLAGS}"
+
 #
 # Look for the NetCDF library
 #

met/docs/Users_Guide/installation.rst

@@ -290,6 +290,14 @@ Enable compilation of the MODE-Graphics tool. Requires $MET_CAIRO and $MET_FREET
 
 Disable use of BLOCK4 in the compilation. Use this if you have trouble using PrepBUFR files.
 
+.. code-block:: none
+
+  --disable-openmp
+
+Disable compilation of OpenMP directives within the code which allows some code
+regions to benefit from thread-parallel execution. Runtime environment variable
+:code:`OMP_NUM_THREADS` controls the number of threads.
+
 Run the configure script with the **-help** argument to see the full list of configuration options.
 
 Make Targets

met/src/basic/vx_util/Makefile.am

@@ -69,6 +69,7 @@ libvx_util_a_SOURCES = ascii_table.cc ascii_table.h \
                GridOffset.h GridOffset.cc \
                observation.h observation.cc \
                stat_column_defs.h \
+	       handle_openmp.h handle_openmp.cc \
                RectangularTemplate.h RectangularTemplate.cc $(OPT_PYTHON_SOURCES)
 
 libvx_util_a_CPPFLAGS = ${MET_CPPFLAGS}

met/src/basic/vx_util/data_plane_util.cc

@@ -16,6 +16,10 @@ using namespace std;
 #include <string.h>
 #include <unistd.h>
 
+#ifdef _OPENMP
+  #include "omp.h"
+#endif
+
 #include "data_plane_util.h"
 #include "interp_util.h"
 #include "two_to_one.h"
@@ -247,83 +251,95 @@ void fractional_coverage(const DataPlane &dp, DataPlane &frac_dp,
       }
    }
 
-   // Build the grid template
-   GridTemplateFactory gtf;
-   GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon);
-
-   mlog << Debug(3)
-        << "Computing fractional coverage field using the "
-        << t.get_str() << " threshold and the "
-        << interpmthd_to_string(InterpMthd_Nbrhd) << "(" << gt->size()
-        << ") " << gt->getClassName() << " interpolation method.\n";
-
-   // Initialize the fractional coverage field
-   frac_dp = dp;
-   frac_dp.set_constant(bad_data_double);
-
-   // Compute the fractional coverage meeting the threshold criteria
-   for(x=0; x<dp.nx(); x++) {
-      for(y=0; y<dp.ny(); y++) {
-
-         // For a new column, reset the grid template and counts.
-         if(y == 0) {
-
-            // Initialize counts
-            n_vld = n_thr = 0;
-
-            // Sum all the points
-            for(gp  = gt->getFirstInGrid(x, y, dp.nx(), dp.ny());
-                gp != NULL;
-                gp  = gt->getNextInGrid()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld++;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
-            }
-         }
-         // Subtract off the bottom edge, shift up, and add the top.
-         else {
-
-            // Subtract points from the the bottom edge
-            for(gp  = gt->getFirstInBotEdge();
-                gp != NULL;
-                gp  = gt->getNextInBotEdge()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld--;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr--;
-            }
-
-            // Increment Y
-            gt->incBaseY(1);
-
-            // Add points from the the top edge
-            for(gp  = gt->getFirstInTopEdge();
-                gp != NULL;
-                gp  = gt->getNextInTopEdge()) {
-               if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
-               n_vld++;
-               if(t.check(v,
-                  (use_climo ? cmn->get(gp->x, gp->y) : bad),
-                  (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
-            }
-         }
-
-         // Check for enough valid data and compute fractional coverage
-         if((double)(n_vld)/gt->size() >= vld_t && n_vld != 0) {
-            frac_dp.set((double) n_thr/n_vld, x, y);
-         }
-
-      } // end for y
-
-      // Increment X
-      if(x < (dp.nx() - 1)) gt->incBaseX(1);
-
-   } // end for x
-
-   delete gt;
+#pragma omp parallel default(none)                 \
+   shared(mlog, dp, frac_dp, width, wrap_lon, t)   \
+   shared(use_climo, cmn, csd, vld_t, bad)         \
+   private(x, y, n_vld, n_thr, gp, v)
+   {
+
+     // Build the grid template
+     GridTemplateFactory gtf;
+     GridTemplate* gt = gtf.buildGT(shape, width, wrap_lon);
+
+#pragma omp single
+     {
+       mlog << Debug(3)
+            << "Computing fractional coverage field using the "
+            << t.get_str() << " threshold and the "
+            << interpmthd_to_string(InterpMthd_Nbrhd) << "(" << gt->size()
+            << ") " << gt->getClassName() << " interpolation method.\n";
+
+       // Initialize the fractional coverage field
+       frac_dp = dp;
+       frac_dp.set_constant(bad_data_double);
+     }
+
+     // Compute the fractional coverage meeting the threshold criteria
+#pragma omp for schedule (static)
+     for(x=0; x<dp.nx(); x++) {
+        for(y=0; y<dp.ny(); y++) {
+
+           // For a new column, reset the grid template and counts.
+           if(y == 0) {
+
+              // Initialize counts
+              n_vld = n_thr = 0;
+
+              // Sum all the points
+              for(gp  = gt->getFirstInGrid(x, y, dp.nx(), dp.ny());
+                  gp != NULL;
+                  gp  = gt->getNextInGrid()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld++;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
+              }
+           }
+           // Subtract off the bottom edge, shift up, and add the top.
+           else {
+
+              // Subtract points from the the bottom edge
+              for(gp  = gt->getFirstInBotEdge();
+                  gp != NULL;
+                  gp  = gt->getNextInBotEdge()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld--;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr--;
+              }
+
+              // Increment Y
+              gt->incBaseY(1);
+
+              // Add points from the the top edge
+              for(gp  = gt->getFirstInTopEdge();
+                  gp != NULL;
+                  gp  = gt->getNextInTopEdge()) {
+                 if(is_bad_data(v = dp.get(gp->x, gp->y))) continue;
+                 n_vld++;
+                 if(t.check(v,
+                    (use_climo ? cmn->get(gp->x, gp->y) : bad),
+                    (use_climo ? csd->get(gp->x, gp->y) : bad))) n_thr++;
+              }
+           }
+
+           // Check for enough valid data and compute fractional coverage
+           if((double)(n_vld)/gt->size() >= vld_t && n_vld != 0) {
+              frac_dp.set((double) n_thr/n_vld, x, y);
+           }
+
+        } // end for y
+
+        // Increment X
+        if(x < (dp.nx() - 1)) gt->incBaseX(1);
+
+     } // end for x
+
+     delete gt;
+
+   } // End of omp parallel
 
    return;
 }

met/src/basic/vx_util/handle_openmp.cc

@@ -0,0 +1,51 @@
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+// ** Copyright UCAR (c) 1992 - 2021
+// ** University Corporation for Atmospheric Research (UCAR)
+// ** National Center for Atmospheric Research (NCAR)
+// ** Research Applications Lab (RAL)
+// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+
+///////////////////////////////////////////////////////////////////////////////
+
+#include <cstdlib>
+
+#ifdef _OPENMP
+  #include "omp.h"
+#endif
+
+#include "vx_log.h"
+#include "handle_openmp.h"
+
+///////////////////////////////////////////////////////////////////////////////
+
+void init_openmp() {
+
+#ifdef _OPENMP
+
+  // If OMP_NUM_THREADS was not set, use the OpenMP API to set the thread count
+  // to 1 thread only.
+  const char* env_omp_num_threads = std::getenv("OMP_NUM_THREADS");
+  if (!env_omp_num_threads) {
+    mlog << Debug(2) << "OMP_NUM_THREADS is not set." 
+     << " Defaulting to 1 thread. \n";
+    omp_set_num_threads(1);
+  }
+
+#pragma omp parallel
+#pragma omp single
+  {
+     mlog << Debug(2) << "OpenMP running on " 
+       << omp_get_num_threads() << " thread(s).\n";
+  }
+
+#else  /* _OPENMP */
+
+  mlog << Debug(2) << "OpenMP disabled.\n";
+
+#endif  /* _OPENMP */
+
+}
+
+///////////////////////////////////////////////////////////////////////////////
+

met/src/basic/vx_util/handle_openmp.h

@@ -0,0 +1,24 @@
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+// ** Copyright UCAR (c) 1992 - 2021
+// ** University Corporation for Atmospheric Research (UCAR)
+// ** National Center for Atmospheric Research (NCAR)
+// ** Research Applications Lab (RAL)
+// ** P.O.Box 3000, Boulder, Colorado, 80307-3000, USA
+// *=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*
+
+///////////////////////////////////////////////////////////////////////////////
+
+#ifndef __HANDLE_OPENMP_H___ 
+#define __HANDLE_OPENMP_H___ 
+
+///////////////////////////////////////////////////////////////////////////////
+
+void init_openmp();
+
+///////////////////////////////////////////////////////////////////////////////
+
+#endif   /*  __HANDLE_OPENMP_H__  */
+
+///////////////////////////////////////////////////////////////////////////////
+
+

met/src/tools/core/ensemble_stat/ensemble_stat.cc

@@ -92,6 +92,8 @@ using namespace std;
 #include "nc_obs_util.h"
 #include "nc_point_obs_in.h"
 
+#include "handle_openmp.h"
+
 ////////////////////////////////////////////////////////////////////////
 
 static void process_command_line  (int, char **);
@@ -171,6 +173,9 @@ static void set_compress(const StringArray &);
 
 int main(int argc, char *argv[]) {
 
+   // Set up OpenMP (if enabled)
+   init_openmp();
+
    // Set handler to be called for memory allocation error
    set_new_handler(oom);
 

met/src/tools/core/grid_stat/grid_stat.cc

@@ -127,6 +127,8 @@ using namespace std;
 #include <sys/types.h>
 #include <unistd.h>
 
+#include "handle_openmp.h"
+
 #include "grid_stat.h"
 
 #include "vx_statistics.h"
@@ -185,6 +187,9 @@ static bool read_data_plane(VarInfo* info, DataPlane& dp, Met2dDataFile* mtddf,
 
 int main(int argc, char *argv[]) {
 
+   // Set up OpenMP (if enabled)
+   init_openmp();
+
    // Set handler to be called for memory allocation error
    set_new_handler(oom);
 

met/src/tools/other/gen_ens_prod/gen_ens_prod.cc

@@ -45,6 +45,8 @@ using namespace std;
 #include "nc_obs_util.h"
 #include "nc_point_obs_in.h"
 
+#include "handle_openmp.h"
+
 ////////////////////////////////////////////////////////////////////////
 
 static void process_command_line(int, char **);
@@ -81,6 +83,9 @@ static void set_ctrl_file  (const StringArray &);
 
 int main(int argc, char *argv[]) {
 
+   // Set up OpenMP (if enabled)
+   init_openmp();
+
    // Set handler to be called for memory allocation error
    set_new_handler(oom);
 

test/config/GridStatConfig_rtma

@@ -111,7 +111,7 @@ interp = {
    type = [
       { method = NEAREST; width  = 1; },
       { method = UW_MEAN; width  = 3; },
-      { method = MIN;     width  = 1; },
+      { method = MIN;     width  = 3; },
       { method = MAX;     width  = 3; }
    ];
 }

test/xml/unit_grid_stat.xml

@@ -75,6 +75,32 @@
     </output>
   </test>
 
+  <test name="grid_stat_GRIB2_NAM_RTMA_NP2">
+    <exec>&MET_BIN;/grid_stat</exec>
+    <env>
+      <pair><name>OMP_NUM_THREADS</name>  <value>2</value></pair>
+      <pair><name>OUTPUT_PREFIX</name>  <value>GRIB2_NAM_RTMA_NP2</value></pair>
+    </env>
+    <param> \
+      &DATA_DIR_MODEL;/grib2/nam/nam_2012040900_F012_gRtma.grib2 \
+      &DATA_DIR_OBS;/rtma/rtma_2012051712_F000.grib2 \
+      &CONFIG_DIR;/GridStatConfig_rtma \
+      -outdir &OUTPUT_DIR;/grid_stat -v 1
+    </param>
+    <output>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V.stat</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_fho.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_ctc.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_cts.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_cnt.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_sl1l2.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_nbrctc.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_nbrcts.txt</stat>
+      <stat>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_nbrcnt.txt</stat>
+      <grid_nc>&OUTPUT_DIR;/grid_stat/grid_stat_GRIB2_NAM_RTMA_NP2_120000L_20120409_120000V_pairs.nc</grid_nc>
+    </output>
+  </test>
+
   <test name="grid_stat_GRIB1_NAM_STAGE4">
     <exec>&MET_BIN;/grid_stat</exec>
     <env>