C and Fortran API for Interpolator

src/pcms/capi/CMakeLists.txt

@@ -16,9 +16,24 @@ set(CAPI_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/client.h
                  ${CMAKE_CURRENT_SOURCE_DIR}/kokkos.h)
 set_target_properties(
   pcms_capi PROPERTIES PUBLIC_HEADER "${CAPI_HEADERS}" OUTPUT_NAME pcmscapi
-                       EXPORT_NAME capi)
+             EXPORT_NAME capi)
+
+add_library(pcms_interpolator_capi mesh.cpp 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
+                     ${CMAKE_CURRENT_SOURCE_DIR}/mesh.h)
+set_target_properties(pcms_interpolator_capi
+  PROPERTIES PUBLIC_HEADERS "${INT_CAPI_HEADERS}"
+             OUTPUT_NAME 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}

src/pcms/capi/interpolator.cpp

@@ -0,0 +1,134 @@
+//
+// Created by hasanm4 on 2/17/25.
+//
+#include <pcms/capi/kokkos.h>
+#include <pcms/capi/interpolator.h>
+#include <pcms/interpolator/interpolation_base.h>
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_mesh.hpp>
+#include <pcms/print.h>
+
+//[[nodiscard]]
+PcmsInterpolatorHandle pcms_create_interpolator(PcmsOmegaHMeshHandle oh_mesh,
+                                                double radius)
+{
+  auto* source_mesh = reinterpret_cast<Omega_h::Mesh*>(oh_mesh.mesh_handle);
+  auto* interpolator = new MLSMeshInterpolation(*source_mesh, radius);
+  return {reinterpret_cast<void*>(interpolator)};
+}
+
+PcmsInterpolatorHandle pcms_create_point_based_interpolator(
+  void* source_points, int source_points_size, void* target_points,
+  int target_points_size, double radius, int degree, int min_req_supports,
+  double lambda, double decay_factor)
+{
+
+  auto source_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
+    reinterpret_cast<double*>(source_points), source_points_size);
+  auto target_points_view = pcms::Rank1View<double, pcms::HostMemorySpace>(
+    reinterpret_cast<double*>(target_points), target_points_size);
+  auto* interpolator = new MLSPointCloudInterpolation(
+    source_points_view, target_points_view, 2, radius, min_req_supports, degree,
+    true, lambda, decay_factor);
+  return {reinterpret_cast<void*>(interpolator)};
+}
+
+Omega_h::HostRead<Omega_h::Real> read_mesh_centroids(const char* mesh_filename,
+                                                     int& num_elements)
+{
+  auto fname = std::string(mesh_filename);
+  fname = fname.erase(fname.find_last_not_of(" \n\r\t") + 1);
+  pcms::printInfo("The interpolator got dg2 mesh file: %s\n", fname.c_str());
+  auto mesh_lib = Omega_h::Library(nullptr, nullptr, MPI_COMM_SELF);
+  auto mesh = Omega_h::binary::read(fname, mesh_lib.world());
+  auto elem_centroids = getCentroids(mesh);
+  num_elements = mesh.nelems();
+  OMEGA_H_CHECK_PRINTF(num_elements * 2 == elem_centroids.size(),
+                       "Mesh element centroids size does not match the number "
+                       "of elements %d != %d\n",
+                       num_elements * 2, elem_centroids.size());
+
+  pcms::printInfo("Number of element centroids: %d\n",
+                  elem_centroids.size() / 2);
+  OMEGA_H_CHECK_PRINTF(mesh.dim() == 2, "Mesh dimension is not 2D %d\n",
+                       mesh.dim());
+
+  return {elem_centroids};
+}
+
+void write_void_int_pointer(void* pointer, int value)
+{
+  if (pointer) {
+    int* dg2_elem_count_int = reinterpret_cast<int*>(pointer);
+    *dg2_elem_count_int = value;
+  } else {
+    pcms::printError("Error: NULL pointer provided to write integer value\n");
+  }
+}
+
+PcmsInterpolatorHandle pcms_create_degas2xgcnode_interpolator(
+  void* target_points, int target_points_size, const char* dg2_mesh_filename,
+  double radius, void* dg2_elem_count, int degree, int min_req_supports,
+  double lambda, double decay_factor)
+{
+  // 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
+
+  int dg2_num_elems = 0;
+  Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host =
+    read_mesh_centroids(dg2_mesh_filename, dg2_num_elems);
+  write_void_int_pointer(dg2_elem_count, dg2_num_elems);
+
+  return pcms_create_point_based_interpolator(
+    (void*)dg2_elem_centroids_host.data(), dg2_elem_centroids_host.size(),
+    target_points, target_points_size, radius, degree, min_req_supports, lambda,
+    decay_factor);
+}
+
+PcmsInterpolatorHandle pcms_create_xgcnodedegas2_interpolator(
+  const char* dg2_mesh_filename, void* source_points, int source_points_size,
+  double radius, void* dg2_elem_count, int degree, int min_req_supports,
+  double lambda, double decay_factor)
+{
+  int dg2_num_elems = 0;
+  Omega_h::HostRead<Omega_h::Real> dg2_elem_centroids_host =
+    read_mesh_centroids(dg2_mesh_filename, dg2_num_elems);
+  write_void_int_pointer(dg2_elem_count, dg2_num_elems);
+
+  return pcms_create_point_based_interpolator(
+    source_points, source_points_size, (void*)dg2_elem_centroids_host.data(),
+    dg2_elem_centroids_host.size(), radius, degree, min_req_supports, lambda,
+    decay_factor);
+}
+
+void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator)
+{
+  if (interpolator.pointer != nullptr) {
+    delete reinterpret_cast<InterpolationBase*>(interpolator.pointer);
+  }
+}
+
+void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input,
+                      int input_size, void* output, int output_size)
+{
+  auto* mls_interpolator =
+    reinterpret_cast<InterpolationBase*>(interpolator.pointer);
+
+  OMEGA_H_CHECK_PRINTF(
+    input_size == mls_interpolator->getSourceSize(),
+    "Input array size does not match the source size %d != %zu\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 != %zu\n", output_size,
+    mls_interpolator->getTargetSize());
+
+  pcms::Rank1View<double, pcms::HostMemorySpace> input_array(
+    reinterpret_cast<double*>(input), input_size);
+  pcms::Rank1View<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,182 @@
+/**
+ * @file interpolator.h
+ * @brief PCMS Interpolator C API Header
+ *
+ * This header defines the C API for the PCMS Interpolator library, which
+ * provides functionality for creating and using interpolators based on Moving
+ * Least Squares (MLS) methods. For now, it only supports 2D interpolation and
+ * hard-coded for RBF, Gaussian weight function.
+ */
+
+#ifndef PCMS_INTERPOLATOR_CAPI_H
+#define PCMS_INTERPOLATOR_CAPI_H
+
+#include <stdint.h>
+#include "mesh.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * @brief Holds a void pointer of InterpolationBase object
+ * @see InterpolationBase, MLSMeshInterpolation, MLSPointCloudInterpolation
+ */
+struct PcmsInterpolatorHandle
+{
+  void* pointer;
+};
+
+/**
+ * @brief Typedef for PcmsInterpolatorHandle struct
+ * @copydetails PcmsInterpolatorHandle
+ * @see PcmsInterpolatorHandle, InterpolationBase, MLSMeshInterpolation,
+ * MLSPointCloudInterpolation
+ */
+typedef struct PcmsInterpolatorHandle PcmsInterpolatorHandle;
+
+/**
+ * @brief Create centroid to node interpolator
+ * @param oh_mesh Omega_h mesh handle
+ * @param radius Starting radius for support search
+ * @return PcmsInterpolatorHandle Handle to the created interpolator
+ *
+ * @details Creates an interpolator that maps values from element centroids to
+ * mesh nodes It uses MLS interpolation and uses the radius to start adaptive
+ * support search.
+ *
+ * @note Uses default parameters of MLSMeshInterpolation.
+ * @see MLSMeshInterpolation
+ */
+PcmsInterpolatorHandle pcms_create_interpolator(PcmsOmegaHMeshHandle oh_mesh,
+                                                double radius);
+
+/**
+ * @brief Create 2D point-based MLS interpolator for RBF interpolation
+ * @param source_points Void pointer to the first source point of double of size
+ * source_points_size
+ * @param source_points_size Size of source points array \f$
+ * Number\ of\ source\ points \times 2 \f$
+ * @param target_points Void pointer to the first target point of double of size
+ * target_points_size
+ * @param target_points_size Size of target points array \f$
+ * Number\ of\ target\ points \times 2 \f$
+ * @param radius Starting radius for support search (always adaptive)
+ * @param degree Degree of the MLS basis functions
+ * @param min_req_supports Minimum required supports for each target point
+ * (maximum allowed supports is three times this value)
+ * @param lambda Regularization parameter
+ * @param decay_factor Decay factor for weight function
+ * @return PcmsInterpolatorHandle Handle to the created interpolator
+ *
+ * @details To better select the interpolation parameters, look at their
+ * explanation in mls_interpolation function documentation. Call
+ * pcms_interpolate function with this interpolator to perform interpolation
+ * from source points to target points. Remember to delete the created
+ * interpolator using pcms_destroy_interpolator function after use to avoid
+ * memory leaks. The user is responsible for ensuring that the source and target
+ * points are sustained in memory during the call of this constructor. They are
+ * no longer needed after the interpolator is created as they are copied
+ * internally.
+ *
+ * @see mls_interpolation, InterpolationBase, MLSPointCloudInterpolation
+ */
+PcmsInterpolatorHandle pcms_create_point_based_interpolator(
+  void* source_points, int source_points_size, void* target_points,
+  int target_points_size, double radius, int degree, int min_req_supports,
+  double lambda, double decay_factor);
+
+/**
+ * @brief Create 2D point-based interpolator to map from Degas2 mesh element
+ * centroids (source) to XGC nodes (target)
+ * @param target_points Void pointer to the first target point of double of size
+ * target_points_size
+ * @param target_points_size Size of target points array \f$
+ * Number\ of\ target\ points \times 2 \f$
+ * @param dg2_mesh_filename C-string of Degas2 mesh filename (Omega_h (.osh)
+ * format)
+ * @param radius Starting radius for support search (always adaptive)
+ * @param dg2_elem_count Void pointer to integer to write the number of Degas2
+ * elements
+ * @param degree Degree of the MLS basis functions
+ * @param min_req_supports Minimum required supports for each target point
+ * (maximum allowed supports is three times this value)
+ * @param lambda Regularization parameter of MLS
+ * @param decay_factor Decay factor for weight function of MLS
+ * @return PcmsInterpolatorHandle Handle to the created interpolator
+ *
+ * @note This is used to interpolate after the Degas2 step in a coupled
+ * Degas2-XGC simulation.
+ * @copydetails pcms_create_point_based_interpolator
+ * @see mls_interpolation
+ */
+PcmsInterpolatorHandle pcms_create_degas2xgcnode_interpolator(
+  void* target_points, int target_points_size, const char* dg2_mesh_filename,
+  double radius, void* dg2_elem_count, int degree, int min_req_supports,
+  double lambda, double decay_factor);
+
+/**
+ * @brief Create 2D point-based interpolator to map from XGC nodes (source) to
+ * Degas2 mesh element centroids (target)
+ * @param dg2_mesh_filename C-string of Degas2 mesh filename (Omega_h (.osh)
+ * format)
+ * @param source_points Void pointer to the first source point of double of size
+ * source_points_size
+ * @param source_points_size Size of source points array \f$
+ * Number\ of\ source\ points \times 2 \f$
+ * @param radius Starting radius for support search (always adaptive)
+ * @param dg2_elem_count Void pointer to integer to write the number of Degas2
+ * elements
+ * @param degree Degree of the MLS basis functions
+ * @param min_req_supports Minimum required supports for each target point
+ * (maximum allowed supports is three times this value)
+ * @param lambda Regularization parameter of MLS
+ * @param decay_factor Decay factor for weight function of MLS
+ * @return PcmsInterpolatorHandle Handle to the created interpolator
+ *
+ * @note This is used to interpolate before the Degas2 step in a coupled
+ * Degas2-XGC simulation.
+ * @copydetails pcms_create_point_based_interpolator
+ * @see mls_interpolation
+ */
+PcmsInterpolatorHandle pcms_create_xgcnodedegas2_interpolator(
+  const char* dg2_mesh_filename, void* source_points, int source_points_size,
+  double radius, void* dg2_elem_count, int degree, int min_req_supports,
+  double lambda, double decay_factor);
+
+/**
+ * @brief Destroy interpolator
+ * @param interpolator Handle to the created interpolator
+ *
+ * @details Call this function to delete the created interpolator using
+ * pcms_create_interpolator or pcms_create_point_based_interpolator functions
+ * after use to avoid memory leaks.
+ */
+void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator);
+
+/**
+ * @brief Perform interpolation
+ * @param interpolator Handle to the created interpolator
+ * @param input Void pointer to the first element of the input data array of
+ * type double and size of input_size
+ * @param input_size Size of the input data array which is the same as the
+ * number of source points/elements
+ * @param output Void pointer to the first element of the output data array of
+ * type double and size of output_size
+ * @param output_size Size of the output data array which is the same as the
+ * number of target points/elements
+ *
+ * @details Uses the given interpolator to perform interpolation using the input
+ * data and writes the results to the output data array.
+ *
+ * @see pcms_create_interpolator, pcms_create_point_based_interpolator,
+ * mls_interpolation, InterpolationBase
+ */
+void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input,
+                      int input_size, void* output, int output_size);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // PCMS_INTERPOLATOR_CAPI_H

