diff --git a/src/basic/vx_util/vx_util.h b/src/basic/vx_util/vx_util.h index 84638028ff..c50e7d16fc 100644 --- a/src/basic/vx_util/vx_util.h +++ b/src/basic/vx_util/vx_util.h @@ -37,6 +37,7 @@ #include "fix_float.h" #include "get_filenames.h" #include "grib_constants.h" +#include "GridTemplate.h" #include "int_array.h" #include "interp_mthd.h" #include "interp_util.h" diff --git a/src/libcode/vx_shapedata/shapedata.cc b/src/libcode/vx_shapedata/shapedata.cc index 9051d00819..cbad2eb497 100644 --- a/src/libcode/vx_shapedata/shapedata.cc +++ b/src/libcode/vx_shapedata/shapedata.cc @@ -15,8 +15,9 @@ // // Mod# Date Name Description // ---- ---- ---- ----------- -// 000 11-05-31 Halley Gotway Adapated from wrfdata.cc. -// 001 14-05-29 Halley Gotway Add ShapeData::n_objects() +// 000 05/31/11 Halley Gotway Adapated from wrfdata.cc +// 001 05/29/14 Halley Gotway Add ShapeData::n_objects() +// 002 11/02/23 Halley Gotway MET #2724 add OpenMP to convolution // /////////////////////////////////////////////////////////////////////////////// @@ -513,160 +514,110 @@ double ShapeData::get_attr(const ConcatString &attr_name, /////////////////////////////////////////////////////////////////////////////// +void ShapeData::conv_filter_circ(int diameter, double vld_thresh) { + const char *method_name = "ShapeData::conv_filter_circ() -> "; + GridPoint *gp = nullptr; + int x, y, n_vld; + double v, v_sum; + DataPlane conv_dp; -void ShapeData::conv_filter_circ(int diameter, double vld_thresh) - -{ - -int x, y, xx, yy, u, v; -int dn, fn, nn; -int vpr, upr; -int count, bd_count; -double center, cur, sum; -double dx, dy, dist; -double vld_ratio; -const int nx = data.nx(); -const int ny = data.ny(); -bool * f = (bool *) 0; -bool center_bad = false; -DataPlane in_data = data; -const bool vld_thresh_one = is_eq(vld_thresh, 1.0); - - -if ( (diameter%2 == 0) || (diameter < 3) ) { - - mlog << Error << "\nShapeData::conv_filter_circ() -> " - << "diameter must be odd and >= 3 ... diameter = " - << diameter << "\n\n"; - - exit(1); - -} - -const int radius = (diameter - 1)/2; - -const vector * in = &(in_data.Data); - vector * out = &(data.Data); - -f = new bool [diameter*diameter]; - - // - // set up the filter - // - -for (y=0; y= 3 ... diameter = " + << diameter << "\n\n"; + exit(1); } -} - - // - // do the convolution - // - -dn = -1; - -for(y=0; y= ny) ) continue; - - vpr = v + radius; - - for(u=-radius; u<=radius; ++u) { - - xx = x + u; - - if ( (xx < 0) || (xx >= nx) ) continue; - - upr = u + radius; - - fn = STANDARD_XY_YO_N(diameter, upr, vpr); - - if ( !(f[fn]) ) continue; - - nn = STANDARD_XY_YO_N(nx, xx, yy) ; - - cur = (*in)[nn]; - - if( ::is_bad_data(cur) ) { bd_count++; continue; } - - sum += cur; - - count++; - - } // for v - - } // for u - - // - // If the center of the convolution contains bad data and the ratio - // of bad data in the convolution area is too high, set the convoled - // value to bad data. - // - - if ( count == 0 ) sum = bad_data_double; - else { +#pragma omp parallel default(none) \ + shared(mlog, data, conv_dp, diameter) \ + shared(vld_thresh, bad_data_double) \ + private(x, y, n_vld, v, v_sum, gp) + { + + // Build the grid template with shape circle and wrap_lon false + GridTemplateFactory gtf; + GridTemplate* gt = gtf.buildGT(GridTemplateFactory::GridTemplate_Circle, + diameter, false); + +#pragma omp single + { + // Initialize the convolved field to bad data + conv_dp = data; + conv_dp.set_constant(bad_data_double); + } - vld_ratio = ((double) count)/(bd_count + count); + // Compute the convolved values +#pragma omp for schedule (static) + for(x=0; xgetFirstInGrid(x, y, data.nx(), data.ny()); + gp != nullptr; + gp = gt->getNextInGrid()) { + v = data.get(gp->x, gp->y); + if(::is_bad_data(v)) continue; + n_vld += 1; + v_sum += v; + } + } + // Subtract off the bottom edge, shift up, and add the top + else { + + // Subtract points from the the bottom edge + for(gp = gt->getFirstInBotEdge(); + gp != nullptr; + gp = gt->getNextInBotEdge()) { + v = data.get(gp->x, gp->y); + if(::is_bad_data(v)) continue; + n_vld -= 1; + v_sum -= v; + } + + // Increment Y + gt->incBaseY(1); + + // Add points from the the top edge + for(gp = gt->getFirstInTopEdge(); + gp != nullptr; + gp = gt->getNextInTopEdge()) { + v = data.get(gp->x, gp->y); + if(::is_bad_data(v)) continue; + n_vld += 1; + v_sum += v; + } + } - if ( center_bad && (vld_ratio < vld_thresh) ) sum = bad_data_double; - else sum /= count; + // Check for enough valid data and compute the mean value + if((double)(n_vld)/gt->size() >= vld_thresh && n_vld != 0) { + conv_dp.set((double) v_sum/n_vld, x, y); + } - } + } // end for y - (*out)[dn] = sum; + // Increment X + if(x < (data.nx() - 1)) gt->incBaseX(1); - } // for y + } // end for x -} // for x + delete gt; - // - // done - // + } // End of omp parallel -if ( f ) { delete [] f; f = (bool *) 0; } - -return; + // Save the result + data = conv_dp; + return; } - //////////////////////////////////////////////////////////////////////// diff --git a/src/tools/core/mode/mode.cc b/src/tools/core/mode/mode.cc index 88e18d0b40..a9277a0884 100644 --- a/src/tools/core/mode/mode.cc +++ b/src/tools/core/mode/mode.cc @@ -56,6 +56,7 @@ // 019 04/01/19 Fillmore Add FCST and OBS units. // 020 07/06/22 Howard Soh METplus-Internal #19 Rename main to met_main // 021 06/09/23 Albo Major changes for multivariate mode +// 022 11/02/23 Halley Gotway MET #2724 add OpenMP to convolution // //////////////////////////////////////////////////////////////////////// @@ -73,6 +74,7 @@ using namespace std; #include #include "main.h" +#include "handle_openmp.h" #include "string_array.h" #include "mode_usage.h" #include "mode_frontend.h" @@ -131,6 +133,9 @@ int met_main(int argc, char * argv []) string s; const char * user_config_filename = 0; + // Set up OpenMP (if enabled) + init_openmp(); + for (j=0,n=0; j