C and Fortran API for Interpolator

src/pcms/capi/CMakeLists.txt

@@ -11,8 +11,21 @@ set(CAPI_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/client.h ${CMAKE_CURRENT_SOURCE_DIR
 set_target_properties(pcms_capi PROPERTIES PUBLIC_HEADER "${CAPI_HEADERS}"
          OUTPUT_NAME pcmscapi
          EXPORT_NAME capi)
+
+add_library(pcms_interpolator_capi interpolator.cpp kokkos.cpp)
+add_library(pcms::interpolator::capi ALIAS pcms_interpolator_capi)
+target_link_libraries(pcms_interpolator_capi PUBLIC MPI::MPI_C PRIVATE pcms::interpolator)
+target_include_directories( pcms_interpolator_capi
+        PUBLIC
+        "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/../..>" # this makes the module path cpms/capi
+        "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}/pcms/capi>")
+set(INT_CAPI_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/interpolator.h)
+set_target_properties(pcms_interpolator_capi PROPERTIES PUBLIC_HEADERS "${INT_CAPI_HEADERS}"
+        OUTPUTNAME pcmsinterpolatorcapi
+        EXPORT_NAME interpolatorcapi)
+
 install(
-        TARGETS pcms_capi
+        TARGETS pcms_capi pcms_interpolator_capi
         EXPORT pcms_capi-targets
         LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
         ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
@@ -23,3 +36,4 @@ install(
         EXPORT pcms_capi-targets
         NAMESPACE pcms::
         DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/pcms)
+

src/pcms/capi/interpolator.cpp

@@ -0,0 +1,151 @@
+//
+// Created by hasanm4 on 2/17/25.
+//
+#include <pcms/capi/kokkos.h>
+#include "interpolator.h"
+#include <interpolation_base.h>
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_mesh.hpp>
+
+
+//[[nodiscard]]
+PcmsInterpolatorHandle pcms_create_interpolator(PcmsInterpolatorOHMeshHandle oh_mesh, double radius)
+{
+  auto* source_mesh = reinterpret_cast<Omega_h::Mesh*>(oh_mesh.mesh_handle);
+  auto* interpolator = new MLSInterpolationHandler(*source_mesh, radius);
+  return {reinterpret_cast<void*>(interpolator)};
+}
+
+
+PcmsPointBasedInterpolatorHandle pcms_create_point_based_interpolator(void* source_points, int source_points_size,
+                                                            void* target_points, int target_points_size, double radius) {
+
+  auto source_points_view = pcms::ScalarArrayView<double, pcms::HostMemorySpace>(reinterpret_cast<double*>(source_points), source_points_size);
+  auto target_points_view = pcms::ScalarArrayView<double, pcms::HostMemorySpace>(reinterpret_cast<double*>(target_points), target_points_size);
+  auto* interpolator = new MLSPointCloudInterpolation(source_points_view, target_points_view, 2, radius, 12, 3, true);
+  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) {
+    delete reinterpret_cast<MLSPointCloudInterpolation*>(interpolator.pointer);
+  }
+}
+
+void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator)
+{
+  if (interpolator.pointer != nullptr) {
+    delete reinterpret_cast<MLSInterpolationHandler*>(interpolator.pointer);
+  }
+}
+
+PcmsInterpolatorOHMeshHandle read_oh_mesh(const char * filename){
+  auto fname = std::string(filename);
+  // trim the filename since it is coming from c or fortran api which may have extra spaces at the end
+  fname.erase(fname.find_last_not_of(" \n\r\t")+1);
+  auto* oh_lib = new Omega_h::Library();
+  auto* mesh = new Omega_h::Mesh(Omega_h::binary::read(fname, oh_lib->world()));
+
+  return {reinterpret_cast<void*>(mesh), reinterpret_cast<void*>(oh_lib)};
+}
+
+void release_oh_mesh(PcmsInterpolatorOHMeshHandle oh_mesh){
+  if (oh_mesh.mesh_handle != nullptr) {
+    assert(oh_mesh.lib_handle != nullptr);
+    delete reinterpret_cast<Omega_h::Mesh*>(oh_mesh.mesh_handle);
+    delete reinterpret_cast<Omega_h::Library*>(oh_mesh.lib_handle);
+  }
+}
+
+void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size) {
+  auto* mls_interpolator = reinterpret_cast<MLSInterpolationHandler*>(interpolator.pointer);
+
+  OMEGA_H_CHECK_PRINTF(input_size == mls_interpolator->getSourceSize(), "Input array size does not match the source size %d != %d\n", input_size, mls_interpolator->getSourceSize());
+  OMEGA_H_CHECK_PRINTF(output_size == mls_interpolator->getTargetSize(), "Output array size does not match the target size %d != %d\n", output_size, mls_interpolator->getTargetSize());
+
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> input_array(reinterpret_cast<double*>(input), input_size);
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> output_array(reinterpret_cast<double*> (output), output_size);
+
+  mls_interpolator->eval(input_array, output_array);
+}
+
+void pcms_interpolate_point_based(PcmsPointBasedInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size) {
+  auto* mls_interpolator = reinterpret_cast<MLSPointCloudInterpolation*>(interpolator.pointer);
+
+  OMEGA_H_CHECK_PRINTF(input_size == mls_interpolator->getSourceSize(), "Input array size does not match the source size %d != %d\n", input_size, mls_interpolator->getSourceSize());
+  OMEGA_H_CHECK_PRINTF(output_size == mls_interpolator->getTargetSize(), "Output array size does not match the target size %d != %d\n", output_size, mls_interpolator->getTargetSize());
+
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> input_array(reinterpret_cast<double*>(input), input_size);
+  pcms::ScalarArrayView<double, pcms::HostMemorySpace> output_array(reinterpret_cast<double*> (output), output_size);
+
+  mls_interpolator->eval(input_array, output_array);
+}