src/pcms/capi/mesh.cpp

@@ -0,0 +1,45 @@
+#include "mesh.h"
+
+#include <Omega_h_file.hpp>
+#include <Omega_h_library.hpp>
+#include <Omega_h_mesh.hpp>
+
+PcmsOmegaHLibraryHandle pcms_create_omega_h_library()
+{
+#ifdef Omega_H_USE_MPI
+  Omega_h::Library* oh_lib =
+    new Omega_h::Library(nullptr, nullptr, MPI_COMM_SELF);
+#else
+  Omega_h::Library* oh_lib = new Omega_h::Library(nullptr, nullptr);
+#endif
+
+  return {reinterpret_cast<void*>(oh_lib)};
+}
+
+PcmsOmegaHMeshHandle pcms_create_omega_h_mesh(
+  const char* filename, const PcmsOmegaHLibraryHandle oh_lib_handle)
+{
+  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);
+  Omega_h::Library* oh_lib =
+    reinterpret_cast<Omega_h::Library*>(oh_lib_handle.lib_handle);
+  auto* mesh = new Omega_h::Mesh(Omega_h::binary::read(fname, oh_lib->world()));
+
+  return {reinterpret_cast<void*>(mesh)};
+}
+
+void pcms_destroy_omega_h_mesh(PcmsOmegaHMeshHandle oh_mesh)
+{
+  if (oh_mesh.mesh_handle != nullptr) {
+    delete reinterpret_cast<Omega_h::Mesh*>(oh_mesh.mesh_handle);
+  }
+}
+
+void pcms_destroy_omega_h_library(PcmsOmegaHLibraryHandle oh_lib_handle)
+{
+  if (oh_lib_handle.lib_handle != nullptr) {
+    delete reinterpret_cast<Omega_h::Library*>(oh_lib_handle.lib_handle);
+  }
+}