From 7c3f0c56661bcc08b2598b3e265fad598c3332eb Mon Sep 17 00:00:00 2001 From: Fuad Hasan Date: Sat, 25 Jul 2026 15:03:33 -0400 Subject: [PATCH] Add conservative projection functionality to interpolator - runs but results are incorrect - followed the python example - the petsc mpi comm changes are temp fix for some error --- src/pcms/capi/interpolator.cpp | 154 +++++++++++ src/pcms/capi/interpolator.h | 77 ++++++ src/pcms/fortranapi/interpolator.i | 22 ++ src/pcms/fortranapi/interpolator_wrap.c | 165 ++++++++++++ src/pcms/fortranapi/pcms_interpolator.f90 | 249 ++++++++++++++++++ .../conservative_projection_solver.cpp | 26 +- .../omega_h_intersection_rhs_integrator.cpp | 10 +- src/pcms/transfer/omega_h_mass_integrator.cpp | 16 +- .../transfer/omega_h_mc_rhs_integrator.cpp | 10 +- 9 files changed, 698 insertions(+), 31 deletions(-) diff --git a/src/pcms/capi/interpolator.cpp b/src/pcms/capi/interpolator.cpp index 9b5703d17..05193ce45 100644 --- a/src/pcms/capi/interpolator.cpp +++ b/src/pcms/capi/interpolator.cpp @@ -133,3 +133,157 @@ void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, mls_interpolator->eval(input_array, output_array); } + +// --------------------------------------------------------------------------- +// Conservative Projection (mesh-intersection-based) +// --------------------------------------------------------------------------- + +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + +// Trims trailing whitespace from a C string (Fortran strings may be padded). +std::string trim_filename(const char* filename) +{ + auto fname = std::string(filename); + fname.erase(fname.find_last_not_of(" \n\r\t") + 1); + return fname; +} + +struct ConservativeProjectionContext +{ + // Library must be declared before meshes (Omega_h::Mesh holds a pointer to it). + Omega_h::Library library; + Omega_h::Mesh source_mesh; + Omega_h::Mesh target_mesh; + pcms::LagrangeFunctionSpace source_space; + pcms::LagrangeFunctionSpace target_space; + std::unique_ptr projection; + pcms::Field source_field; + pcms::Field target_field; + + ConservativeProjectionContext(const char* source_mesh_name, + const char* target_mesh_name, int src_order, + int tgt_order) + : library(nullptr, nullptr, MPI_COMM_SELF), + source_mesh(Omega_h::binary::read(trim_filename(source_mesh_name), + library.world())), + target_mesh(Omega_h::binary::read(trim_filename(target_mesh_name), + library.world())), + source_space(pcms::LagrangeFunctionSpace::FromMesh( + source_mesh, src_order, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH)), + target_space(pcms::LagrangeFunctionSpace::FromMesh( + target_mesh, tgt_order, 1, pcms::CoordinateSystem::Cartesian, "global", + pcms::LagrangeFunctionSpace::Backend::OmegaH)), + projection(std::make_unique( + source_space, target_space)), + source_field(source_space.CreateField()), + target_field(target_space.CreateField()) + { + } +}; + +} // namespace + +PcmsConservativeProjectionHandle pcms_create_conservative_projection( + const char* source_mesh_name, int source_order, + const char* target_mesh_name, int target_order) +{ + auto* ctx = new ConservativeProjectionContext(source_mesh_name, + target_mesh_name, source_order, + target_order); + return {reinterpret_cast(ctx)}; +} + +int pcms_conservative_projection_get_source_size( + PcmsConservativeProjectionHandle projection) +{ + auto* ctx = + reinterpret_cast(projection.pointer); + return ctx->source_space.GetLayout()->GetNumOwnedDofHolder(); +} + +int pcms_conservative_projection_get_target_size( + PcmsConservativeProjectionHandle projection) +{ + auto* ctx = + reinterpret_cast(projection.pointer); + return ctx->target_space.GetLayout()->GetNumOwnedDofHolder(); +} + +void pcms_conservative_projection_apply( + PcmsConservativeProjectionHandle projection, void* source_data, + int source_size, void* target_data, int target_size) +{ + auto* ctx = + reinterpret_cast(projection.pointer); + + // Copy source data into internal field + auto source_view = pcms::Rank2View( + reinterpret_cast(source_data), source_size, 1); + ctx->source_field.SetDOFHolderDataHost(source_view); + + // Apply conservative projection + ctx->projection->Apply(ctx->source_field, ctx->target_field); + + // Copy result back to user buffer + auto target_view = ctx->target_field.GetDOFHolderDataHost(); + auto flat = pcms::FlattenToRank1View(target_view); + std::memcpy(target_data, flat.data_handle(), + static_cast(target_size) * sizeof(double)); +} + +void pcms_destroy_conservative_projection( + PcmsConservativeProjectionHandle projection) +{ + if (projection.pointer != nullptr) { + delete reinterpret_cast( + projection.pointer); + } +} + +#else // !(PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS) + +// Stub implementations when PETSc or MeshFields are unavailable + +PcmsConservativeProjectionHandle pcms_create_conservative_projection( + const char*, int, const char*, int) +{ + pcms::printError("Conservative projection requires PCMS_ENABLE_PETSC and " + "PCMS_ENABLE_MESHFIELDS\n"); + return {nullptr}; +} + +int pcms_conservative_projection_get_source_size( + PcmsConservativeProjectionHandle) +{ + return 0; +} + +int pcms_conservative_projection_get_target_size( + PcmsConservativeProjectionHandle) +{ + return 0; +} + +void pcms_conservative_projection_apply(PcmsConservativeProjectionHandle, + void*, int, void*, int) +{ + pcms::printError("Conservative projection requires PCMS_ENABLE_PETSC and " + "PCMS_ENABLE_MESHFIELDS\n"); +} + +void pcms_destroy_conservative_projection(PcmsConservativeProjectionHandle) {} + +#endif diff --git a/src/pcms/capi/interpolator.h b/src/pcms/capi/interpolator.h index c557e1a29..bf67713cf 100644 --- a/src/pcms/capi/interpolator.h +++ b/src/pcms/capi/interpolator.h @@ -175,6 +175,83 @@ void pcms_destroy_interpolator(PcmsInterpolatorHandle interpolator); void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size); +/** + * @brief Holds a void pointer to ConservativeProjectionContext + */ +struct PcmsConservativeProjectionHandle +{ + void* pointer; +}; + +/** + * @brief Typedef for PcmsConservativeProjectionHandle struct + * @copydetails PcmsConservativeProjectionHandle + */ +typedef struct PcmsConservativeProjectionHandle PcmsConservativeProjectionHandle; + +/** + * @brief Create a mesh-intersection-based conservative projection + * @param source_mesh_name C-string path to source Omega_h mesh file (.osh) + * @param source_order Order of source Lagrange space (0 or 1) + * @param target_mesh_name C-string path to target Omega_h mesh file (.osh) + * @param target_order Order of target Lagrange space (0 or 1) + * @return Handle to the created conservative projection + * + * @details Loads both meshes internally, creates LagrangeFunctionSpace objects + * (Backend::OmegaH), performs mesh intersection via intersectTargets(), builds + * quadrature data, assembles and factors the target mass matrix (PETSc KSP), + * and caches everything for repeated Apply calls. The meshes are owned by + * the returned handle and live until pcms_destroy_conservative_projection. + * + * @note Requires PCMS_ENABLE_PETSC and PCMS_ENABLE_MESHFIELDS. Returns a + * handle with null pointer if either is unavailable. + * @note Only scalar fields (1 component) are supported. + * @note Source and target orders can differ (e.g., P0<->P1). + */ +PcmsConservativeProjectionHandle pcms_create_conservative_projection( + const char* source_mesh_name, int source_order, + const char* target_mesh_name, int target_order); + +/** + * @brief Get the number of source DOF holders + * @param projection Handle to the conservative projection + * @return Number of source DOF holders + */ +int pcms_conservative_projection_get_source_size( + PcmsConservativeProjectionHandle projection); + +/** + * @brief Get the number of target DOF holders + * @param projection Handle to the conservative projection + * @return Number of target DOF holders + */ +int pcms_conservative_projection_get_target_size( + PcmsConservativeProjectionHandle projection); + +/** + * @brief Apply conservative projection + * @param projection Handle to the projection + * @param source_data Flat array of source field values (size = source_size) + * @param source_size Number of source DOF holders + * @param target_data Flat array to receive target field values (size = target_size) + * @param target_size Number of target DOF holders + * + * @details Copies source data into the internal field, calls + * OmegaHConservativeProjection::Apply(), and copies the result back. + * This is the cheap per-time-step call: no mesh intersection, + * no matrix assembly, no refactorization. + */ +void pcms_conservative_projection_apply( + PcmsConservativeProjectionHandle projection, void* source_data, + int source_size, void* target_data, int target_size); + +/** + * @brief Destroy conservative projection and free all internal resources + * @param projection Handle to the conservative projection to destroy + */ +void pcms_destroy_conservative_projection( + PcmsConservativeProjectionHandle projection); + #ifdef __cplusplus } #endif diff --git a/src/pcms/fortranapi/interpolator.i b/src/pcms/fortranapi/interpolator.i index 970f3eec0..626dbdf84 100644 --- a/src/pcms/fortranapi/interpolator.i +++ b/src/pcms/fortranapi/interpolator.i @@ -29,3 +29,25 @@ void pcms_kokkos_initialize_without_args(); void pcms_kokkos_finalize(); void pcms_interpolate(PcmsInterpolatorHandle interpolator, void* input, int input_size, void* output, int output_size); + + +// --- Conservative Projection (mesh-intersection-based) --- + +struct PcmsConservativeProjectionHandle { + void* pointer; +}; +typedef struct PcmsConservativeProjectionHandle PcmsConservativeProjectionHandle; + +PcmsConservativeProjectionHandle pcms_create_conservative_projection( + const char* source_mesh_name, int source_order, + const char* target_mesh_name, int target_order); +int pcms_conservative_projection_get_source_size( + PcmsConservativeProjectionHandle projection); +int pcms_conservative_projection_get_target_size( + PcmsConservativeProjectionHandle projection); +void pcms_conservative_projection_apply( + PcmsConservativeProjectionHandle projection, + void* source_data, int source_size, + void* target_data, int target_size); +void pcms_destroy_conservative_projection( + PcmsConservativeProjectionHandle projection); diff --git a/src/pcms/fortranapi/interpolator_wrap.c b/src/pcms/fortranapi/interpolator_wrap.c index d0d7e4ccb..012eb02ac 100644 --- a/src/pcms/fortranapi/interpolator_wrap.c +++ b/src/pcms/fortranapi/interpolator_wrap.c @@ -531,3 +531,168 @@ SWIGEXPORT void _wrap_pcms_interpolate(SwigClassWrapper* farg1, arg5 = (int)(*farg5); pcms_interpolate(arg1, arg2, arg3, arg4, arg5); } + +SWIGEXPORT void _wrap_PcmsConservativeProjectionHandle_pointer_set( + SwigClassWrapper* farg1, void const** farg2) +{ + struct PcmsConservativeProjectionHandle* arg1 = + (struct PcmsConservativeProjectionHandle*)0; + void* arg2 = (void*)0; + + SWIG_check_nonnull(farg1->cptr, "struct PcmsConservativeProjectionHandle *", + "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle::pointer", return); + arg1 = (struct PcmsConservativeProjectionHandle*)farg1->cptr; + arg2 = (void*)(*farg2); + if (arg1) + (arg1)->pointer = arg2; +} + +SWIGEXPORT void* _wrap_PcmsConservativeProjectionHandle_pointer_get( + SwigClassWrapper* farg1) +{ + void* fresult; + struct PcmsConservativeProjectionHandle* arg1 = + (struct PcmsConservativeProjectionHandle*)0; + void* result = 0; + + SWIG_check_nonnull(farg1->cptr, "struct PcmsConservativeProjectionHandle *", + "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle::pointer", return 0); + arg1 = (struct PcmsConservativeProjectionHandle*)farg1->cptr; + result = (void*)((arg1)->pointer); + fresult = (void*)(result); + return fresult; +} + +SWIGEXPORT SwigClassWrapper _wrap_new_PcmsConservativeProjectionHandle() +{ + SwigClassWrapper fresult; + struct PcmsConservativeProjectionHandle* result = 0; + + result = (struct PcmsConservativeProjectionHandle*)calloc( + 1, sizeof(struct PcmsConservativeProjectionHandle)); + fresult.cptr = (void*)result; + fresult.cmemflags = SWIG_MEM_RVALUE | (1 ? SWIG_MEM_OWN : 0); + return fresult; +} + +SWIGEXPORT void _wrap_delete_PcmsConservativeProjectionHandle( + SwigClassWrapper* farg1) +{ + struct PcmsConservativeProjectionHandle* arg1 = + (struct PcmsConservativeProjectionHandle*)0; + + arg1 = (struct PcmsConservativeProjectionHandle*)farg1->cptr; + free((char*)arg1); +} + +SWIGEXPORT void _wrap_PcmsConservativeProjectionHandle_op_assign__( + SwigClassWrapper* farg1, SwigClassWrapper* farg2) +{ + struct PcmsConservativeProjectionHandle* arg1 = + (struct PcmsConservativeProjectionHandle*)0; + struct PcmsConservativeProjectionHandle* arg2 = 0; + + (void)sizeof(arg1); + (void)sizeof(arg2); + SWIG_assign(farg1, *farg2); +} + +SWIGEXPORT SwigClassWrapper _wrap_pcms_create_conservative_projection( + SwigArrayWrapper* farg1, int const* farg2, SwigArrayWrapper* farg3, + int const* farg4) +{ + SwigClassWrapper fresult; + char* arg1 = (char*)0; + int arg2; + char* arg3 = (char*)0; + int arg4; + PcmsConservativeProjectionHandle result; + + arg1 = (char*)(farg1->data); + arg2 = (int)(*farg2); + arg3 = (char*)(farg3->data); + arg4 = (int)(*farg4); + result = pcms_create_conservative_projection((char const*)arg1, arg2, + (char const*)arg3, arg4); + fresult.cptr = (PcmsConservativeProjectionHandle*)memcpy( + (PcmsConservativeProjectionHandle*)calloc( + 1, sizeof(PcmsConservativeProjectionHandle)), + &result, sizeof(PcmsConservativeProjectionHandle)); + fresult.cmemflags = SWIG_MEM_RVALUE | SWIG_MEM_OWN; + return fresult; +} + +SWIGEXPORT int _wrap_pcms_conservative_projection_get_source_size( + SwigClassWrapper* farg1) +{ + int fresult; + PcmsConservativeProjectionHandle arg1; + int result; + + SWIG_check_nonnull(farg1->cptr, "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle", + "pcms_conservative_projection_get_source_size(" + "PcmsConservativeProjectionHandle)", + return 0); + arg1 = *((PcmsConservativeProjectionHandle*)(farg1->cptr)); + result = (int)pcms_conservative_projection_get_source_size(arg1); + fresult = (int)(result); + return fresult; +} + +SWIGEXPORT int _wrap_pcms_conservative_projection_get_target_size( + SwigClassWrapper* farg1) +{ + int fresult; + PcmsConservativeProjectionHandle arg1; + int result; + + SWIG_check_nonnull(farg1->cptr, "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle", + "pcms_conservative_projection_get_target_size(" + "PcmsConservativeProjectionHandle)", + return 0); + arg1 = *((PcmsConservativeProjectionHandle*)(farg1->cptr)); + result = (int)pcms_conservative_projection_get_target_size(arg1); + fresult = (int)(result); + return fresult; +} + +SWIGEXPORT void _wrap_pcms_conservative_projection_apply( + SwigClassWrapper* farg1, void const** farg2, int const* farg3, + void const** farg4, int const* farg5) +{ + PcmsConservativeProjectionHandle arg1; + void* arg2 = (void*)0; + int arg3; + void* arg4 = (void*)0; + int arg5; + + SWIG_check_nonnull(farg1->cptr, "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle", + "pcms_conservative_projection_apply(" + "PcmsConservativeProjectionHandle,void *,int,void *,int)", + return); + arg1 = *((PcmsConservativeProjectionHandle*)(farg1->cptr)); + arg2 = (void*)(*farg2); + arg3 = (int)(*farg3); + arg4 = (void*)(*farg4); + arg5 = (int)(*farg5); + pcms_conservative_projection_apply(arg1, arg2, arg3, arg4, arg5); +} + +SWIGEXPORT void _wrap_pcms_destroy_conservative_projection( + SwigClassWrapper* farg1) +{ + PcmsConservativeProjectionHandle arg1; + + SWIG_check_nonnull( + farg1->cptr, "PcmsConservativeProjectionHandle", + "PcmsConservativeProjectionHandle", + "pcms_destroy_conservative_projection(PcmsConservativeProjectionHandle)", + return); + arg1 = *((PcmsConservativeProjectionHandle*)(farg1->cptr)); + pcms_destroy_conservative_projection(arg1); +} diff --git a/src/pcms/fortranapi/pcms_interpolator.f90 b/src/pcms/fortranapi/pcms_interpolator.f90 index 2e8bdaf56..a09b56a8c 100644 --- a/src/pcms/fortranapi/pcms_interpolator.f90 +++ b/src/pcms/fortranapi/pcms_interpolator.f90 @@ -39,6 +39,24 @@ module pcms_interpolator public :: pcms_kokkos_initialize_without_args public :: pcms_kokkos_finalize public :: pcms_interpolate + ! struct struct PcmsConservativeProjectionHandle + type, public :: PcmsConservativeProjectionHandle + type(SwigClassWrapper), public :: swigdata + contains + procedure :: set_pointer => swigf_PcmsConservativeProjectionHandle_pointer_set + procedure :: get_pointer => swigf_PcmsConservativeProjectionHandle_pointer_get + procedure :: release => swigf_PcmsConservativeProjectionHandle_release + procedure, private :: swigf_PcmsConservativeProjectionHandle_op_assign__ + generic :: assignment(=) => swigf_PcmsConservativeProjectionHandle_op_assign__ + end type PcmsConservativeProjectionHandle + public :: pcms_create_conservative_projection + public :: pcms_conservative_projection_get_source_size + public :: pcms_conservative_projection_get_target_size + public :: pcms_conservative_projection_apply + public :: pcms_destroy_conservative_projection + interface PcmsConservativeProjectionHandle + module procedure swigf_new_PcmsConservativeProjectionHandle + end interface interface PcmsInterpolatorHandle module procedure swigf_new_PcmsInterpolatorHandle end interface @@ -176,6 +194,95 @@ subroutine swigc_pcms_interpolate(farg1, farg2, farg3, farg4, farg5) & integer(C_INT), intent(in) :: farg5 end subroutine +subroutine swigc_PcmsConservativeProjectionHandle_pointer_set(farg1, farg2) & +bind(C, name="_wrap_PcmsConservativeProjectionHandle_pointer_set") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +type(C_PTR), intent(in) :: farg2 +end subroutine + +function swigc_PcmsConservativeProjectionHandle_pointer_get(farg1) & +bind(C, name="_wrap_PcmsConservativeProjectionHandle_pointer_get") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +type(C_PTR) :: fresult +end function + +function swigc_new_PcmsConservativeProjectionHandle() & +bind(C, name="_wrap_new_PcmsConservativeProjectionHandle") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper) :: fresult +end function + +subroutine swigc_delete_PcmsConservativeProjectionHandle(farg1) & +bind(C, name="_wrap_delete_PcmsConservativeProjectionHandle") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(inout) :: farg1 +end subroutine + +subroutine swigc_PcmsConservativeProjectionHandle_op_assign__(farg1, farg2) & +bind(C, name="_wrap_PcmsConservativeProjectionHandle_op_assign__") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(inout) :: farg1 +type(SwigClassWrapper), intent(in) :: farg2 +end subroutine + +function swigc_pcms_create_conservative_projection(farg1, farg2, farg3, farg4) & +bind(C, name="_wrap_pcms_create_conservative_projection") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigarraywrapper +import :: swigclasswrapper +type(SwigArrayWrapper) :: farg1 +integer(C_INT), intent(in) :: farg2 +type(SwigArrayWrapper) :: farg3 +integer(C_INT), intent(in) :: farg4 +type(SwigClassWrapper) :: fresult +end function + +function swigc_pcms_conservative_projection_get_source_size(farg1) & +bind(C, name="_wrap_pcms_conservative_projection_get_source_size") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +integer(C_INT) :: fresult +end function + +function swigc_pcms_conservative_projection_get_target_size(farg1) & +bind(C, name="_wrap_pcms_conservative_projection_get_target_size") & +result(fresult) +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +integer(C_INT) :: fresult +end function + +subroutine swigc_pcms_conservative_projection_apply(farg1, farg2, farg3, farg4, farg5) & +bind(C, name="_wrap_pcms_conservative_projection_apply") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +type(C_PTR), intent(in) :: farg2 +integer(C_INT), intent(in) :: farg3 +type(C_PTR), intent(in) :: farg4 +integer(C_INT), intent(in) :: farg5 +end subroutine + +subroutine swigc_pcms_destroy_conservative_projection(farg1) & +bind(C, name="_wrap_pcms_destroy_conservative_projection") +use, intrinsic :: ISO_C_BINDING +import :: swigclasswrapper +type(SwigClassWrapper), intent(in) :: farg1 +end subroutine + end interface @@ -435,5 +542,147 @@ subroutine pcms_interpolate(interpolator, input, input_size, output, output_size call swigc_pcms_interpolate(farg1, farg2, farg3, farg4, farg5) end subroutine +subroutine swigf_PcmsConservativeProjectionHandle_pointer_set(self, pointer) +use, intrinsic :: ISO_C_BINDING +class(PcmsConservativeProjectionHandle), intent(in) :: self +type(C_PTR), intent(in) :: pointer +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 + +farg1 = self%swigdata +farg2 = pointer +call swigc_PcmsConservativeProjectionHandle_pointer_set(farg1, farg2) +end subroutine + +function swigf_PcmsConservativeProjectionHandle_pointer_get(self) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(C_PTR) :: swig_result +class(PcmsConservativeProjectionHandle), intent(in) :: self +type(C_PTR) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +fresult = swigc_PcmsConservativeProjectionHandle_pointer_get(farg1) +swig_result = fresult +end function + +function swigf_new_PcmsConservativeProjectionHandle() & +result(self) +use, intrinsic :: ISO_C_BINDING +type(PcmsConservativeProjectionHandle) :: self +type(SwigClassWrapper) :: fresult + +fresult = swigc_new_PcmsConservativeProjectionHandle() +self%swigdata = fresult +end function + +subroutine swigf_PcmsConservativeProjectionHandle_release(self) +use, intrinsic :: ISO_C_BINDING +class(PcmsConservativeProjectionHandle), intent(inout) :: self +type(SwigClassWrapper) :: farg1 + +farg1 = self%swigdata +if (btest(farg1%cmemflags, swig_cmem_own_bit)) then +call swigc_delete_PcmsConservativeProjectionHandle(farg1) +endif +farg1%cptr = C_NULL_PTR +farg1%cmemflags = 0 +self%swigdata = farg1 +end subroutine + +subroutine swigf_PcmsConservativeProjectionHandle_op_assign__(self, other) +use, intrinsic :: ISO_C_BINDING +class(PcmsConservativeProjectionHandle), intent(inout) :: self +type(PcmsConservativeProjectionHandle), intent(in) :: other +type(SwigClassWrapper) :: farg1 +type(SwigClassWrapper) :: farg2 + +farg1 = self%swigdata +farg2 = other%swigdata +call swigc_PcmsConservativeProjectionHandle_op_assign__(farg1, farg2) +self%swigdata = farg1 +end subroutine + +function pcms_create_conservative_projection(source_mesh_name, source_order, target_mesh_name, target_order) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +type(PcmsConservativeProjectionHandle) :: swig_result +character(len=*), intent(in) :: source_mesh_name +integer(C_INT), intent(in) :: source_order +character(len=*), intent(in) :: target_mesh_name +integer(C_INT), intent(in) :: target_order +type(SwigClassWrapper) :: fresult +character(kind=C_CHAR), dimension(:), allocatable, target :: farg1_temp +type(SwigArrayWrapper) :: farg1 +integer(C_INT) :: farg2 +character(kind=C_CHAR), dimension(:), allocatable, target :: farg3_temp +type(SwigArrayWrapper) :: farg3 +integer(C_INT) :: farg4 + +call SWIGTM_fin_char_Sm_(source_mesh_name, farg1, farg1_temp) +farg2 = source_order +call SWIGTM_fin_char_Sm_(target_mesh_name, farg3, farg3_temp) +farg4 = target_order +fresult = swigc_pcms_create_conservative_projection(farg1, farg2, farg3, farg4) +swig_result%swigdata = fresult +end function + +function pcms_conservative_projection_get_source_size(projection) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(PcmsConservativeProjectionHandle), intent(in) :: projection +integer(C_INT) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = projection%swigdata +fresult = swigc_pcms_conservative_projection_get_source_size(farg1) +swig_result = fresult +end function + +function pcms_conservative_projection_get_target_size(projection) & +result(swig_result) +use, intrinsic :: ISO_C_BINDING +integer(C_INT) :: swig_result +type(PcmsConservativeProjectionHandle), intent(in) :: projection +integer(C_INT) :: fresult +type(SwigClassWrapper) :: farg1 + +farg1 = projection%swigdata +fresult = swigc_pcms_conservative_projection_get_target_size(farg1) +swig_result = fresult +end function + +subroutine pcms_conservative_projection_apply(projection, source_data, source_size, target_data, target_size) +use, intrinsic :: ISO_C_BINDING +type(PcmsConservativeProjectionHandle), intent(in) :: projection +type(C_PTR), intent(in) :: source_data +integer(C_INT), intent(in) :: source_size +type(C_PTR), intent(in) :: target_data +integer(C_INT), intent(in) :: target_size +type(SwigClassWrapper) :: farg1 +type(C_PTR) :: farg2 +integer(C_INT) :: farg3 +type(C_PTR) :: farg4 +integer(C_INT) :: farg5 + +farg1 = projection%swigdata +farg2 = source_data +farg3 = source_size +farg4 = target_data +farg5 = target_size +call swigc_pcms_conservative_projection_apply(farg1, farg2, farg3, farg4, farg5) +end subroutine + +subroutine pcms_destroy_conservative_projection(projection) +use, intrinsic :: ISO_C_BINDING +type(PcmsConservativeProjectionHandle), intent(in) :: projection +type(SwigClassWrapper) :: farg1 + +farg1 = projection%swigdata +call swigc_pcms_destroy_conservative_projection(farg1) +end subroutine + end module diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index e6955580a..e58486b2b 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -12,11 +12,11 @@ static Omega_h::Reals vecToOmegaHReals(Vec vec) { PetscInt n = 0; PetscErrorCode ierr = VecGetSize(vec, &n); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); const PetscScalar* array = nullptr; ierr = VecGetArrayRead(vec, &array); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto values_host = Omega_h::HostWrite(n); for (PetscInt i = 0; i < n; ++i) { @@ -24,7 +24,7 @@ static Omega_h::Reals vecToOmegaHReals(Vec vec) } ierr = VecRestoreArrayRead(vec, &array); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); return Omega_h::Reals(values_host); } @@ -42,7 +42,7 @@ GalerkinProjectionSolver::GalerkinProjectionSolver( PetscInt m = 0, n = 0; PetscErrorCode ierr = MatGetSize(A, &m, &n); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); nverts_ = m; const std::size_t num_pts = @@ -50,14 +50,14 @@ GalerkinProjectionSolver::GalerkinProjectionSolver( sampled_values_ = Kokkos::View("rhs_sampled", num_pts, 1); - ierr = KSPCreate(PETSC_COMM_WORLD, &ksp_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + ierr = KSPCreate(PETSC_COMM_SELF, &ksp_); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSetOperators(ksp_, A, A); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSetFromOptions(ksp_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSetUp(ksp_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } GalerkinProjectionSolver::~GalerkinProjectionSolver() @@ -82,16 +82,16 @@ Omega_h::Reals GalerkinProjectionSolver::Solve( Vec rhs_vector = rhs_integrator_->GetVector(); Vec solution = nullptr; - PetscErrorCode ierr = createSeqVec(PETSC_COMM_WORLD, nverts_, &solution); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + PetscErrorCode ierr = createSeqVec(PETSC_COMM_SELF, nverts_, &solution); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = KSPSolve(ksp_, rhs_vector, solution); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto result = vecToOmegaHReals(solution); ierr = VecDestroy(&solution); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); return result; } diff --git a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp index ae9b293f2..b5baad53c 100644 --- a/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_intersection_rhs_integrator.cpp @@ -194,15 +194,15 @@ OmegaHIntersectionRHSIntegrator::OmegaHIntersectionRHSIntegrator( const PetscInt nnz = static_cast(node_gids_.extent(0)); PetscErrorCode ierr = - createSeqVec(PETSC_COMM_WORLD, data.num_target_dofs, &vec_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + createSeqVec(PETSC_COMM_SELF, data.num_target_dofs, &vec_); + CHKERRABORT(PETSC_COMM_SELF, ierr); // VecSetPreallocationCOO takes the COO indices on the host // TODO: ask Todd/PETSc folks if there is a better way to do this for GPU // support auto node_gids_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, node_gids_); ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_host.data()); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } OmegaHIntersectionRHSIntegrator::~OmegaHIntersectionRHSIntegrator() @@ -239,7 +239,7 @@ void OmegaHIntersectionRHSIntegrator::Assemble( PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); PetscErrorCode ierr = VecZeroEntries(vec_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto sv = Kokkos::View>( @@ -258,7 +258,7 @@ void OmegaHIntersectionRHSIntegrator::Assemble( }); ierr = VecSetValuesCOO(vec_, coo_vals.data(), ADD_VALUES); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } // --------------------------------------------------------------------------- diff --git a/src/pcms/transfer/omega_h_mass_integrator.cpp b/src/pcms/transfer/omega_h_mass_integrator.cpp index 29cd90e2e..0775b741f 100644 --- a/src/pcms/transfer/omega_h_mass_integrator.cpp +++ b/src/pcms/transfer/omega_h_mass_integrator.cpp @@ -101,17 +101,17 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( coo_cols, vals); PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + createSeqAIJMat(PETSC_COMM_SELF, num_dofs, num_dofs, 0, nullptr, &mat_); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto coo_rows_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); auto coo_cols_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), coo_cols_host.data()); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = MatSetValuesCOO(mat_, vals.data(), INSERT_VALUES); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); return; } @@ -139,17 +139,17 @@ OmegaHMassIntegrator::OmegaHMassIntegrator( // elm_mass_dev is in the same element-major order as coo_rows/coo_cols, so // it can be passed directly to MatSetValuesCOO — no host copy needed. PetscErrorCode ierr = - createSeqAIJMat(PETSC_COMM_WORLD, num_dofs, num_dofs, 0, nullptr, &mat_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + createSeqAIJMat(PETSC_COMM_SELF, num_dofs, num_dofs, 0, nullptr, &mat_); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto coo_rows_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_rows); auto coo_cols_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coo_cols); ierr = MatSetPreallocationCOO(mat_, nnz, coo_rows_host.data(), coo_cols_host.data()); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); ierr = MatSetValuesCOO(mat_, elm_mass_dev.data(), INSERT_VALUES); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } OmegaHMassIntegrator::~OmegaHMassIntegrator() diff --git a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp index c57ce9f1f..1684ebbbc 100644 --- a/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp +++ b/src/pcms/transfer/omega_h_mc_rhs_integrator.cpp @@ -136,15 +136,15 @@ OmegaHMonteCarloRHSIntegrator::OmegaHMonteCarloRHSIntegrator( const PetscInt nnz = static_cast(node_gids_.extent(0)); PetscErrorCode ierr = - createSeqVec(PETSC_COMM_WORLD, static_cast(mesh.nverts()), &vec_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + createSeqVec(PETSC_COMM_SELF, static_cast(mesh.nverts()), &vec_); + CHKERRABORT(PETSC_COMM_SELF, ierr); // VecSetPreallocationCOO takes the COO indices on the host // TODO: ask Todd/PETSc folks if there is a better way to do this for GPU // support auto node_gids_host = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, node_gids_); ierr = VecSetPreallocationCOO(vec_, nnz, node_gids_host.data()); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } OmegaHMonteCarloRHSIntegrator::~OmegaHMonteCarloRHSIntegrator() @@ -180,7 +180,7 @@ void OmegaHMonteCarloRHSIntegrator::Assemble( PCMS_ALWAYS_ASSERT(sampled_values.extent(1) >= 1); PetscErrorCode ierr = VecZeroEntries(vec_); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); auto sv = Kokkos::View>( @@ -198,7 +198,7 @@ void OmegaHMonteCarloRHSIntegrator::Assemble( }); ierr = VecSetValuesCOO(vec_, coo_vals.data(), ADD_VALUES); - CHKERRABORT(PETSC_COMM_WORLD, ierr); + CHKERRABORT(PETSC_COMM_SELF, ierr); } } // namespace pcms