Merge pull request #244 from ngs333/cpp_dev

create_xgrid_2Dx2D_order2 C++ routine working for GPU.
NOAA-GFDL · Sep 6, 2023 · 84af6b4 · 84af6b4
2 parents f33420f + 25aaecf
commit 84af6b4
Show file tree

Hide file tree

Showing 42 changed files with 1,103 additions and 679 deletions.
diff --git a/.gitmodules b/.gitmodules
diff --git a/Makefile.am b/Makefile.am
@@ -52,10 +52,6 @@ SUBDIRS = postprocessing/combine_restarts \
           tools/simple_hydrog/rmvpr \
           tools/nc_null_check
 
-if ENABLE_OCEAN_MODEL_GRID_GENERATOR
-SUBDIRS += tools/ocean_model_grid_generator
-endif
-
 if WITH_CHECK_PROGS
 SUBDIRS += t
 endif
diff --git a/README.md b/README.md
@@ -8,7 +8,7 @@ These tools were largely written by members of the GFDL
 primarily for use in the
 [Flexible Modeling System](https://www.gfdl.noaa.gov/fms) (FMS)
 [Runtime Environment](https://www.gfdl.noaa.gov/fre) (FRE) supporting the
-work of the 
+work of the
 [Geophysical Fluid Dynamics Laboratory](https://www.gfdl.noaa.gov)
 (GFDL).
 
@@ -58,16 +58,8 @@ The tools available in FRE-NCtools are:
 
 
 ### Other Tools
-The [Ocean Model Grid Generator](https://github.com/NOAA-GFDL/ocean_model_grid_generator)
-is a collection of tools for creating finite element spherical tripolar grids for
-GFDL's MOM based ocean models. Unlike the other tools, NCTools includes it as a submodule,
-and also it is a Python3 project. Because of the former attribute, recursive
-cloning (see below) is recommended. Because of the latter attribute,
-the users python environment may need modification and/or options to the
-autotools configure command may need to be specified. This includes the
-``--disable-ocean-model-grid-generator`` option (default is enable)
-and the ```--enable-venv``` option to build a Python venv containing the
-Ocean Model Grid Generator script and all python dependencies.
+The [Ocean Model Grid Generator](https://github.com/NOAA-GFDL/ocean_model_grid_generator) can be copied or cloned from its GFDL homepage.
+
 
 ### User Documentation
 Documentation on using individual tools may be obtained by running
@@ -98,24 +90,13 @@ Contribution requirements include :
  * The passing of existing unit tests when they are run in the CI system.
  * The (potential) addition of unit tests when adding new functionality.
  * For external projects, unit tests may be required and the unit tests
-   should be run on the external CI system. 
+   should be run on the external CI system.
 
  Additionally, since NCTools is distributable via the Spack package manager,
  the NCTools team will need to be able to compile and distribute via Spack any
  submodules and their dependencies. Contributors are encouraged to provide
  Spack recipes for their projects.
 
-## Cloning and submodules
-The NCTools github repository contains Ocean Model Grid Generator's repository
-as a submodule. After cloning NCTools, it must be initialized and updated for
-its submodules:
-
-```
-git clone --recursive https://github.com/NOAA-GFDL/FRE-NCtools
-cd FRE-NCtools
-git submodule update --init --recursive
-```
-
 ## Building and Installation - General Information
 FRE-NCtools has a collection of C and Fortran sources. Within GFDL, FRE-NCtools
 is built using a recent version of the GNU and Intel C and Fortran compilers.
@@ -164,26 +145,13 @@ configure --help=recursive
 
 It is common to compile into a build directory (e.g. named `build`) and
 install into an installation directory (e.g. with full path `<install path>`).
-If the ocean_model_grid_generator is desired, it may be convenient to allow
-the build system to set up a Python venv. These three choices can be
-done with these steps:
-```
-cd FRE-NCtools
-autoreconf -i
-mkdir build && cd build
-../configure --prefix=<install path> --enable-venv
-make
-make install
-```
-
-If the ocean_model_grid_generator is not desired, a similar configuration would
-be achieved by :
+For example:
 
 ```
 cd FRE-NCtools
 autoreconf -i
 mkdir build && cd build
-../configure --prefix=<install path> --disable-ocean-model-grid-generator
+../configure --prefix=<install path>
 make
 make install
 ```

diff --git a/configure.ac b/configure.ac
@@ -22,7 +22,7 @@
 AC_PREREQ([2.63])
 AC_INIT(
   [FRE NC Tools],
-  [2022.02],
+  [2023.01],
   [oar.gfdl.help@noaa.gov],
   [fre-nctools],
   [https://github.com/NOAA-GFDL/FRE-NCtools])
@@ -53,15 +53,6 @@ AC_ARG_ENABLE([quad-precision],
     [],
     [enable_quad_precision=no])
 
-AC_ARG_ENABLE([ocean-model-grid-generator],
-  [AS_HELP_STRING([--enable-ocean-model-grid-generator],
-    [build and install the ocean model grid generator tool])])
-AS_IF([test ${enable_ocean_model_grid_generator:-yes} = yes],
-  [enable_ocean_model_grid_generator=yes],
-  [enable_ocean_model_grid_generator=no])
-AM_CONDITIONAL([ENABLE_OCEAN_MODEL_GRID_GENERATOR],
-  [test $enable_ocean_model_grid_generator = yes])
-
 AC_PROG_CC([icc gcc])
 AM_PROG_CC_C_O
 ## When autoconf v2.70 is more available, this can be replaced with:
@@ -139,6 +130,9 @@ if test "$with_netcdf_fortran" = "no"; then
 fi
 AC_LANG_POP([Fortran])
 
+# Check if Linux sched_getaffinity is available
+AC_CHECK_FUNCS([sched_getaffinity], [], [])
+
 # Check for typedefs, structures, and compiler characteristics
 # Check if the C compiler supports a working `long double` with more range
 # or precision than `double`.
@@ -160,9 +154,6 @@ AC_DEFINE([GIT_REVISION], "git_revision",
 AC_DEFINE([GIT_HEADHASH], "git_hashval",
           [Holds the 'git rev-parse HEAD' information if configure ran within a git working directory])
 
-AS_IF([test $enable_ocean_model_grid_generator = yes],
-  AC_CONFIG_SUBDIRS([tools/ocean_model_grid_generator])
-)
 
 #Code for setting rpath based ncview's configure.in code.
 echo "ac_computer_gnu: $ac_compiler_gnu"

diff --git a/cpp/libfrencutils/create_xgrid.C b/cpp/libfrencutils/create_xgrid.C
@@ -37,6 +37,7 @@
 #include "DITree.h"
 #include "mosaic_util.h"
 #include "create_xgrid.h"
+#include "create_xgrid_aux.h"
 #include "constant.h"
 
 #define AREA_RATIO_THRESH (1.e-6)
@@ -45,6 +46,15 @@
 #define EPSLN30           (1.0e-30)
 #define EPSLN10           (1.0e-10)
 
+using namespace std;
+using namespace nct;
+using BBox_t = nct::BBox3D;
+using BPair_t = nct::BoxAndId;
+using Poly_t = nct::MeshPolygon<double>;
+using Point_t = nct::Point3D<double>;
+using std::vector;
+using std::string;
+
 double grid_box_radius(const double *x, const double *y, const double *z, int n);
 double dist_between_boxes(const double *x1, const double *y1, const double *z1, int n1,
 			  const double *x2, const double *y2, const double *z2, int n2);
@@ -2341,13 +2351,6 @@ int inside_edge(double x0, double y0, double x1, double y1, double x, double y)
 
 // Code below for using search algorithms in create_xgrid_...
 
-using namespace std;
-using namespace nct;
-using BBox_t = nct::BBox3D;
-using BPair_t = nct::BoxAndId;
-using Poly_t = nct::MeshPolygon<double>;
-using Point_t = nct::Point3D<double>;
-
 inline
 size_t pt_idx(const size_t i, const size_t j,  const size_t nx) {
   return ( j * nx + i);
@@ -2436,8 +2439,9 @@ size_t  latlons_outside_ccd_domain(const unsigned int NV4, const double *yv, dou
  * @param debugf
  * @return
  */
+
 BBox_t getBoxForSphericalPolygon(const double lat_m[], const double lon_m[],
-                                  const array<size_t, 4> &is, bool debugf) {
+                                 const array<size_t, 4> &is, bool debugf) {
   constexpr unsigned int NV4{4};
   // xlons are the longitudes that define the X=0 and Y=0 planes(see ll2xyz function)
   // where its possible to have an extrema in X or Y when an edge crosses them.
@@ -2521,6 +2525,7 @@ BBox_t getBoxForSphericalPolygon(const double lat_m[], const double lon_m[],
 }
 
 
+
 /**
  * Generate the exchange grid between two grids for the 2nd order
   * conservative interpolation. This version explicitly uses a search
@@ -2634,6 +2639,7 @@ void  create_xgrid_2dx2d_order2_ws(const int nlon_in, const int nlat_in, const i
           for (auto l = 0; l < n2_in; l++) x2_in[l] -= TPI;
         }
 
+        //Call the 2D_by_2D clipping algorithm
         auto n_out = clip_2dx2d(x1_in, y1_in, n1_in, x2_in,
                                 y2_in, n2_in, x_out, y_out);
         if (n_out > 0) {
@@ -2687,7 +2693,12 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl
 
   }else{
 
-    create_xgrid_2dx2d_order2_ws(nlon_in, nlat_in, nlon_out, nlat_out,
+    /*create_xgrid_2dx2d_order2_ws(nlon_in, nlat_in, nlon_out, nlat_out,
+                                 lon_in, lat_in, lon_out, lat_out, mask_in,
+                                 i_in_r, j_in_r, i_out_r, j_out_r,
+                                 xgrid_area_r, xgrid_clon_r, xgrid_clat_r);*/
+
+    create_xgrid_2dx2d_order2_bfwbb(nlon_in, nlat_in, nlon_out, nlat_out,
                                  lon_in, lat_in, lon_out, lat_out, mask_in,
                                  i_in_r, j_in_r, i_out_r, j_out_r,
                                  xgrid_area_r, xgrid_clon_r, xgrid_clat_r);
@@ -2730,6 +2741,7 @@ int create_xgrid_2dx2d_order2(const int nlon_in, const int nlat_in, const int nl
   }
 }
 
+//TODO: this is to be removed:
 void create_xgrid_2dx2d_order2_check(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out,
                                   const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
                                   const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,

diff --git a/cpp/libfrencutils/create_xgrid.h b/cpp/libfrencutils/create_xgrid.h
@@ -23,14 +23,22 @@
 #define MAXXGRID 1e7
 #endif
 
-#include <array>
-#include "BBox3D.h"
-#define MV 50
 #include <vector>
+#include <array>
 #include <span>
 
+#define MV 50
+
+#include "BBox3D.h"
+using std::vector;
+
+
+
 /* this value is small compare to earth area */
 
+
+
+
 double poly_ctrlon(const double lon[], const double lat[], int n, double clon);
 double poly_ctrlat(const double lon[], const double lat[], int n);
 double box_ctrlon(double ll_lon, double ll_lat, double ur_lon, double ur_lat, double clon);
@@ -79,15 +87,20 @@ int create_xgrid_great_circle(const int *nlon_in, const int *nlat_in, const int
 			      double *xgrid_area, double *xgrid_clon, double *xgrid_clat);
 
 void latlon2xyz(const double lat, const double lon,  std::array<double,3> &  v);
+std::array<size_t, 4>
+get_cell_idxs_ccw_4(const size_t i, const size_t j, const size_t nx);
 nct::BBox3D getBoxForSphericalPolygon(const double lat_m[], const double lon_m[],
-                                 const std::array<size_t, 4> &is, bool debugf  = false);
+                               const std::array<size_t, 4> &is, bool debugf  = false);
 int search_grids(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out,
                  const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
                  const double *mask_in, std::vector<std::vector<size_t>> & results1) ;
-int create_xgrid_2dx2d_order2_ws(const int *nlon_in, const int *nlat_in, const int *nlon_out, const int *nlat_out,
-                              const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
-                              const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out,
-                              double *xgrid_area, double *xgrid_clon, double *xgrid_clat);
+
+void  create_xgrid_2dx2d_order2_ws(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out,
+                                   const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
+                                   const double *mask_in, std::vector<size_t>& i_in, std::vector<size_t>& j_in,
+                                   std::vector<size_t>& i_out, std::vector<size_t>& j_out, std::vector<double>& xgrid_area,
+                                   std::vector<double>& xgrid_clon, std::vector<double>& xgrid_clat);
+
 void create_xgrid_2dx2d_order2_ws_check(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out,
                                         const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
                                         const double *mask_in, int *i_in, int *j_in, int *i_out, int *j_out, int nxgrid1);
@@ -103,5 +116,11 @@ void create_xgrid_2dx2d_order2_check(const int nlon_in, const int nlat_in, const
 template<class T>
 void printPolygon(std::ostream &os, std::span<T> lonv, std::span<T> latv) ;
 
+void  create_xgrid_2dx2d_order2_bfwbb(const int nlon_in, const int nlat_in, const int nlon_out, const int nlat_out,
+                                   const double *lon_in, const double *lat_in, const double *lon_out, const double *lat_out,
+                                   const double *mask_in, std::vector<size_t>& i_in, std::vector<size_t>& j_in,
+                                   std::vector<size_t>& i_out, std::vector<size_t>& j_out, std::vector<double>& xgrid_area,
+                                   std::vector<double>& xgrid_clon, std::vector<double>& xgrid_clat);
+
 
 #endif
diff --git a/cpp/libfrencutils/create_xgrid_aux.h b/cpp/libfrencutils/create_xgrid_aux.h
@@ -8,10 +8,17 @@
 #include <algorithm>
 #include <span>
 #include <source_location>
+
+#include "constant.h"
+#include "mpp.h"
+#include "create_xgrid.h"
+
+
 #include "BBox3D.h"
+#include "BoxedObj.h"
 #include "Polygon.h"
 #include "mosaic_util.h"
-#include "create_xgrid.h"
+
 
 bool checkBBoxViaPolySamples(std::span<double> yv, std::span<double> xv,
                               nct::BBox3D & box, unsigned int npoints1D = 5 ) {
@@ -35,11 +42,11 @@ bool checkBBoxViaPolySamples(std::span<double> yv, std::span<double> xv,
       auto contains_point = nct::BBox3D::contains(box, pt);
       if (!contains_point) {
         passed = false;
-        //TODO: the polygon and the box can be saved to strstream
-         printPolygon<double>(std::cout, xv, yv);
-         std::cout << box << std::endl;
-         auto str = std::format("< {:16.10e}, {:16.10e}, {:16.10e}>", pt[0],pt[1],pt[2]);
-         std::cout << str <<std::endl;
+        //TODO: DEBUG or CTest option ? : saved to strstream the polygon and the box ?
+       //  printPolygon<double>(std::cout, xv, yv);
+       //  std::cout << box << std::endl;
+      //   auto str = std::format("< {:16.10e}, {:16.10e}, {:16.10e}>", pt[0],pt[1],pt[2]);
+      //   std::cout << str <<std::endl;
       }
     }
   }