src/pcms/capi/interpolator.h

@@ -0,0 +1,52 @@
+//
+// Created by hasanm4 on 2/17/25.
+//
+
+#ifndef PCMS_INTERPOLATOR_CAPI_H
+#define PCMS_INTERPOLATOR_CAPI_H
+
+#include <stdint.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+struct PcmsInterpolatorHandle { void* pointer; };
+typedef struct PcmsInterpolatorHandle PcmsInterpolatorHandle;
+
+struct PcmsPointBasedInterpolatorHandle {void * pointer;};
+typedef struct PcmsPointBasedInterpolatorHandle PcmsPointBasedInterpolatorHandle;
+
+struct PcmsInterpolatorOHMeshHandle {
+  void* mesh_handle;
+  void* lib_handle;
+};
+typedef struct PcmsInterpolatorOHMeshHandle PcmsInterpolatorOHMeshHandle;
+
+/*
+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);
+
+PcmsInterpolatorOHMeshHandle read_oh_mesh(const char* filename);
+void release_oh_mesh(PcmsInterpolatorOHMeshHandle oh_mesh_handle);
+
+
+void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size);
+void pcms_interpolate_point_based(PcmsPointBasedInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size);
+
+
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // PCMS_INTERPOLATOR_CAPI_H

src/pcms/fortranapi/CMakeLists.txt

@@ -7,8 +7,22 @@ set_target_properties(pcms_fortranapi PROPERTIES Fortran_MODULE_DIRECTORY ${CMAK
 target_include_directories(pcms_fortranapi PUBLIC $<BUILD_INTERFACE:$<TARGET_PROPERTY:pcms::fortranapi,Fortran_MODULE_DIRECTORY>>
                                                     $<INSTALL_INTERFACE:${CMAKE_INSTALL_LIBDIR}/fortran>)
 target_link_libraries(pcms_fortranapi PRIVATE pcms::capi PUBLIC Kokkos::kokkos)
+
+add_library(pcms_interpolator_fortranapi pcms_interpolator.f90 interpolator_wrap.c)
+add_library(pcms::interpolator::fortranapi ALIAS pcms_interpolator_fortranapi)
+set_target_properties(pcms_interpolator_fortranapi PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/fortran
+        OUTPUT_NAME pcmsinterpolatorfortranapi
+        EXPORT_NAME interpolatorfortranapi)
+target_include_directories(pcms_interpolator_fortranapi PUBLIC $<BUILD_INTERFACE:$<TARGET_PROPERTY:pcms::interpolator::fortranapi,Fortran_MODULE_DIRECTORY>>
+        $<INSTALL_INTERFACE:${CMAKE_INSTALL_LIBDIR}/fortran>)
+target_link_libraries(pcms_interpolator_fortranapi PRIVATE pcms::interpolator::capi PUBLIC Kokkos::kokkos)
+
 INSTALL(DIRECTORY $<TARGET_PROPERTY:pcms::fortranapi,Fortran_MODULE_DIRECTORY>
         DESTINATION ${CMAKE_INSTALL_LIBDIR})
+
+INSTALL(DIRECTORY $<TARGET_PROPERTY:pcms::interpolator::fortranapi,Fortran_MODULE_DIRECTORY>
+        DESTINATION ${CMAKE_INSTALL_LIBDIR})
+
 install(
         TARGETS pcms_fortranapi
         EXPORT pcms_fortranapi-targets
@@ -17,6 +31,16 @@ install(
         RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
         INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
         PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pcms/fortranapi/)
+
+install(
+        TARGETS pcms_interpolator_fortranapi
+        EXPORT pcms_fortranapi-targets
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+        INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
+        PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/pcms/fortranapi/)
+
 install(
         EXPORT pcms_fortranapi-targets
         NAMESPACE pcms::

src/pcms/fortranapi/interpolator.i

@@ -0,0 +1,44 @@
+%module pcms_interpolator
+%{
+#include "pcms/capi/interpolator.h"
+#include "pcms/capi/kokkos.h"
+%}
+%include <../external/flibhpc/include/mpi.i>
+%include <stdint.i>
+%include <typemaps.i>
+
+struct PcmsInterpolatorOHMeshHandle
+{
+  void* mesh_handle;
+  void* lib_handle;
+};
+typedef struct PcmsInterpolatorOHMeshHandle PcmsInterpolatorOHMeshHandle;
+
+struct PcmsInterpolatorHandle {
+  void* pointer;
+};
+typedef struct PcmsInterpolatorHandle PcmsInterpolatorHandle;
+
+struct PcmsPointBasedInterpolatorHandle {void * pointer;};
+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);
+
+void pcms_kokkos_initialize_without_args();
+void pcms_kokkos_finalize();
+
+PcmsInterpolatorOHMeshHandle read_oh_mesh(const char* filename);
+void release_oh_mesh(PcmsInterpolatorOHMeshHandle oh_mesh);
+
+void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size);
+void pcms_interpolate_point_based(PcmsPointBasedInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size);
+