C and Fortran API for Interpolator

src/pcms/capi/interpolator.cpp

@@ -27,6 +27,73 @@ PcmsPointBasedInterpolatorHandle pcms_create_point_based_interpolator(void* sour
   return {reinterpret_cast<void*>(interpolator)};
 }
 
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgc_interpolator(const char* xgc_mesh_filename, const char* dg2_mesh_filename, double radius)
+{
+  // use the xgc_mesh_nodes as source points and
+  // dg2 element centroids as target points
+  // then, create a point based interpolator like above
+  auto xgc_fname = std::string(xgc_mesh_filename);
+  auto dg2_fname = std::string(dg2_mesh_filename);
+
+  // trim the filenames since they are coming from c or fortran api which may have extra spaces at the end
+  xgc_fname = xgc_fname.erase(xgc_fname.find_last_not_of(" \n\r\t")+1);
+  dg2_fname = dg2_fname.erase(dg2_fname.find_last_not_of(" \n\r\t")+1);
+  printf("The interpolator got xgc mesh file: %s and dg2 mesh file: %s\n", xgc_fname.c_str(), dg2_fname.c_str());
+
+  // read the meshes
+  auto xgc_mesh_lib = Omega_h::Library();
+  auto xgc_mesh = Omega_h::binary::read(xgc_fname, xgc_mesh_lib.world());
+  auto dg2_mesh_lib = Omega_h::Library();
+  auto dg2_mesh = Omega_h::binary::read(dg2_fname, dg2_mesh_lib.world());
+
+  auto xgc_nodes = xgc_mesh.coords();
+  auto xgc_num_nodes = xgc_mesh.nverts();
+  OMEGA_H_CHECK_PRINTF(xgc_mesh.dim() == 2, "XGC mesh dimension is not 2D %d\n", xgc_mesh.dim());
+  OMEGA_H_CHECK_PRINTF(xgc_num_nodes * 2 == xgc_nodes.size(),
+    "XGC mesh nodes size does not match the number of vertices %d != %d\n", xgc_num_nodes * 2, xgc_nodes.size());
+  OMEGA_H_CHECK_PRINTF(dg2_mesh.dim() == 2, "DG2 mesh dimension is not 2D %d\n", dg2_mesh.dim());
+
+  auto dg2_num_elems = dg2_mesh.nelems();
+  auto dg2_elem_centroids = getCentroids(dg2_mesh);
+  OMEGA_H_CHECK_PRINTF(dg2_num_elems * 2 == dg2_elem_centroids.size(),
+    "DG2 mesh element centroids size does not match the number of elements %d != %d\n", dg2_num_elems * 2, dg2_elem_centroids.size());
+
+  Omega_h::HostRead<Omega_h::Real> xgc_nodes_host(xgc_nodes);
+  Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host(dg2_elem_centroids);
+
+  return pcms_create_point_based_interpolator(
+    (void*) dg2_elem_centroids_host.data(), dg2_elem_centroids.size(),
+        (void*) xgc_nodes_host.data(), xgc_nodes.size(),
+    radius);
+}
+
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgcnode_interpolator(void* target_points, int target_points_size,
+                                                                const char* dg2_mesh_filename, double radius)
+{
+  // same as above pcms_create_degas2xgc_interpolator but the target points are provided by the user
+  // this is useful when the corresponding xgc mesh is not available
+
+  auto dg2_fname = std::string(dg2_mesh_filename);
+  dg2_fname = dg2_fname.erase(dg2_fname.find_last_not_of(" \n\r\t")+1);
+  printf("The interpolator got dg2 mesh file: %s\n", dg2_fname.c_str());
+  auto dg2_mesh_lib = Omega_h::Library(nullptr, nullptr, MPI_COMM_SELF);
+  auto dg2_mesh = Omega_h::binary::read(dg2_fname, dg2_mesh_lib.world());
+
+  OMEGA_H_CHECK_PRINTF(dg2_mesh.dim() == 2, "DG2 mesh dimension is not 2D %d\n", dg2_mesh.dim());
+  auto dg2_num_elems = dg2_mesh.nelems();
+  auto dg2_elem_centroids = getCentroids(dg2_mesh);
+  OMEGA_H_CHECK_PRINTF(dg2_num_elems * 2 == dg2_elem_centroids.size(),
+    "DG2 mesh element centroids size does not match the number of elements %d != %d\n", dg2_num_elems * 2, dg2_elem_centroids.size());
+
+  Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host(dg2_elem_centroids);
+  printf("Number of DG2 element centroids: %d\n", dg2_elem_centroids.size()/2);
+
+  return pcms_create_point_based_interpolator(
+    (void*) dg2_elem_centroids_host.data(), dg2_elem_centroids_host.size(),
+    target_points, target_points_size,
+    radius);
+}
+
 void pcms_destroy_point_based_interpolator(PcmsPointBasedInterpolatorHandle interpolator)
 {
   if (interpolator.pointer != nullptr) {

src/pcms/capi/interpolator.h

@@ -30,6 +30,9 @@ enum for interpolation type
 PcmsInterpolatorHandle pcms_create_interpolator(PcmsInterpolatorOHMeshHandle oh_mesh, double radius);
 PcmsPointBasedInterpolatorHandle pcms_create_point_based_interpolator(void* source_points, int source_points_size,
                                                             void* target_points, int target_points_size, double radius);
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgc_interpolator(const char* xgc_mesh_filename, const char* dg2_mesh_filename, double radius);
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgcnode_interpolator(void* target_points, int target_points_size,
+                                                                    const char* dg2_mesh_filename, double radius);
 void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator);
 void pcms_destroy_point_based_interpolator(PcmsPointBasedInterpolatorHandle interpolator);
 

src/pcms/fortranapi/interpolator.i

@@ -26,6 +26,9 @@ typedef struct PcmsPointBasedInterpolatorHandle PcmsPointBasedInterpolatorHandle
 
 PcmsPointBasedInterpolatorHandle pcms_create_point_based_interpolator(void* source_points, int source_points_size,
                                                                       void* target_points, int target_points_size, double radius);
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgc_interpolator(const char* xgc_mesh_filename, const char* dg2_mesh_filename, double radius);
+PcmsPointBasedInterpolatorHandle pcms_create_degas2xgcnode_interpolator(void* target_points, int target_points_size,
+                                                                const char* dg2_mesh_filename, double radius);
 PcmsInterpolatorHandle pcms_create_interpolator(PcmsInterpolatorOHMeshHandle oh_mesh, double radius);
 void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator);
 void pcms_destroy_point_based_interpolator(PcmsPointBasedInterpolatorHandle interpolator);

src/pcms/fortranapi/interpolator_wrap.c

@@ -503,6 +503,42 @@ SWIGEXPORT SwigClassWrapper _wrap_pcms_create_point_based_interpolator(void cons
 }
 
 
+SWIGEXPORT SwigClassWrapper _wrap_pcms_create_degas2xgc_interpolator(SwigArrayWrapper *farg1, SwigArrayWrapper *farg2, double const *farg3) {
+  SwigClassWrapper fresult ;
+  char *arg1 = (char *) 0 ;
+  char *arg2 = (char *) 0 ;
+  double arg3 ;
+  PcmsPointBasedInterpolatorHandle result;
+
+  arg1 = (char *)(farg1->data);
+  arg2 = (char *)(farg2->data);
+  arg3 = (double)(*farg3);
+  result = pcms_create_degas2xgc_interpolator((char const *)arg1,(char const *)arg2,arg3);
+  fresult.cptr = (PcmsPointBasedInterpolatorHandle *)memcpy((PcmsPointBasedInterpolatorHandle *)calloc(1,sizeof(PcmsPointBasedInterpolatorHandle)),&result,sizeof(PcmsPointBasedInterpolatorHandle));
+  fresult.cmemflags = SWIG_MEM_RVALUE | SWIG_MEM_OWN;
+  return fresult;
+}
+
+
+SWIGEXPORT SwigClassWrapper _wrap_pcms_create_degas2xgcnode_interpolator(void const **farg1, int const *farg2, SwigArrayWrapper *farg3, double const *farg4) {
+  SwigClassWrapper fresult ;
+  void *arg1 = (void *) 0 ;
+  int arg2 ;
+  char *arg3 = (char *) 0 ;
+  double arg4 ;
+  PcmsPointBasedInterpolatorHandle result;
+
+  arg1 = (void *)(*farg1);
+  arg2 = (int)(*farg2);
+  arg3 = (char *)(farg3->data);
+  arg4 = (double)(*farg4);
+  result = pcms_create_degas2xgcnode_interpolator(arg1,arg2,(char const *)arg3,arg4);
+  fresult.cptr = (PcmsPointBasedInterpolatorHandle *)memcpy((PcmsPointBasedInterpolatorHandle *)calloc(1,sizeof(PcmsPointBasedInterpolatorHandle)),&result,sizeof(PcmsPointBasedInterpolatorHandle));
+  fresult.cmemflags = SWIG_MEM_RVALUE | SWIG_MEM_OWN;
+  return fresult;
+}
+
+
 SWIGEXPORT SwigClassWrapper _wrap_pcms_create_interpolator(SwigClassWrapper *farg1, double const *farg2) {
   SwigClassWrapper fresult ;
   PcmsInterpolatorOHMeshHandle arg1 ;

src/pcms/fortranapi/pcms_interpolator.f90

@@ -49,15 +49,17 @@ module pcms_interpolator
   generic :: assignment(=) => swigf_PcmsPointBasedInterpolatorHandle_op_assign__
  end type PcmsPointBasedInterpolatorHandle
  public :: pcms_create_point_based_interpolator
+ type, bind(C) :: SwigArrayWrapper
+  type(C_PTR), public :: data = C_NULL_PTR
+  integer(C_SIZE_T), public :: size = 0
+ end type
+ public :: pcms_create_degas2xgc_interpolator
+ public :: pcms_create_degas2xgcnode_interpolator
  public :: pcms_create_interpolator
  public :: pcms_destroy_interpolator
  public :: pcms_destroy_point_based_interpolator
  public :: pcms_kokkos_initialize_without_args
  public :: pcms_kokkos_finalize
- type, bind(C) :: SwigArrayWrapper
-  type(C_PTR), public :: data = C_NULL_PTR
-  integer(C_SIZE_T), public :: size = 0
- end type
  public :: read_oh_mesh
  public :: release_oh_mesh
  public :: pcms_interpolate
@@ -224,6 +226,31 @@ function swigc_pcms_create_point_based_interpolator(farg1, farg2, farg3, farg4,
 type(SwigClassWrapper) :: fresult
 end function
 
+function swigc_pcms_create_degas2xgc_interpolator(farg1, farg2, farg3) &
+bind(C, name="_wrap_pcms_create_degas2xgc_interpolator") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigarraywrapper
+import :: swigclasswrapper
+type(SwigArrayWrapper) :: farg1
+type(SwigArrayWrapper) :: farg2
+real(C_DOUBLE), intent(in) :: farg3
+type(SwigClassWrapper) :: fresult
+end function
+
+function swigc_pcms_create_degas2xgcnode_interpolator(farg1, farg2, farg3, farg4) &
+bind(C, name="_wrap_pcms_create_degas2xgcnode_interpolator") &
+result(fresult)
+use, intrinsic :: ISO_C_BINDING
+import :: swigarraywrapper
+import :: swigclasswrapper
+type(C_PTR), intent(in) :: farg1
+integer(C_INT), intent(in) :: farg2
+type(SwigArrayWrapper) :: farg3
+real(C_DOUBLE), intent(in) :: farg4
+type(SwigClassWrapper) :: fresult
+end function
+
 function swigc_pcms_create_interpolator(farg1, farg2) &
 bind(C, name="_wrap_pcms_create_interpolator") &
 result(fresult)
@@ -538,6 +565,68 @@ function pcms_create_point_based_interpolator(source_points, source_points_size,
 swig_result%swigdata = fresult
 end function
 
+
+subroutine SWIGTM_fin_char_Sm_(finp, iminp, temp)
+  use, intrinsic :: ISO_C_BINDING
+  character(len=*), intent(in) :: finp
+  type(SwigArrayWrapper), intent(out) :: iminp
+  character(kind=C_CHAR), dimension(:), target, allocatable, intent(out) :: temp
+  integer :: i
+
+  allocate(character(kind=C_CHAR) :: temp(len(finp) + 1))
+  do i=1,len(finp)
+    temp(i) = char(ichar(finp(i:i)), kind=C_CHAR)
+  end do
+  i = len(finp) + 1
+  temp(i) = C_NULL_CHAR ! C finp compatibility
+  iminp%data = c_loc(temp)
+  iminp%size = len(finp, kind=C_SIZE_T)
+end subroutine
+
+function pcms_create_degas2xgc_interpolator(xgc_mesh_filename, dg2_mesh_filename, radius) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+type(PcmsPointBasedInterpolatorHandle) :: swig_result
+character(len=*), intent(in) :: xgc_mesh_filename
+character(len=*), intent(in) :: dg2_mesh_filename
+real(C_DOUBLE), intent(in) :: radius
+type(SwigClassWrapper) :: fresult 
+character(kind=C_CHAR), dimension(:), allocatable, target :: farg1_temp 
+type(SwigArrayWrapper) :: farg1 
+character(kind=C_CHAR), dimension(:), allocatable, target :: farg2_temp 
+type(SwigArrayWrapper) :: farg2 
+real(C_DOUBLE) :: farg3 
+
+call SWIGTM_fin_char_Sm_(xgc_mesh_filename, farg1, farg1_temp)
+call SWIGTM_fin_char_Sm_(dg2_mesh_filename, farg2, farg2_temp)
+farg3 = radius
+fresult = swigc_pcms_create_degas2xgc_interpolator(farg1, farg2, farg3)
+swig_result%swigdata = fresult
+end function
+
+function pcms_create_degas2xgcnode_interpolator(target_points, target_points_size, dg2_mesh_filename, radius) &
+result(swig_result)
+use, intrinsic :: ISO_C_BINDING
+type(PcmsPointBasedInterpolatorHandle) :: swig_result
+type(C_PTR), intent(in) :: target_points
+integer(C_INT), intent(in) :: target_points_size
+character(len=*), intent(in) :: dg2_mesh_filename
+real(C_DOUBLE), intent(in) :: radius
+type(SwigClassWrapper) :: fresult 
+type(C_PTR) :: farg1 
+integer(C_INT) :: farg2 
+character(kind=C_CHAR), dimension(:), allocatable, target :: farg3_temp 
+type(SwigArrayWrapper) :: farg3 
+real(C_DOUBLE) :: farg4 
+
+farg1 = target_points
+farg2 = target_points_size
+call SWIGTM_fin_char_Sm_(dg2_mesh_filename, farg3, farg3_temp)
+farg4 = radius
+fresult = swigc_pcms_create_degas2xgcnode_interpolator(farg1, farg2, farg3, farg4)
+swig_result%swigdata = fresult
+end function
+
 function pcms_create_interpolator(oh_mesh, radius) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING
@@ -584,24 +673,6 @@ subroutine pcms_kokkos_finalize()
 call swigc_pcms_kokkos_finalize()
 end subroutine
 
-
-subroutine SWIGTM_fin_char_Sm_(finp, iminp, temp)
-  use, intrinsic :: ISO_C_BINDING
-  character(len=*), intent(in) :: finp
-  type(SwigArrayWrapper), intent(out) :: iminp
-  character(kind=C_CHAR), dimension(:), target, allocatable, intent(out) :: temp
-  integer :: i
-
-  allocate(character(kind=C_CHAR) :: temp(len(finp) + 1))
-  do i=1,len(finp)
-    temp(i) = char(ichar(finp(i:i)), kind=C_CHAR)
-  end do
-  i = len(finp) + 1
-  temp(i) = C_NULL_CHAR ! C finp compatibility
-  iminp%data = c_loc(temp)
-  iminp%size = len(finp, kind=C_SIZE_T)
-end subroutine
-
 function read_oh_mesh(filename) &
 result(swig_result)
 use, intrinsic :: ISO_C_BINDING

src/pcms/interpolator/interpolation_base.cpp

@@ -153,6 +153,8 @@ void MLSPointCloudInterpolation::find_supports(uint min_req_supports)
   uint min_supports_found = 0;
   uint max_supports_found = 0;
   // radius adjustment loop
+  int loop_count = 0;
+  int max_count = 100;
   while (min_supports_found < min_req_supports ||
          max_supports_found > 3 * min_req_supports) {
     // n^2 search, compare each point with all other points
@@ -184,6 +186,15 @@ void MLSPointCloudInterpolation::find_supports(uint min_req_supports)
       break;
     }
 
+    loop_count++;
+    Kokkos::fence();
+    if (loop_count > 100) {
+      printf("Loop count exceeded 100 and still not converged.\n"
+        "Manually check if the number of minimum and maximum supports are reasonable.\n"
+        "There are situations when it may not converge.\n");
+      break;
+    }
+
     Kokkos::Min<uint> min_reducer(min_supports_found);
     Kokkos::parallel_reduce(
       "find number of supports", num_supports.size(),
@@ -203,8 +214,8 @@ void MLSPointCloudInterpolation::find_supports(uint min_req_supports)
 
     // increase radius if not enough supports or too many supports
     if (min_supports_found < min_req_supports || max_supports_found > 3 * min_req_supports) {
-      printf("Adjusting radius:(min: %d max: %d) min_req_supports: %d\n",
-             min_supports_found, max_supports_found, min_req_supports);
+      printf("Adjusting radius iter %d:(min: %d max: %d) min_req_supports: %d\n",
+             loop_count, min_supports_found, max_supports_found, min_req_supports);
 
       Kokkos::fence();
       Omega_h::parallel_for(
@@ -227,7 +238,7 @@ void MLSPointCloudInterpolation::find_supports(uint min_req_supports)
         });
     }
   }
-  printf("Supports found: min: %d max: %d\n", min_supports_found, max_supports_found);
+  printf("Searched %d times and supports found: min: %d max: %d\n", loop_count, min_supports_found, max_supports_found);
 
   // parallel scan for fill the support index with cumulative sum
   auto support_ptr_l = Omega_h::Write<Omega_h::LO>(n_targets + 1, 0);