From 4bf033af26f3ece8d836b891deeb7f71646ed059 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 13:43:00 -0400 Subject: [PATCH 01/44] add mc conservative projection file --- src/pcms/transfer/pcms_mc_proj_utils.cpp | 257 +++++++++++++++++++++++ src/pcms/transfer/pcms_mc_proj_utils.hpp | 228 ++++++++++++++++++++ 2 files changed, 485 insertions(+) create mode 100644 src/pcms/transfer/pcms_mc_proj_utils.cpp create mode 100644 src/pcms/transfer/pcms_mc_proj_utils.hpp diff --git a/src/pcms/transfer/pcms_mc_proj_utils.cpp b/src/pcms/transfer/pcms_mc_proj_utils.cpp new file mode 100644 index 00000000..152ff9c9 --- /dev/null +++ b/src/pcms/transfer/pcms_mc_proj_utils.cpp @@ -0,0 +1,257 @@ +#include + +constexpr int MAX_POINTS = 1600; + +// TODO:: create a function that can sample sobol sequences instead of +// generating in python and reading and using +Kokkos::View read_sobol_barycentric_samples_from_file( + std::string file_path) +{ + std::ifstream infile(file_path); + if (!infile) { + throw std::runtime_error("Could not open sobol sample file : "; + } + + std::vector buffer; + size_t nrows = 0, ncols = 0; + std::string line; + + // skip header row + std::getline(infile, line); + + // Read the file and store values in a flat vector + while (std::getline(infile, line)) { + if (line.empty()) + continue; + + std::istringstream iss(line); + MeshField::Real b0, b1, b2; + + if (!(iss >> b0 >> b1 >> b2)) { + throw std::runtime_error("Invalid barycentric row in file: " + file_path); + } + + buffer.push_back(b0); + buffer.push_back(b1); + buffer.push_back(b2); + ++col; + ++nrows; + } + + Kokkos::View device_samples("device_data", nrows); + + auto host_samples = Kokkos::create_mirror_view(device_sobol_samples); + + // Fill host view from buffer + for (size_t i = 0; i < nrows; ++i) { + for (size_t j = 0; j < ncols; ++j) { + host_samples(i, j) = buffer[i * ncols + j]; + } + } + + Kokkos::deep_copy(device_samples, host_samples); + + return device_samples; +} + +Kokkos::View generate_uniform_random_barycentric_coords( + const int npoints_each_tri) +{ + + Kokkos::Random_XorShift64_Pool<> rand_pool(1200); + Kokkos::View barycentric_coords( + "uniform random barycentric coordinates", npoints_each_tri); + Omega_h::parallel_for( + npoints_each_tri, OMEGA_H_LAMBDA(int i) { + // taken from Shape Distributions (ACM Transactions on Graphics, Vol. + // 21, No. 4, October 2002.) page 814 Eq 1 + auto rand_gen = rand_pool.get_state(); + Omega_h::Real r1 = rand_gen.drand(); + Omega_h::Real r2 = rand_gen.drand(); + Omega_h::Real sqrt_r1 = Kokkos::sqrt(r1); + barycentric_coords(i, 0) = 1.0 - sqrt_r1; + barycentric_coords(i, 1) = sqrt_r1 * (1.0 - r2); + barycentric_coords(i, 2) = sqrt_r1 * r2; + rand_pool.free_state(rand_gen); + }); + + return barycentric_coords; +} + +Kokkos::View global_coords_from_ref_barycentric_coords( + Omega_h::Mesh& target_mesh, + const Kokkos::View& ref_barycentric_coords) +{ + + int nelems = target_mesh.nelems(); + + const auto& faces2nodes = + target_mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& coordinates = target_mesh.coords(); + const int npoints_each_tri = ref_barycentric_coords.extent(0); + + Kokkos::View global_coords( + "stores global coordinates of reference samples in each element", + nelems * npoints_each_tri); + + Omega_h::parallel_for( + nelems, OMEGA_H_LAMBDA(int id) { + const auto elem_verts = Omega_h::gather_verts<3>(faces2nodes, id); + const Omega_h::Few, 3> verts_coords = + Omega_h::gather_vectors<3, 2>(coordinates, elem_verts); + + int base_idx = id * npoints_each_tri; + for (int i = 0; i < npoints_each_tri; ++i) { + + global_coords(base_idx + i, 0) = + verts_coords[0][0] * ref_barycentric_coords(i, 0) + + verts_coords[1][0] * ref_barycentric_coords(i, 1) + + verts_coords[2][0] * ref_barycentric_coords(i, 2); + + global_coords(base_idx + i, 1) = + verts_coords[0][1] * ref_barycentric_coords(i, 0) + + verts_coords[1][1] * ref_barycentric_coords(i, 1) + + verts_coords[2][1] * ref_barycentric_coords(i, 2); + } + }); + + return global_coords; +} + +Kokkos::View localize_points_in_mesh( + Omega_h::Mesh& mesh, const Kokkos::View& points) +{ + + PCMS_ALWAYS_ASSERT(mesh.dim() == 2); + const auto mesh_bbox = Omega_h::get_bounding_box<2>(&mesh); + + const auto npoints = points.extent(0); + + Kokkos::Array edge_length; + + edge_length = {mesh_bbox.max[0] - mesh_bbox.min[0], + mesh_bbox.max[1] - mesh_bbox.min[1]}; + + int nx = edge_length[0] * Kokkos::sqrt(mesh.nelems()); + int ny = edge_length[1] * Kokkos::sqrt(mesh.nelems()); + + pcms::GridPointSearch2D search_element(mesh, nx, ny); + auto results = search_element(points); + + return results; +} + +Omega_h::Reals evaluate_field_from_point_localization( + Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, + const Kokkos::View& results) +{ + + const auto& npoints = results.extent(0); + + const auto& faces2nodes = mesh.ask_down(Omega_h::FACE, Omega_h::VERT).ab2b; + const auto& coordinates = mesh.coords(); + + Omega_h::Write field_values_at_points( + npoints, 0.0, "stores field values at the given points"); + + Omega_h::parallel_for( + npoints, OMEGA_H_LAMBDA(int id) { + field_values_at_points[id] = 0.0; + + const int tri = results[id].tri_id; + if (tri < 0) + return; + const auto el_verts = Omega_h::gather_verts<3>(faces2nodes, tri); + + for (int j = 0; j < 3; ++j) { + const auto node_id = el_verts[j]; + field_values_at_points[id] += + nodal_field_values[node_id] * results[id].parametric_coords[j]; + } + }); + Kokkos::fence(); + return Omega_h::read(field_values_at_points); +} + +KOKKOS_INLINE_FUNCTION +double montecarlo_integration(const Omega_h::Real* shape_func_values_at_points, + const Omega_h::Real* src_values_at_points, + const int npoints_each_tri, const double volume) +{ + + double sum = 0; + + for (int i = 0; i < npoints_each_tri; ++i) { + sum += shape_func_values_at_points[i] * src_values_at_points[i]; + } + + sum *= volume; + return sum / npoints_each_tri; +} + +// TODO:: for black box coupling we wouldn't be provided with source mesh info +// and source field this works when source mesh and field are known +Kokkos::View loadVectorMCIntegrator( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, + const std::string& sobol_filename) +{ + + Kokkos::View ref_barycentric_coords; + if (method == SamplingMethod::SOBOL) { + if (sobol_filename.empty()) { + throw std::runtime_error("Could not open sobol sample file : "; + } + ref_barycentric_coords = + read_sobol_barycentric_samples_from_file(sobol_filename); + } else { + ref_barycentric_coords = + generate_uniform_random_barycentric_coords(npoints_each_tri); + } + + if (ref_barycentric_coords.extent(0) != + static_cast(npoints_each_tri)) { + throw std::runtime_error( + "Sobol sample count does not match npoints_each_tri."); + } + + int subVectorSize = 3; + + Omega_h::Reals elementsArea; + elementsArea = Omega_h::measure_elements_real(&target_mesh); + + Kokkos::View elmLoadVector("elmLoadVector", + target_mesh.nelems() * 3); + + Kokkos::parallel_for( + "eval load vector using MC", target_mesh.nelems(), + KOKKOS_LAMBDA(const int elm) { + Omega_h::Real N0[MAX_POINTS] = {}; + Omega_h::Real N1[MAX_POINTS] = {}; + Omega_h::Real N2[MAX_POINTS] = {}; + Omega_h::Real src_field_values[MAX_POINTS] = {}; + + int base_idx_src_values = elm * npoints_each_tri; + + for (int i = 0; i < npoints_each_tri; ++i) { + N0[i] = ref_barycentric_coords(i, 0); + N1[i] = ref_barycentric_coords(i, 1); + N2[i] = sampled_barycentric_coords(i, 2); + src_field_values[i] = field_values_at_points[base_idx_src_values + i]; + } + + Omega_h::Vector<3> result; + result[0] = montecarlo_integration(N0, src_field_values, npoints_each_tri, + elementsArea[elm]); + result[1] = montecarlo_integration(N1, src_field_values, npoints_each_tri, + elementsArea[elm]); + result[2] = montecarlo_integration(N2, src_field_values, npoints_each_tri, + elementsArea[elm]); + + for (int i = 0; i < 3; ++i) { + elmLoadVector(elm * subVectorSize + i) = result[i]; + } + }); + + return elmLoadVector; +} diff --git a/src/pcms/transfer/pcms_mc_proj_utils.hpp b/src/pcms/transfer/pcms_mc_proj_utils.hpp new file mode 100644 index 00000000..afd934f9 --- /dev/null +++ b/src/pcms/transfer/pcms_mc_proj_utils.hpp @@ -0,0 +1,228 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +constexpr int MAX_POINTS = 1600; +/** + * @brief Read precomputed Sobol barycentric samples from file. + * + * For now this routine reads barycentric Sobol samples from a text file with a + * header row. TODO: replace this with a function that directly generates Sobol + * sequence samples and maps them to barycentric coordinates. + * + * @param[in] file_path Path to the sample file. + * @return Device view of shape `(nsamples, 3)` containing barycentric samples. + */ +Kokkos::View read_sobol_barycentric_samples_from_file( + std::string file_path); + +/** + * @brief Generate uniform random barycentric coordinates for triangle sampling. + * + * Uses the standard square-to-triangle transform to produce barycentric + * coordinates uniformly distributed over triangle area. + * + * @param[in] npoints_each_tri Number of samples to generate. + * @return Device view of shape `(npoints_each_tri, 3)`. + */ +Kokkos::View generate_uniform_random_barycentric_coords( + const int npoints_each_tri); + +/** + * @brief Compute global sample coordinates in each element from reference + * barycentric sample coordinates. + * + * Reuses the same reference-triangle barycentric sample pattern in every target + * element and maps each barycentric sample to its corresponding physical + * coordinate using the element vertex coordinates. + * + * @param[in] target_mesh Target triangular mesh. + * @param[in] ref_barycentric_coords Reference barycentric coordinates of shape + * `(npoints_each_tri, 3)`. + * @return Device view of shape `(nelems * npoints_each_tri, 2)` containing the + * global sample coordinates in all elements. + */ +Kokkos::View global_coords_from_ref_barycentric_coords( + Omega_h::Mesh& target_mesh, + const Kokkos::View& ref_barycentric_coords); + +/** + * @brief Locate query points in a 2D triangular mesh. + * + * This routine performs point localization for a set of query points in the + * given 2D mesh using a structured grid search. For each query + * point, it returns the containing element id together with the associated + * parametric coordinates stored in `pcms::GridPointSearch::Result`. + * + * @param[in] mesh Input 2D mesh in which the query points are to be localized. + * @param[in] points Query point coordinates of shape `(npoints, 2)`. + * + * @return Search results for all query points. + * + * @note This routine assumes that `mesh` is two-dimensional. + */ +Kokkos::View localize_points_in_mesh( + Omega_h::Mesh& mesh, const Kokkos::View& points); + +/** + * @brief Evaluate a nodal field at localized query points. + * + * Uses the containing element id and parametric or barycentric coordinates from + * a point-localization step to interpolate the nodal field values at the query + * points. + * + * @param[in] mesh Input triangular mesh. + * @param[in] nodal_field_values Field values stored at mesh vertices. + * @param[in] results Point-localization results for the query points. + * @return Field values evaluated at the query points. + * + * @note Points with invalid element ids retain the default value `0.0`. + */ +Omega_h::Reals evaluate_field_from_point_localization( + Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, + const Kokkos::View& points, + const Kokkos::View& results); + +/** + * @brief Compute a Monte Carlo estimate of an element integral. + * + * Approximates the integral of the product of shape-function values and source + * values over an element using uniformly distributed sample points. + * + * @param[in] shape_func_values_at_points Shape-function values at sample + * points. + * @param[in] src_values_at_points Source-field values at the same sample + * points. + * @param[in] npoints_each_tri Number of sample points. + * @param[in] volume Element volume. + * @return Estimated integral value. + */ +KOKKOS_INLINE_FUNCTION +double monte_carlo_integral(const Omega_h::Real* shape_func_values_at_points, + const Omega_h::Real* src_values_at_points, + const int npoints_each_tri, + const Omega_h::Real volume) +{ + + if (npoints_each_tri <= 0) + return 0.0; + Omega_h::Real sum = 0; + + for (int i = 0; i < npoints_each_tri; ++i) { + sum += shape_func_values_at_points[i] * src_values_at_points[i]; + } + + sum *= volume; + return sum / npoints_each_tri; +} + +/** + * @brief Compute element-local load vectors using Monte Carlo integration. + * + * For each target triangle \f$\Omega_t\f$, this routine approximates the local + * load-vector entries by uniform Monte Carlo sampling: + * \f[ + * \widehat b_k^{\,t} + * = \frac{|\Omega_t|}{N}\sum_{i=1}^N + * f^s(\mathbf X_{t,i})\,\psi_k(\mathbf X_{t,i}), + * \f] + * where \f$N\f$ is the number of sample points in the element, + * \f$f^s(\mathbf X_{t,i})\f$ is the source-field value at the sampled point, + * and \f$\psi_k\f$ is the local target basis function. + * + * The sampling points are generated on the reference triangle in barycentric + * coordinates and the same reference sample set is reused for all target + * elements. For linear triangular elements, the barycentric coordinates are + * equal to the local shape-function values, so they are used directly in the + * Monte Carlo estimator. + * + * @param[in] target_mesh Target 2D triangular mesh. + * @param[in] field_values_at_points Source-field values evaluated at the + * sampled physical points, stored element-by-element. + * @param[in] npoints_each_tri Number of Monte Carlo sample points per target + * triangle. + * @param[in] method Sampling method used to generate the reference barycentric + * coordinates. + * @param[in] sobol_filename File containing precomputed Sobol barycentric + * samples when Sobol sampling is selected. + * + * @return Flattened element-local load vectors of size + * `3 * target_mesh.nelems()`. + * + * @note This routine computes element-local contributions only; it does not + * assemble a global load vector. + * @note This routine assumes that `field_values_at_points` has already been + * evaluated at the sampled physical points. + */ + +Kokkos::View loadVectorMCIntegrator( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, + const std::string& sobol_filename) +{ + + Kokkos::View ref_barycentric_coords; + if (method == SamplingMethod::SOBOL) { + if (sobol_filename.empty()) { + throw std::runtime_error("Could not open sobol sample file : "; + } + ref_barycentric_coords = + read_sobol_barycentric_samples_from_file(sobol_filename); + } else { + ref_barycentric_coords = + generate_uniform_random_barycentric_coords(npoints_each_tri); + } + + int subVectorSize = 3; + + Omega_h::Reals elementsArea; + elementsArea = Omega_h::measure_elements_real(&target_mesh); + + Kokkos::View elmLoadVector("elmLoadVector", + target_mesh.nelems() * 3); + + Kokkos::parallel_for( + "eval load vector using MC", target_mesh.nelems(), + KOKKOS_LAMBDA(const int elm) { + Omega_h::Real N0[MAX_POINTS] = {}; + Omega_h::Real N1[MAX_POINTS] = {}; + Omega_h::Real N2[MAX_POINTS] = {}; + Omega_h::Real src_field_values[MAX_POINTS] = {}; + + int base_idx_src_values = elm * npoints_each_tri; + + for (int i = 0; i < npoints_each_tri; ++i) { + N0[i] = ref_barycentric_coords(i, 0); + N1[i] = ref_barycentric_coords(i, 1); + N2[i] = sampled_barycentric_coords(i, 2); + src_field_values[i] = field_values_at_points[base_idx_src_values + i]; + } + + Omega_h::Vector<3> result; + result[0] = montecarlo_integration(N0, src_field_values, npoints_each_tri, + elementsArea[elm]); + result[1] = montecarlo_integration(N1, src_field_values, npoints_each_tri, + elementsArea[elm]); + result[2] = montecarlo_integration(N2, src_field_values, npoints_each_tri, + elementsArea[elm]); + + for (int i = 0; i < 3; ++i) { + elmLoadVector(elm * subVectorSize + i) = result[i]; + } + }); + + return elmLoadVector; +} From d2eb0af6d9471cb68e71bd35fa3142e1d3388de9 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:06:27 -0400 Subject: [PATCH 02/44] add mc helper functions declaration --- src/pcms/transfer/load_vector_integrator.hpp | 82 ++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index 4e1d2fad..eb7f12be 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -139,6 +139,7 @@ struct IntegrationData Kokkos::View buildLoadVector( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values); + /// Holds projection and conservation error metrics returned by /// evaluate_pro_and_cons_errors(). struct Errors @@ -215,6 +216,87 @@ Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, const IntersectionResults& intersection, const Omega_h::Reals& target_values, const Omega_h::Reals& source_values); + +/** + * @brief Read precomputed Sobol barycentric samples from file. + * + * For now this routine reads barycentric Sobol samples from a text file with a + * header row. TODO: replace this with a function that directly generates Sobol + * sequence samples and maps them to barycentric coordinates. + * + * @param[in] file_path Path to the sample file. + * @return Device view of shape `(nsamples, 3)` containing barycentric samples. + */ +Kokkos::View read_sobol_barycentric_samples_from_file( + std::string file_path); + +/** + * @brief Generate uniform random barycentric coordinates for triangle sampling. + * + * Uses the standard square-to-triangle transform to produce barycentric + * coordinates uniformly distributed over triangle area. + * + * @param[in] npoints_each_tri Number of samples to generate. + * @return Device view of shape `(npoints_each_tri, 3)`. + */ +Kokkos::View generate_uniform_random_barycentric_coords( + const int npoints_each_tri); + +/** + * @brief Compute global sample coordinates in each element from reference + * barycentric sample coordinates. + * + * Reuses the same reference-triangle barycentric sample pattern in every target + * element and maps each barycentric sample to its corresponding physical + * coordinate using the element vertex coordinates. + * + * @param[in] target_mesh Target triangular mesh. + * @param[in] ref_barycentric_coords Reference barycentric coordinates of shape + * `(npoints_each_tri, 3)`. + * @return Device view of shape `(nelems * npoints_each_tri, 2)` containing the + * global sample coordinates in all elements. + */ +Kokkos::View global_coords_from_ref_barycentric_coords( + Omega_h::Mesh& target_mesh, + const Kokkos::View& ref_barycentric_coords); + +/** + * @brief Locate query points in a 2D triangular mesh. + * + * This routine performs point localization for a set of query points in the + * given 2D mesh using a structured grid search. For each query + * point, it returns the containing element id together with the associated + * parametric coordinates stored in `pcms::GridPointSearch::Result`. + * + * @param[in] mesh Input 2D mesh in which the query points are to be localized. + * @param[in] points Query point coordinates of shape `(npoints, 2)`. + * + * @return Search results for all query points. + * + * @note This routine assumes that `mesh` is two-dimensional. + */ +Kokkos::View localize_points_in_mesh( + Omega_h::Mesh& mesh, const Kokkos::View& points); + +/** + * @brief Evaluate a nodal field at localized query points. + * + * Uses the containing element id and parametric or barycentric coordinates from + * a point-localization step to interpolate the nodal field values at the query + * points. + * + * @param[in] mesh Input triangular mesh. + * @param[in] nodal_field_values Field values stored at mesh vertices. + * @param[in] results Point-localization results for the query points. + * @return Field values evaluated at the query points. + * + * @note Points with invalid element ids retain the default value `0.0`. + */ +Omega_h::Reals evaluate_field_from_point_localization( + Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, + const Kokkos::View& points, + const Kokkos::View& results); + } // namespace pcms #endif // PCMS_TRANSFER_LOAD_VECTOR_INTEGRATOR_HPP From 8f5148af40f7dc623c795bd00a9488146748079c Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:10:14 -0400 Subject: [PATCH 03/44] change function name from buildLoadVector to buildLoadVectorMI --- src/pcms/transfer/load_vector_integrator.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator.cpp b/src/pcms/transfer/load_vector_integrator.cpp index 477f744b..efb25969 100644 --- a/src/pcms/transfer/load_vector_integrator.cpp +++ b/src/pcms/transfer/load_vector_integrator.cpp @@ -297,7 +297,7 @@ OMEGA_H_INLINE void for_each_intersection_subtriangle( } } -Kokkos::View buildLoadVector( +Kokkos::View buildLoadVectorMI( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values) { @@ -357,6 +357,7 @@ Kokkos::View buildLoadVector( return elmLoadVector; } + Errors evaluate_proj_and_cons_errors(Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, From 9094916f589a45bc50a8b967c9d1ee4a908ae4c5 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:11:10 -0400 Subject: [PATCH 04/44] fix the load vector integrator name --- src/pcms/transfer/load_vector_integrator.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index eb7f12be..e4b6a331 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -136,7 +136,7 @@ struct IntegrationData * @see IntersectionResults */ -Kokkos::View buildLoadVector( +Kokkos::View buildLoadVectorMI( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values); From a95460ba2ad0613894c8e76e939d89fefe562c14 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:15:26 -0400 Subject: [PATCH 05/44] add mc load vector build declaration --- src/pcms/transfer/load_vector_integrator.hpp | 55 ++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index e4b6a331..aeb416a0 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -23,7 +23,18 @@ #include #include #include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace pcms { @@ -297,6 +308,50 @@ Omega_h::Reals evaluate_field_from_point_localization( const Kokkos::View& points, const Kokkos::View& results); +/** + * @brief Compute element-local load vectors using Monte Carlo integration. + * + * For each target triangle \f$\Omega_t\f$, this routine approximates the local + * load-vector entries by uniform Monte Carlo sampling: + * \f[ + * \widehat b_k^{\,t} + * = \frac{|\Omega_t|}{N}\sum_{i=1}^N + * f^s(\mathbf X_{t,i})\,\psi_k(\mathbf X_{t,i}), + * \f] + * where \f$N\f$ is the number of sample points in the element, + * \f$f^s(\mathbf X_{t,i})\f$ is the source-field value at the sampled point, + * and \f$\psi_k\f$ is the local target basis function. + * + * The sampling points are generated on the reference triangle in barycentric + * coordinates and the same reference sample set is reused for all target + * elements. For linear triangular elements, the barycentric coordinates are + * equal to the local shape-function values, so they are used directly in the + * Monte Carlo estimator. + * + * @param[in] target_mesh Target 2D triangular mesh. + * @param[in] field_values_at_points Source-field values evaluated at the + * sampled physical points, stored element-by-element. + * @param[in] npoints_each_tri Number of Monte Carlo sample points per target + * triangle. + * @param[in] method Sampling method used to generate the reference barycentric + * coordinates. + * @param[in] sobol_filename File containing precomputed Sobol barycentric + * samples when Sobol sampling is selected. + * + * @return Flattened element-local load vectors of size + * `3 * target_mesh.nelems()`. + * + * @note This routine computes element-local contributions only; it does not + * assemble a global load vector. + * @note This routine assumes that `field_values_at_points` has already been + * evaluated at the sampled physical points. + */ + +Kokkos::View buildLoadVectorMC( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, + const std::string& sobol_filename) + } // namespace pcms #endif // PCMS_TRANSFER_LOAD_VECTOR_INTEGRATOR_HPP From e0b9fcd3fd5872c6973125f502d3f492e379b519 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:17:43 -0400 Subject: [PATCH 06/44] add header --- src/pcms/transfer/pcms_mc_proj_utils.cpp | 53 ++++++++---------------- 1 file changed, 18 insertions(+), 35 deletions(-) diff --git a/src/pcms/transfer/pcms_mc_proj_utils.cpp b/src/pcms/transfer/pcms_mc_proj_utils.cpp index 152ff9c9..9a7c9320 100644 --- a/src/pcms/transfer/pcms_mc_proj_utils.cpp +++ b/src/pcms/transfer/pcms_mc_proj_utils.cpp @@ -1,6 +1,5 @@ #include - -constexpr int MAX_POINTS = 1600; +#include "pcms/transfer/load_vector_integrator.hpp" // TODO:: create a function that can sample sobol sequences instead of // generating in python and reading and using @@ -174,24 +173,27 @@ Omega_h::Reals evaluate_field_from_point_localization( } KOKKOS_INLINE_FUNCTION -double montecarlo_integration(const Omega_h::Real* shape_func_values_at_points, - const Omega_h::Real* src_values_at_points, - const int npoints_each_tri, const double volume) +Omega_h::Vector<3> compute_elm_load_vector_mc( + const Kokkos::View& ref_barycentric_coords, + const Omega_h::Reals& field_values_at_points, + const Omega_h::Real& elmVolume const int elm) { - double sum = 0; - + const int npoints_each_tri = ref_barycentric_coords.extent(0); + Omega_h::Real sum0 = 0.0; + Omega_h::Real sum1 = 0.0; + Omega_h::Real sum2 = 0.0; + const int base = elm * npoints_each_tri; for (int i = 0; i < npoints_each_tri; ++i) { - sum += shape_func_values_at_points[i] * src_values_at_points[i]; + const auto f = field_values_at_points[base + i]; + sum0 += ref_barycentric_coords(i, 0) * f; + sum1 += ref_barycentric_coords(i, 1) * f; + sum2 += ref_barycentric_coords(i, 2) * f; } - - sum *= volume; - return sum / npoints_each_tri; + return elmVolume / npoints_each_tri * {sum0, sum1, sum2}; } -// TODO:: for black box coupling we wouldn't be provided with source mesh info -// and source field this works when source mesh and field are known -Kokkos::View loadVectorMCIntegrator( +Kokkos::View buildLoadVectorMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, const int npoints_each_tri, SamplingMethod method, const std::string& sobol_filename) @@ -226,27 +228,8 @@ Kokkos::View loadVectorMCIntegrator( Kokkos::parallel_for( "eval load vector using MC", target_mesh.nelems(), KOKKOS_LAMBDA(const int elm) { - Omega_h::Real N0[MAX_POINTS] = {}; - Omega_h::Real N1[MAX_POINTS] = {}; - Omega_h::Real N2[MAX_POINTS] = {}; - Omega_h::Real src_field_values[MAX_POINTS] = {}; - - int base_idx_src_values = elm * npoints_each_tri; - - for (int i = 0; i < npoints_each_tri; ++i) { - N0[i] = ref_barycentric_coords(i, 0); - N1[i] = ref_barycentric_coords(i, 1); - N2[i] = sampled_barycentric_coords(i, 2); - src_field_values[i] = field_values_at_points[base_idx_src_values + i]; - } - - Omega_h::Vector<3> result; - result[0] = montecarlo_integration(N0, src_field_values, npoints_each_tri, - elementsArea[elm]); - result[1] = montecarlo_integration(N1, src_field_values, npoints_each_tri, - elementsArea[elm]); - result[2] = montecarlo_integration(N2, src_field_values, npoints_each_tri, - elementsArea[elm]); + auto result = compute_elm_load_vector_mc( + ref_barycentric_coords, field_values_at_points, elementsArea[elm], elm); for (int i = 0; i < 3; ++i) { elmLoadVector(elm * subVectorSize + i) = result[i]; From 26f1a4fc31ce48ea67942f4f231070b5133093d2 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:19:38 -0400 Subject: [PATCH 07/44] delete pcms_mc_proj_utils.hpp --- src/pcms/transfer/pcms_mc_proj_utils.hpp | 228 ----------------------- 1 file changed, 228 deletions(-) delete mode 100644 src/pcms/transfer/pcms_mc_proj_utils.hpp diff --git a/src/pcms/transfer/pcms_mc_proj_utils.hpp b/src/pcms/transfer/pcms_mc_proj_utils.hpp deleted file mode 100644 index afd934f9..00000000 --- a/src/pcms/transfer/pcms_mc_proj_utils.hpp +++ /dev/null @@ -1,228 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -constexpr int MAX_POINTS = 1600; -/** - * @brief Read precomputed Sobol barycentric samples from file. - * - * For now this routine reads barycentric Sobol samples from a text file with a - * header row. TODO: replace this with a function that directly generates Sobol - * sequence samples and maps them to barycentric coordinates. - * - * @param[in] file_path Path to the sample file. - * @return Device view of shape `(nsamples, 3)` containing barycentric samples. - */ -Kokkos::View read_sobol_barycentric_samples_from_file( - std::string file_path); - -/** - * @brief Generate uniform random barycentric coordinates for triangle sampling. - * - * Uses the standard square-to-triangle transform to produce barycentric - * coordinates uniformly distributed over triangle area. - * - * @param[in] npoints_each_tri Number of samples to generate. - * @return Device view of shape `(npoints_each_tri, 3)`. - */ -Kokkos::View generate_uniform_random_barycentric_coords( - const int npoints_each_tri); - -/** - * @brief Compute global sample coordinates in each element from reference - * barycentric sample coordinates. - * - * Reuses the same reference-triangle barycentric sample pattern in every target - * element and maps each barycentric sample to its corresponding physical - * coordinate using the element vertex coordinates. - * - * @param[in] target_mesh Target triangular mesh. - * @param[in] ref_barycentric_coords Reference barycentric coordinates of shape - * `(npoints_each_tri, 3)`. - * @return Device view of shape `(nelems * npoints_each_tri, 2)` containing the - * global sample coordinates in all elements. - */ -Kokkos::View global_coords_from_ref_barycentric_coords( - Omega_h::Mesh& target_mesh, - const Kokkos::View& ref_barycentric_coords); - -/** - * @brief Locate query points in a 2D triangular mesh. - * - * This routine performs point localization for a set of query points in the - * given 2D mesh using a structured grid search. For each query - * point, it returns the containing element id together with the associated - * parametric coordinates stored in `pcms::GridPointSearch::Result`. - * - * @param[in] mesh Input 2D mesh in which the query points are to be localized. - * @param[in] points Query point coordinates of shape `(npoints, 2)`. - * - * @return Search results for all query points. - * - * @note This routine assumes that `mesh` is two-dimensional. - */ -Kokkos::View localize_points_in_mesh( - Omega_h::Mesh& mesh, const Kokkos::View& points); - -/** - * @brief Evaluate a nodal field at localized query points. - * - * Uses the containing element id and parametric or barycentric coordinates from - * a point-localization step to interpolate the nodal field values at the query - * points. - * - * @param[in] mesh Input triangular mesh. - * @param[in] nodal_field_values Field values stored at mesh vertices. - * @param[in] results Point-localization results for the query points. - * @return Field values evaluated at the query points. - * - * @note Points with invalid element ids retain the default value `0.0`. - */ -Omega_h::Reals evaluate_field_from_point_localization( - Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, - const Kokkos::View& points, - const Kokkos::View& results); - -/** - * @brief Compute a Monte Carlo estimate of an element integral. - * - * Approximates the integral of the product of shape-function values and source - * values over an element using uniformly distributed sample points. - * - * @param[in] shape_func_values_at_points Shape-function values at sample - * points. - * @param[in] src_values_at_points Source-field values at the same sample - * points. - * @param[in] npoints_each_tri Number of sample points. - * @param[in] volume Element volume. - * @return Estimated integral value. - */ -KOKKOS_INLINE_FUNCTION -double monte_carlo_integral(const Omega_h::Real* shape_func_values_at_points, - const Omega_h::Real* src_values_at_points, - const int npoints_each_tri, - const Omega_h::Real volume) -{ - - if (npoints_each_tri <= 0) - return 0.0; - Omega_h::Real sum = 0; - - for (int i = 0; i < npoints_each_tri; ++i) { - sum += shape_func_values_at_points[i] * src_values_at_points[i]; - } - - sum *= volume; - return sum / npoints_each_tri; -} - -/** - * @brief Compute element-local load vectors using Monte Carlo integration. - * - * For each target triangle \f$\Omega_t\f$, this routine approximates the local - * load-vector entries by uniform Monte Carlo sampling: - * \f[ - * \widehat b_k^{\,t} - * = \frac{|\Omega_t|}{N}\sum_{i=1}^N - * f^s(\mathbf X_{t,i})\,\psi_k(\mathbf X_{t,i}), - * \f] - * where \f$N\f$ is the number of sample points in the element, - * \f$f^s(\mathbf X_{t,i})\f$ is the source-field value at the sampled point, - * and \f$\psi_k\f$ is the local target basis function. - * - * The sampling points are generated on the reference triangle in barycentric - * coordinates and the same reference sample set is reused for all target - * elements. For linear triangular elements, the barycentric coordinates are - * equal to the local shape-function values, so they are used directly in the - * Monte Carlo estimator. - * - * @param[in] target_mesh Target 2D triangular mesh. - * @param[in] field_values_at_points Source-field values evaluated at the - * sampled physical points, stored element-by-element. - * @param[in] npoints_each_tri Number of Monte Carlo sample points per target - * triangle. - * @param[in] method Sampling method used to generate the reference barycentric - * coordinates. - * @param[in] sobol_filename File containing precomputed Sobol barycentric - * samples when Sobol sampling is selected. - * - * @return Flattened element-local load vectors of size - * `3 * target_mesh.nelems()`. - * - * @note This routine computes element-local contributions only; it does not - * assemble a global load vector. - * @note This routine assumes that `field_values_at_points` has already been - * evaluated at the sampled physical points. - */ - -Kokkos::View loadVectorMCIntegrator( - Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, - const int npoints_each_tri, SamplingMethod method, - const std::string& sobol_filename) -{ - - Kokkos::View ref_barycentric_coords; - if (method == SamplingMethod::SOBOL) { - if (sobol_filename.empty()) { - throw std::runtime_error("Could not open sobol sample file : "; - } - ref_barycentric_coords = - read_sobol_barycentric_samples_from_file(sobol_filename); - } else { - ref_barycentric_coords = - generate_uniform_random_barycentric_coords(npoints_each_tri); - } - - int subVectorSize = 3; - - Omega_h::Reals elementsArea; - elementsArea = Omega_h::measure_elements_real(&target_mesh); - - Kokkos::View elmLoadVector("elmLoadVector", - target_mesh.nelems() * 3); - - Kokkos::parallel_for( - "eval load vector using MC", target_mesh.nelems(), - KOKKOS_LAMBDA(const int elm) { - Omega_h::Real N0[MAX_POINTS] = {}; - Omega_h::Real N1[MAX_POINTS] = {}; - Omega_h::Real N2[MAX_POINTS] = {}; - Omega_h::Real src_field_values[MAX_POINTS] = {}; - - int base_idx_src_values = elm * npoints_each_tri; - - for (int i = 0; i < npoints_each_tri; ++i) { - N0[i] = ref_barycentric_coords(i, 0); - N1[i] = ref_barycentric_coords(i, 1); - N2[i] = sampled_barycentric_coords(i, 2); - src_field_values[i] = field_values_at_points[base_idx_src_values + i]; - } - - Omega_h::Vector<3> result; - result[0] = montecarlo_integration(N0, src_field_values, npoints_each_tri, - elementsArea[elm]); - result[1] = montecarlo_integration(N1, src_field_values, npoints_each_tri, - elementsArea[elm]); - result[2] = montecarlo_integration(N2, src_field_values, npoints_each_tri, - elementsArea[elm]); - - for (int i = 0; i < 3; ++i) { - elmLoadVector(elm * subVectorSize + i) = result[i]; - } - }); - - return elmLoadVector; -} From 12c50c7e850f08d9f5444b72ce37bf5ff36bf239 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:22:01 -0400 Subject: [PATCH 08/44] rename pcms_mc_proj_utils.cpp -> load_vector_integrator_mc.cpp --- .../{pcms_mc_proj_utils.cpp => load_vector_integrator_mc.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/pcms/transfer/{pcms_mc_proj_utils.cpp => load_vector_integrator_mc.cpp} (100%) diff --git a/src/pcms/transfer/pcms_mc_proj_utils.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp similarity index 100% rename from src/pcms/transfer/pcms_mc_proj_utils.cpp rename to src/pcms/transfer/load_vector_integrator_mc.cpp From 21fdabad8ec088dedbed7ac4b856808f8a22f04b Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:25:32 -0400 Subject: [PATCH 09/44] rename load_vector_integrator.cpp -> load_vector_integrator_mi.cpp --- .../{load_vector_integrator.cpp => load_vector_integrator_mi.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/pcms/transfer/{load_vector_integrator.cpp => load_vector_integrator_mi.cpp} (100%) diff --git a/src/pcms/transfer/load_vector_integrator.cpp b/src/pcms/transfer/load_vector_integrator_mi.cpp similarity index 100% rename from src/pcms/transfer/load_vector_integrator.cpp rename to src/pcms/transfer/load_vector_integrator_mi.cpp From ab31a90d2dcea9271ed3713c8e0b33f951175b6b Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 16:53:16 -0400 Subject: [PATCH 10/44] add mc load vector build --- src/pcms/transfer/calculate_load_vector.hpp | 75 +++++++++++++++++---- 1 file changed, 62 insertions(+), 13 deletions(-) diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp index b1968338..43f70a3b 100644 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ b/src/pcms/transfer/calculate_load_vector.hpp @@ -15,12 +15,15 @@ #include #include +namespace pcms +{ + /** - * @brief Assembles the global load vector. + * @brief Assembles the global load vector using mesh intersection method. * * This function computes the unassembled local load vector contributions for - * each triangular element in the target mesh using `buildLoadVector()` and then - * assembles them into a global PETSc vector in COO format. + * each triangular element in the target mesh using `buildLoadVectorMI()` and + * then assembles them into a global PETSc vector in COO format. * * * @param target_mesh The target Omega_h mesh to which the scalar field is being @@ -40,20 +43,66 @@ * @note * - Works for 2D linear triangular elements. * - Uses COO-style preallocation and insertion into the PETSc vector. - * - Internally calls `buildLoadVector()` to compute per-element contributions. + * - Internally calls `buildLoadVectorMI()` to compute per-element + * contributions. * - The resulting vector is used as the right-hand side (RHS) in a projection * solve. * - * @see buildLoadVector,IntersectionResults + * @see buildLoadVectorMI,IntersectionResults */ - -namespace pcms -{ // FIXME use PCMS error handling rather than returning a PETSC error code -PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values, - Vec* loadVec_out); +PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& source_values, + Vec* loadVec_out); + +/** + * @brief Assembles the global load vector using Monte Carlo integration. + * + * This function computes the unassembled local load vector contributions for + * each triangular element in the target mesh using Monte Carlo integration and + * then assembles them into a global PETSc vector in COO format. + * + * The element-local contributions are computed from source-field values already + * evaluated at the sampled physical points in each target element. The sampling + * pattern is defined on the reference triangle and may be generated either by + * uniform random sampling or by precomputed Sobol barycentric samples. + * + * @param target_mesh The target Omega_h mesh to which the scalar field is being + * projected. + * @param field_values_at_points Source-field values evaluated at the sampled + * physical points in the target mesh. The data is expected to be stored + * element-by-element, with total size + * `target_mesh.nelems() * npoints_each_tri`. + * @param npoints_each_tri Number of sample points used in each target triangle + * for the Monte Carlo integration. + * @param method Sampling method used to define the reference barycentric sample + * coordinates. Supported options include uniform random sampling and Sobol + * sampling. + * @param sobol_filename Path to the file containing precomputed Sobol + * barycentric samples when `method` is `SamplingMethod::SOBOL`. + * @param[out] loadVec_out Pointer to a PETSc Vec where the assembled load + * vector will be stored. + * + * @return PetscErrorCode Returns PETSC_SUCCESS if successful, or an appropriate + * PETSc error code otherwise. + * + * @note + * - Works for 2D linear triangular elements. + * - Uses COO-style preallocation and insertion into the PETSc vector. + * - Internally computes per-element Monte Carlo load vector contributions + * before assembling the global vector. + * - For linear triangular elements, the barycentric coordinates are equal to + * the local shape-function values used in the Monte Carlo estimator. + * - The resulting vector is used as the right-hand side (RHS) in a projection + * solve. + * + * @see buildLoadVectorMC, SamplingMethod + */ +PetscErrorCode calculateLoadVectorMC( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, + const std::string sobol_filename = "", Vec* loadVec_out); } // namespace pcms #endif // PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP From 54213bb7824bb6d4aa49f3d462b9029f563eb7b4 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 17:22:29 -0400 Subject: [PATCH 11/44] add build load vector with mc approach --- src/pcms/transfer/calculate_load_vector.cpp | 47 ++++++++++++++++++++- 1 file changed, 46 insertions(+), 1 deletion(-) diff --git a/src/pcms/transfer/calculate_load_vector.cpp b/src/pcms/transfer/calculate_load_vector.cpp index f351a4ef..6da01206 100644 --- a/src/pcms/transfer/calculate_load_vector.cpp +++ b/src/pcms/transfer/calculate_load_vector.cpp @@ -21,7 +21,52 @@ PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh, // Fill COO values auto elmLoadVector = - buildLoadVector(target_mesh, source_mesh, intersection, source_values); + buildLoadVectorMI(target_mesh, source_mesh, intersection, source_values); + + auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector); + Kokkos::deep_copy(hostElmLoadVector, elmLoadVector); + PetscCheck(static_cast(hostElmLoadVector.extent(0)) == nnz, + PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, + "Element load vector size (%d) does not match COO nnz (%d)", + static_cast(hostElmLoadVector.extent(0)), nnz); + + for (PetscInt e = 0; e < nnz; ++e) { + coo_vals[e] = hostElmLoadVector(e); + } + + // create vector with preallocated COO structure + Vec vec; + PetscCall(createSeqVec(PETSC_COMM_WORLD, target_mesh.nverts(), &vec)); + PetscCall(VecSetPreallocationCOO(vec, nnz, coo_i)); + PetscCall(VecSetValuesCOO(vec, coo_vals, ADD_VALUES)); + PetscCall(PetscFree(coo_i)); + PetscCall(PetscFree(coo_vals)); + + if (target_mesh.nelems() < 10) { + PetscCall(VecView(vec, PETSC_VIEWER_STDOUT_WORLD)); + } + + *loadVec_out = vec; + PetscFunctionReturn(PETSC_SUCCESS); +} + +PetscErrorCode calculateLoadVectorMC( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out, + const std::string sobol_filename) +{ + + PetscFunctionBeginUser; + PetscInt nnz = 0; + PetscInt* coo_i = nullptr; + PetscCall(build_linear_triangle_coo_rows(target_mesh, &coo_i, &nnz)); + PetscScalar* coo_vals = nullptr; + PetscCall(PetscMalloc1(nnz, &coo_vals)); + + // Fill COO values + auto elmLoadVector = + buildLoadVectorMC(target_mesh, field_values_at_points, npoints_each_tri, + method, sobol_filename); auto hostElmLoadVector = Kokkos::create_mirror_view(elmLoadVector); Kokkos::deep_copy(hostElmLoadVector, elmLoadVector); From ecd82348d069e86de97ac3d8c88deffb0725ad9d Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 17:43:56 -0400 Subject: [PATCH 12/44] change the parameter sequence --- src/pcms/transfer/calculate_load_vector.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp index 43f70a3b..0533b72e 100644 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ b/src/pcms/transfer/calculate_load_vector.hpp @@ -102,7 +102,7 @@ PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh, */ PetscErrorCode calculateLoadVectorMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, - const int npoints_each_tri, SamplingMethod method, - const std::string sobol_filename = "", Vec* loadVec_out); + const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out, + const std::string sobol_filename = ""); } // namespace pcms #endif // PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP From 56c01605fec973c8d801a9687171cec7af8eae8b Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 17:45:18 -0400 Subject: [PATCH 13/44] add mc solver function --- src/pcms/transfer/conservative_projection_solver.hpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index d3e46a37..c8f9868c 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -58,10 +58,14 @@ namespace pcms * */ -Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values); +Omega_h::Reals solveGalerkinProjectionMI( + Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, const Omega_h::Reals& source_values); + +Omega_h::Reals solveGalerkinProjectionMC( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out, + const std::string sobol_filename = ""); Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, From e489bea65d41bad7a2c587a376d2afd9152b1794 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Tue, 17 Mar 2026 17:45:38 -0400 Subject: [PATCH 14/44] change parameter sequence --- src/pcms/transfer/load_vector_integrator_mc.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 9a7c9320..1109dac3 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -1,4 +1,3 @@ -#include #include "pcms/transfer/load_vector_integrator.hpp" // TODO:: create a function that can sample sobol sequences instead of From 17644e3f1fae1d47a6dd945fd2e16c7cf1aff715 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:12:38 -0400 Subject: [PATCH 15/44] update variants of galerkin projection solve --- .../conservative_projection_solver.cpp | 83 ++++++++- .../conservative_projection_solver.hpp | 172 +++++++++++++++--- 2 files changed, 216 insertions(+), 39 deletions(-) diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index 148454f8..35d788e3 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -86,10 +86,9 @@ static Omega_h::Reals vecToOmegaHReals(Vec vec) return Omega_h::Reals(values_host); } -Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values) +Omega_h::Reals solveGalerkinProjectionMI( + Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, const Omega_h::Reals& source_values) { if ((PetscInt)source_values.size() != source_mesh.coords().size() / source_mesh.dim()) { @@ -105,6 +104,46 @@ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, CHKERRABORT(PETSC_COMM_WORLD, ierr); Vec vec; + ierr = calculateLoadVectorMI(target_mesh, source_mesh, intersection, + source_values, &vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Vec x = solveLinearSystem(mass, vec); + auto solution_vector = vecToOmegaHReals(x); + + ierr = VecDestroy(&x); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = MatDestroy(&mass); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + ierr = VecDestroy(&vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + return solution_vector; +} + +Omega_h::Reals solveGalerkinProjectionMC( + Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, SamplingMethod method, + const std::string sobol_filename); +{ + if ((PetscInt)source_values.size() != + source_mesh.coords().size() / source_mesh.dim()) { + std::cerr << "ERROR: source_values size (" << source_values.size() + << ") doesn't match expected size (" + << source_mesh.coords().size() / source_mesh.dim() << ")" + << std::endl; + throw std::runtime_error("source_values length mismatch"); + } + + Mat mass; + PetscErrorCode ierr = calculateMassMatrix(target_mesh, &mass); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + Vec vec; + ierr = calculateLoadVectorMC(target_mesh, field_values_at_points, + npoints_each_tri, method, &vec, sobol_filename); ierr = calculateLoadVector(target_mesh, source_mesh, intersection, source_values, &vec); CHKERRABORT(PETSC_COMM_WORLD, ierr); @@ -123,15 +162,16 @@ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, return solution_vector; } -Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values) + +Omega_h::Reals computeRhsVectorMI(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& source_values) { Vec vec; PetscErrorCode ierr; - ierr = calculateLoadVector(target_mesh, source_mesh, intersection, - source_values, &vec); + ierr = calculateLoadVectorMI(target_mesh, source_mesh, intersection, + source_values, &vec); CHKERRABORT(PETSC_COMM_WORLD, ierr); auto rhsvector = vecToOmegaHReals(vec); @@ -141,4 +181,27 @@ Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, return rhsvector; } + +Omega_h::Reals computeRhsVectorMC(Omega_h::Mesh& target_mesh, + const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, + SamplingMethod method, + const std::string& sobol_filename = "") +{ + + { + Vec vec; + PetscErrorCode ierr; + ierr = + calculateLoadVectorMC(target_mesh, field_values_at_points, + npoints_each_tri, method, &vec, sobol_filename); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + auto rhsvector = vecToOmegaHReals(vec); + ierr = VecDestroy(&vec); + CHKERRABORT(PETSC_COMM_WORLD, ierr); + + return rhsvector; + } + } // namespace pcms diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index c8f9868c..2d3feee8 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -25,52 +25,166 @@ namespace pcms { /** - * @brief Solves a conservative galerkin projection problem to transfer scalar - * field values onto a target mesh. + * @brief Solves a conservative Galerkin projection problem to transfer a scalar + * field onto a target mesh. * - * This function assembles and solves a linear system of the form: + * This routine assembles and solves a linear system of the form * \f[ - * M \cdot x = f + * M \cdot x = b * \f] - * where: - * - \f$M\f$ is the mass matrix on the target mesh (based on P1 finite - * elements), - * - \f$f\f$ is the load vector computed on the supermesh, - * - \f$x\f$ is the unknown nodal field on the target mesh (solution). - * - * The method computes the conservative field transfer between two non-matching - * meshes using mesh intersections (supermesh). - * - * ### Algorithm Steps: - * 1. Compute and assemble mass matrix and load vector - * 2. Solve the linear system using PETSc. - * 3. Return the solution as a nodal field on the target mesh. - * - * @param target_mesh The Omega_h mesh where the field is projected. - * @param source_mesh The Omega_h mesh containing the original field data. - * @param intersection Precomputed intersection information between source and - * target meshes. - * @param source_values Nodal scalar field values on the source mesh. + * where + * - \f$M\f$ is the mass matrix on the target mesh based on linear (\f$P_1\f$) + * finite elements, + * - \f$b\f$ is the load vector associated with the projection, + * - \f$x\f$ is the unknown nodal field on the target mesh. * - * @return A vector of nodal values on the target mesh after projection - * (Omega_h::Reals). + * The load vector is computed by integrating the source field against the + * target basis functions. Depending on the solver variant, these integrals are + * evaluated either exactly using mesh-intersection method or approximately + * using Monte Carlo sampling. * + * ### Algorithm steps + * 1. Assemble the target mass matrix. + * 2. Assemble the load vector corresponding to the chosen projection method. + * 3. Solve the resulting linear system using PETSc. + * 4. Return the projected nodal field on the target mesh. + * + * @return A vector of nodal values on the target mesh after projection + * (`Omega_h::Reals`). * + * @note The specific inputs required to assemble the load vector depend on the + * projection method used by the solver variant. */ +/** + * @brief Solves the conservative Galerkin projection using exact integration + * over mesh-intersection regions. + * + * @param target_mesh The target Omega_h mesh where the field is projected. + * @param source_mesh The source Omega_h mesh containing the original field. + * @param intersection Precomputed intersection information between source and + * target meshes. + * @param source_values Nodal scalar field values on the source mesh. + * @return Projected nodal values on the target mesh. + */ Omega_h::Reals solveGalerkinProjectionMI( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values); +/** + * @brief Solves the conservative Galerkin projection using Monte Carlo + * integration of the load vector. + * + * @param target_mesh The target Omega_h mesh where the field is projected. + * @param field_values_at_points Source-field values evaluated at the sampled + * physical points in the target mesh. + * @param npoints_each_tri Number of sample points used in each target triangle. + * @param method Sampling method used to define the reference sample pattern. + * @param loadVec_out PETSc load vector associated with the sampled projection + * data, if required by the workflow. + * @param sobol_filename Path to the file containing precomputed Sobol samples + * when Sobol sampling is selected. + * @return Projected nodal values on the target mesh. + */ Omega_h::Reals solveGalerkinProjectionMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out, const std::string sobol_filename = ""); -Omega_h::Reals rhsVectorMI(Omega_h::Mesh& target_mesh, - Omega_h::Mesh& source_mesh, - const IntersectionResults& intersection, - const Omega_h::Reals& source_values); +/** + * @brief Computes the Galerkin right-hand-side vector using mesh-intersection + * method. + * + * This routine computes the global load vector associated with the conservative + * Galerkin projection of a scalar source field onto the target mesh. The load + * vector entries are evaluated by exact integration over the intersection + * regions between the source and target meshes. + * + * The returned vector corresponds to the right-hand side + * \f[ + * b_i = \int_{\Omega} f^s(\mathbf{x}) \, \psi_i(\mathbf{x}) \, + * d\mathbf{\Omega}, + * \f] + * where \f$f^s\f$ is the source field and \f$\psi_i\f$ are the target basis + * functions. + * + * This function computes only the RHS/load vector. It does not assemble or + * solve the full projection system. + * + * @param[in] target_mesh Target Omega_h mesh on which the projected field is + * represented. + * @param[in] source_mesh Source Omega_h mesh containing the original scalar + * field. + * @param[in] intersection Precomputed intersection data between source and + * target meshes. + * @param[in] source_values Nodal scalar field values defined on the source + * mesh. + * + * @return Global Galerkin RHS/load vector as `Omega_h::Reals`. + * + * @note This routine is intended for exact conservative transfer based on mesh + * intersections. + * @note The returned vector may be used independently for diagnostics, + * verification, comparison against approximate methods, or as input to a + * subsequent projection solve. + * + * @see calculateLoadVectorMI + */ +Omega_h::Reals computeRhsVectorMI(Omega_h::Mesh& target_mesh, + Omega_h::Mesh& source_mesh, + const IntersectionResults& intersection, + const Omega_h::Reals& source_values); + +/** + * @brief Computes the Galerkin right-hand-side vector using Monte Carlo + * integration. + * + * This routine computes the global load vector associated with the conservative + * Galerkin projection of a scalar field onto the target mesh using Monte Carlo + * integration on each target element. + * + * The returned vector corresponds to the right-hand side + * \f[ + * b_i = \int_{\Omega} f^s(\mathbf{x}) \, \psi_i(\mathbf{x}) \, d\mathbf{x}, + * \f] + * approximated from sampled source-field values at physical points in the + * target mesh. + * + * The sampling pattern is defined on the reference triangle and may be + * generated either by uniform random sampling or by Sobol-based sampling. + * + * This function computes only the RHS/load vector. It does not assemble or + * solve the full projection system. + * + * @param[in] target_mesh Target Omega_h mesh on which the projected field is + * represented. + * @param[in] field_values_at_points Source-field values evaluated at the sample + * points in the target elements. The expected + * layout is element-wise flattened, with + * total size + * `target_mesh.nelems() * npoints_each_tri`. + * @param[in] npoints_each_tri Number of sample points used in each target + * triangle. + * @param[in] method Sampling method used to define the reference barycentric + * sampling method. + * @param[in] sobol_filename Path to the file containing precomputed Sobol + * barycentric samples when Sobol sampling is used. + * + * @return Global Galerkin RHS/load vector as `Omega_h::Reals`. + * + * @note This routine is intended for approximate conservative transfer using + * Monte Carlo integration. + * @note The returned vector may be used independently for diagnostics, + * verification, comparison against exact mesh-intersection integration, + * or as input to a subsequent projection solve. + * + * @see calculateLoadVectorMC, SamplingMethod + */ +Omega_h::Reals computeRhsVectorMC(Omega_h::Mesh& target_mesh, + const Omega_h::Reals& field_values_at_points, + const int npoints_each_tri, + SamplingMethod method, + const std::string& sobol_filename = ""); } // namespace pcms #endif // PCMS_TRANSFER_CONSERVATIVE_PROJECTION_SOLVER_HPP From d542524b72bfe3addb048da7703302706437d339 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:15:48 -0400 Subject: [PATCH 16/44] add namespace --- src/pcms/transfer/load_vector_integrator_mc.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 1109dac3..6d62842f 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -1,7 +1,9 @@ #include "pcms/transfer/load_vector_integrator.hpp" +namespace pcms +{ // TODO:: create a function that can sample sobol sequences instead of -// generating in python and reading and using +// generating in python and using Kokkos::View read_sobol_barycentric_samples_from_file( std::string file_path) { @@ -237,3 +239,4 @@ Kokkos::View buildLoadVectorMC( return elmLoadVector; } +} // namespace pcms From aea8cb6cc73537fba6e0d35ed7b588bb4de64b10 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:18:31 -0400 Subject: [PATCH 17/44] update filenames in CMakeLists.txt --- src/pcms/transfer/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pcms/transfer/CMakeLists.txt b/src/pcms/transfer/CMakeLists.txt index cbd8889a..01810d29 100644 --- a/src/pcms/transfer/CMakeLists.txt +++ b/src/pcms/transfer/CMakeLists.txt @@ -21,7 +21,8 @@ set(PCMS_FIELD_TRANSFER_HEADERS set(PCMS_FIELD_TRANSFER_SOURCES adj_search.cpp - load_vector_integrator.cpp + load_vector_integrator_mi.cpp + load_vector_integrator_mc.cpp mesh_intersection.cpp mls_interpolation.cpp interpolation_base.cpp) From 15897c13878322d35c620133ccd70e06be3b43da Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:21:44 -0400 Subject: [PATCH 18/44] update solveGalerkinProjection to solveGalerkinProjectionMI --- test/test_mesh_intersection_field_transfer.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index 339e8c93..3dcf4a3b 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -73,8 +73,8 @@ TEST_CASE("mesh intersection linear/constant conservation", Omega_h::parallel_for( source_const.size(), OMEGA_H_LAMBDA(int i) { source_const[i] = c; }); - auto projected = pcms::solveGalerkinProjection(target_mesh, source_mesh, - intersections, source_const); + auto projected = pcms::solveGalerkinProjectionMI( + target_mesh, source_mesh, intersections, source_const); auto projected_h = Omega_h::HostRead(projected); REQUIRE(static_cast(projected.size()) == target_mesh.nverts()); From 5643ec2fd41526deec68a6831585006e03a7a967 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:25:29 -0400 Subject: [PATCH 19/44] update buildLoadVector to buildLoadVectorMI --- test/test_load_vector.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp index 1a38df50..b164b557 100644 --- a/test/test_load_vector.cpp +++ b/test/test_load_vector.cpp @@ -88,7 +88,7 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write values(source_mesh.nverts(), 1.0); auto load_vector = - pcms::buildLoadVector(target_mesh, source_mesh, intersection, values); + pcms::buildLoadVectorMIMI(target_mesh, source_mesh, intersection, values); auto load_vector_host = Kokkos::create_mirror(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); @@ -104,8 +104,8 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") { Omega_h::Write zero_field(source_mesh.nverts(), 0.0); - auto load_vector = - pcms::buildLoadVector(target_mesh, source_mesh, intersection, zero_field); + auto load_vector = pcms::buildLoadVectorMI(target_mesh, source_mesh, + intersection, zero_field); auto load_vector_host = Kokkos::create_mirror_view(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); @@ -124,8 +124,8 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write constant_field(source_mesh.nverts(), 2.0); - auto load_vector = pcms::buildLoadVector(target_mesh, source_mesh, - intersection, constant_field); + auto load_vector = pcms::buildLoadVectorMI(target_mesh, source_mesh, + intersection, constant_field); auto load_vector_host = Kokkos::create_mirror_view(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); From 6356c81924c67c09ac77996b6a2ec79044e233da Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:39:14 -0400 Subject: [PATCH 20/44] add mc test functions --- .../test_monte_carlo_projection_functions.cpp | 330 ++++++++++++++++++ 1 file changed, 330 insertions(+) create mode 100644 test/test_monte_carlo_projection_functions.cpp diff --git a/test/test_monte_carlo_projection_functions.cpp b/test/test_monte_carlo_projection_functions.cpp new file mode 100644 index 00000000..e45ed793 --- /dev/null +++ b/test/test_monte_carlo_projection_functions.cpp @@ -0,0 +1,330 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +constexpr double tol = 1e-12; + +double linear_field(double x, double y) +{ + return 2.0 * x + 3.0 * y + 1.0; +} + +Omega_h::Mesh make_single_triangle_mesh(Omega_h::Library& lib) +{ + // Triangle with vertices: + // v0 = (0,0), v1 = (1,0), v2 = (0,1) + Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 0.0, 1.0}); + + Omega_h::LOs ev2v({0, 1, 2}); + + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + + return mesh; +} + +Omega_h::Mesh make_unit_square_two_tri_mesh(Omega_h::Library& lib) +{ + // Square: + // v0 = (0,0), v1 = (1,0), v2 = (1,1), v3 = (0,1) + // + // Two CCW triangles: + // T0 = (0,1,3) + // T1 = (1,2,3) + Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); + + Omega_h::LOs ev2v({0, 1, 3, 1, 2, 3}); + + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + + return mesh; +} + +Kokkos::View make_barycentric_view( + const std::vector>& vals) +{ + Kokkos::View view("barycentric_coords", vals.size()); + auto host = Kokkos::create_mirror_view(view); + + for (std::size_t i = 0; i < vals.size(); ++i) { + host(i, 0) = vals[i][0]; + host(i, 1) = vals[i][1]; + host(i, 2) = vals[i][2]; + } + + Kokkos::deep_copy(view, host); + return view; +} + +Kokkos::View make_points_view( + const std::vector>& vals) +{ + Kokkos::View view("points", vals.size()); + auto host = Kokkos::create_mirror_view(view); + + for (std::size_t i = 0; i < vals.size(); ++i) { + host(i, 0) = vals[i][0]; + host(i, 1) = vals[i][1]; + } + + Kokkos::deep_copy(view, host); + return view; +} + +Omega_h::Reals make_linear_nodal_field_values(Omega_h::Mesh& mesh) +{ + const int dim = mesh.dim(); + const int nverts = mesh.nverts(); + auto coords = mesh.coords(); + + Omega_h::Write values(nverts, 0.0); + + Omega_h::parallel_for( + nverts, OMEGA_H_LAMBDA(const int i) { + const Omega_h::Real x = coords[i * dim + 0]; + const Omega_h::Real y = coords[i * dim + 1]; + values[i] = linear_field(x, y); + }); + + return Omega_h::read(values); +} + +TEST_CASE("read_sobol_barycentric_samples_from_file reads header and values", + "[mc_projection][sobol_io]") +{ + const std::string filename = "tmp_sobol_barycentric_samples.txt"; + std::ofstream out(filename); + REQUIRE(out.is_open()); + out << "l0 l1 l2\n"; + out << "1.0 0.0 0.0\n"; + out << "0.0 1.0 0.0\n"; + out << "0.2 0.3 0.5\n"; + + auto samples = pcms::read_sobol_barycentric_samples_from_file(filename); + auto host = Kokkos::create_mirror_view(samples); + Kokkos::deep_copy(host, samples); + + REQUIRE(host.extent(0) == 3); + REQUIRE(host.extent(1) == 3); + + CHECK(host(0, 0) == Catch::Approx(1.0).margin(tol)); + CHECK(host(0, 1) == Catch::Approx(0.0).margin(tol)); + CHECK(host(0, 2) == Catch::Approx(0.0).margin(tol)); + + CHECK(host(1, 0) == Catch::Approx(0.0).margin(tol)); + CHECK(host(1, 1) == Catch::Approx(1.0).margin(tol)); + CHECK(host(1, 2) == Catch::Approx(0.0).margin(tol)); + + CHECK(host(2, 0) == Catch::Approx(0.2).margin(tol)); + CHECK(host(2, 1) == Catch::Approx(0.3).margin(tol)); + CHECK(host(2, 2) == Catch::Approx(0.5).margin(tol)); + + std::remove(filename.c_str()); +} + +TEST_CASE( + "generate_uniform_random_barycentric_coords returns valid barycentric rows", + "[mc_projection][sampling]") +{ + const int npoints_each_tri = 280; + + auto bary = + pcms::generate_uniform_random_barycentric_coords(npoints_each_tri); + auto host = Kokkos::create_mirror_view(bary); + Kokkos::deep_copy(host, bary); + + REQUIRE(host.extent(0) == npoints_each_tri); + REQUIRE(host.extent(1) == 3); + + for (int i = 0; i < npoints_each_tri; ++i) { + const double l0 = host(i, 0); + const double l1 = host(i, 1); + const double l2 = host(i, 2); + + CHECK(l0 >= 0.0); + CHECK(l1 >= 0.0); + CHECK(l2 >= 0.0); + + CHECK(l0 <= 1.0); + CHECK(l1 <= 1.0); + CHECK(l2 <= 1.0); + + CHECK(l0 + l1 + l2 == Catch::Approx(1.0).margin(1e-12)); + } +} + +TEST_CASE( + "global_coords_from_ref_barycentric_coords maps reference samples correctly", + "[mc_projection][mapping]") +{ + Omega_h::Library lib; + auto mesh = make_single_triangle_mesh(lib); + + auto ref_barycentric_coords = make_barycentric_view({{1.0, 0.0, 0.0}, + {0.0, 1.0, 0.0}, + {0.0, 0.0, 1.0}, + {0.5, 0.5, 0.0}, + {0.25, 0.25, 0.50}}); + + auto global_points = pcms::global_coords_from_ref_barycentric_coords( + mesh, ref_barycentric_coords); + + auto host_points = Kokkos::create_mirror_view(global_points); + Kokkos::deep_copy(host_points, global_points); + + REQUIRE(host_points.extent(0) == 5); + REQUIRE(host_points.extent(1) == 2); + + // (1,0,0) -> v0 = (0,0) + CHECK(host_points(0, 0) == Catch::Approx(0.0).margin(tol)); + CHECK(host_points(0, 1) == Catch::Approx(0.0).margin(tol)); + + // (0,1,0) -> v1 = (1,0) + CHECK(host_points(1, 0) == Catch::Approx(1.0).margin(tol)); + CHECK(host_points(1, 1) == Catch::Approx(0.0).margin(tol)); + + // (0,0,1) -> v2 = (0,1) + CHECK(host_points(2, 0) == Catch::Approx(0.0).margin(tol)); + CHECK(host_points(2, 1) == Catch::Approx(1.0).margin(tol)); + + // (0.5,0.5,0) -> (0.5,0) + CHECK(host_points(3, 0) == Catch::Approx(0.5).margin(tol)); + CHECK(host_points(3, 1) == Catch::Approx(0.0).margin(tol)); + + // (0.25,0.25,0.5) -> (0.25,0.5) + CHECK(host_points(4, 0) == Catch::Approx(0.25).margin(tol)); + CHECK(host_points(4, 1) == Catch::Approx(0.50).margin(tol)); +} + +TEST_CASE("localize_points_in_mesh identifies inside and outside points", + "[mc_projection][localization]") +{ + Omega_h::Library lib; + auto mesh = make_unit_square_two_tri_mesh(lib); + + auto points = make_points_view({ + {0.2, 0.2}, // inside + {0.8, 0.8}, // inside + {1.2, 0.5} // outside + }); + + auto results = pcms::localize_points_in_mesh(mesh, points); + auto host_results = Kokkos::create_mirror_view(results); + Kokkos::deep_copy(host_results, results); + + REQUIRE(host_results.extent(0) == 3); + + CHECK(host_results(0).tri_id >= 0); + CHECK(host_results(1).tri_id >= 0); + CHECK(host_results(2).tri_id < 0); +} + +TEST_CASE( + "evaluate_field_from_point_localization reproduces linear field exactly", + "[mc_projection][field_eval]") +{ + Omega_h::Library lib; + auto mesh = make_single_triangle_mesh(lib); + + // For f(x,y) = 2x + 3y + 1: + // at (0,0) -> 1 + // at (1,0) -> 3 + // at (0,1) -> 4 + Omega_h::Reals nodal_field_values({1.0, 3.0, 4.0}); + + auto points = make_points_view({{0.2, 0.2}, {0.1, 0.7}, {0.3, 0.1}}); + + auto results = pcms::localize_points_in_mesh(mesh, points); + + auto field_values = pcms::evaluate_field_from_point_localization( + mesh, nodal_field_values, points, results); + + auto host_points = Kokkos::create_mirror_view(points); + Kokkos::deep_copy(host_points, points); + + auto host_values = Omega_h::HostRead(field_values); + + for (int i = 0; i < 3; ++i) { + const double x = host_points(i, 0); + const double y = host_points(i, 1); + const double exact = linear_field(x, y); + + CHECK(host_values[i] == Catch::Approx(exact).margin(1e-12)); + } +} + +TEST_CASE("evaluate_field_from_point_localization works on multi-element mesh", + "[mc_projection][field_eval][multi_elem]") +{ + Omega_h::Library lib; + auto mesh = make_unit_square_two_tri_mesh(lib); + + auto nodal_field_values = make_linear_nodal_field_values(mesh); + + auto points = make_points_view({{0.2, 0.2}, {0.8, 0.8}, {0.75, 0.25}}); + + auto results = pcms::localize_points_in_mesh(mesh, points); + + auto field_values = pcms::evaluate_field_from_point_localization( + mesh, nodal_field_values, points, results); + + auto host_points = Kokkos::create_mirror_view(points); + Kokkos::deep_copy(host_points, points); + + auto host_values = Omega_h::HostRead(field_values); + + for (int i = 0; i < 3; ++i) { + const double x = host_points(i, 0); + const double y = host_points(i, 1); + const double exact = linear_field(x, y); + + CHECK(host_values[i] == Catch::Approx(exact).margin(1e-12)); + } +} + +TEST_CASE("sampling mapping localization and field evaluation pipeline works", + "[mc_projection][integration]") +{ + Omega_h::Library lib; + auto source_mesh = make_unit_square_two_tri_mesh(lib); + auto target_mesh = make_single_triangle_mesh(lib); + + auto ref_barycentric_coords = make_barycentric_view( + {{0.6, 0.2, 0.2}, {0.2, 0.6, 0.2}, {0.2, 0.2, 0.6}, {0.3, 0.3, 0.4}}); + + auto sampled_points = pcms::global_coords_from_ref_barycentric_coords( + target_mesh, ref_barycentric_coords); + + auto results = pcms::localize_points_in_mesh(source_mesh, sampled_points); + + auto nodal_field_values = make_linear_nodal_field_values(source_mesh); + + auto sampled_field_values = pcms::evaluate_field_from_point_localization( + source_mesh, nodal_field_values, sampled_points, results); + + auto host_points = Kokkos::create_mirror_view(sampled_points); + Kokkos::deep_copy(host_points, sampled_points); + + auto host_field_values = + Omega_h::HostRead(sampled_field_values); + + for (int i = 0; i < 4; ++i) { + const double x = host_points(i, 0); + const double y = host_points(i, 1); + const double exact = linear_field(x, y); + + CHECK(host_field_values[i] == Catch::Approx(exact).margin(1e-12)); + } +} From bb20610812533dda78ac0d43b15e06d5a78b40b6 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Wed, 18 Mar 2026 14:39:38 -0400 Subject: [PATCH 21/44] update mc test file in CMakeLists.txt --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c27a6de2..0312c77f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -398,6 +398,7 @@ if(Catch2_FOUND) test_uniform_grid_field.cpp test_interpolation_class.cpp test_mesh_geometry.cpp + test_monte_carlo_projection_functions.cpp test_omega_h_field2_outofbounds.cpp) endif() From 6d6e7e1460875a89bb911bb087bafc5bccd13625 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 10:41:19 -0400 Subject: [PATCH 22/44] add semicolon --- src/pcms/transfer/load_vector_integrator.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index aeb416a0..2c771331 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -25,7 +25,7 @@ #include #include #include -#include +#include #include #include #include @@ -38,6 +38,8 @@ namespace pcms { + +enum class SamplingMethod { SOBOL, UNIFORM }; /** * @brief Computes the load vector for each target element in the conservative * field transfer. @@ -350,7 +352,7 @@ Omega_h::Reals evaluate_field_from_point_localization( Kokkos::View buildLoadVectorMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, const int npoints_each_tri, SamplingMethod method, - const std::string& sobol_filename) + const std::string& sobol_filename); } // namespace pcms From 65a09c98c60be01b521706de9b6748923b722528 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 10:47:14 -0400 Subject: [PATCH 23/44] correct the class name --- src/pcms/transfer/load_vector_integrator.hpp | 4 ++-- src/pcms/transfer/load_vector_integrator_mc.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index 2c771331..a634b0ee 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -288,7 +288,7 @@ Kokkos::View global_coords_from_ref_barycentric_coords( * * @note This routine assumes that `mesh` is two-dimensional. */ -Kokkos::View localize_points_in_mesh( +Kokkos::View localize_points_in_mesh( Omega_h::Mesh& mesh, const Kokkos::View& points); /** @@ -308,7 +308,7 @@ Kokkos::View localize_points_in_mesh( Omega_h::Reals evaluate_field_from_point_localization( Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, const Kokkos::View& points, - const Kokkos::View& results); + const Kokkos::View& results); /** * @brief Compute element-local load vectors using Monte Carlo integration. diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 6d62842f..18772fe8 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -118,7 +118,7 @@ Kokkos::View global_coords_from_ref_barycentric_coords( return global_coords; } -Kokkos::View localize_points_in_mesh( +Kokkos::View localize_points_in_mesh( Omega_h::Mesh& mesh, const Kokkos::View& points) { @@ -143,7 +143,7 @@ Kokkos::View localize_points_in_mesh( Omega_h::Reals evaluate_field_from_point_localization( Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, - const Kokkos::View& results) + const Kokkos::View& results) { const auto& npoints = results.extent(0); From 4f664bca3a20d366bdeb72422ded5557caeff569 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 11:07:00 -0400 Subject: [PATCH 24/44] update variable --- src/pcms/transfer/load_vector_integrator_mc.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 18772fe8..89903d9d 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -9,7 +9,7 @@ Kokkos::View read_sobol_barycentric_samples_from_file( { std::ifstream infile(file_path); if (!infile) { - throw std::runtime_error("Could not open sobol sample file : "; + throw std::runtime_error("Could not open sobol sample file"); } std::vector buffer; @@ -34,13 +34,12 @@ Kokkos::View read_sobol_barycentric_samples_from_file( buffer.push_back(b0); buffer.push_back(b1); buffer.push_back(b2); - ++col; ++nrows; } Kokkos::View device_samples("device_data", nrows); - auto host_samples = Kokkos::create_mirror_view(device_sobol_samples); + auto host_samples = Kokkos::create_mirror_view(device_samples); // Fill host view from buffer for (size_t i = 0; i < nrows; ++i) { From af6e7285e271e09dd1fe4661a7d9cbb532be2bcf Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 11:08:28 -0400 Subject: [PATCH 25/44] correct class name --- src/pcms/transfer/load_vector_integrator_mc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 89903d9d..3c019bce 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -117,7 +117,7 @@ Kokkos::View global_coords_from_ref_barycentric_coords( return global_coords; } -Kokkos::View localize_points_in_mesh( +Kokkos::View localize_points_in_mesh( Omega_h::Mesh& mesh, const Kokkos::View& points) { From a6207e738df612f70a9372f86bf21d8ce902d308 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 11:14:23 -0400 Subject: [PATCH 26/44] change tri_id to element_id --- src/pcms/transfer/load_vector_integrator_mc.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 3c019bce..c30329ea 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -1,4 +1,5 @@ #include "pcms/transfer/load_vector_integrator.hpp" +#include "pcms/utility/assert.h" namespace pcms { @@ -157,7 +158,7 @@ Omega_h::Reals evaluate_field_from_point_localization( npoints, OMEGA_H_LAMBDA(int id) { field_values_at_points[id] = 0.0; - const int tri = results[id].tri_id; + const int tri = results[id].element_id; if (tri < 0) return; const auto el_verts = Omega_h::gather_verts<3>(faces2nodes, tri); From d51b3cd0c06ef1eb71bdf0dae338c8cf4078fb31 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 12:04:45 -0400 Subject: [PATCH 27/44] change tri_id to element_id --- test/test_monte_carlo_projection_functions.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_monte_carlo_projection_functions.cpp b/test/test_monte_carlo_projection_functions.cpp index e45ed793..72acf0ef 100644 --- a/test/test_monte_carlo_projection_functions.cpp +++ b/test/test_monte_carlo_projection_functions.cpp @@ -226,9 +226,9 @@ TEST_CASE("localize_points_in_mesh identifies inside and outside points", REQUIRE(host_results.extent(0) == 3); - CHECK(host_results(0).tri_id >= 0); - CHECK(host_results(1).tri_id >= 0); - CHECK(host_results(2).tri_id < 0); + CHECK(host_results(0).element_id >= 0); + CHECK(host_results(1).element_id >= 0); + CHECK(host_results(2).element_id < 0); } TEST_CASE( From 73d21430021a44eeba61a505040579080acc0633 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 12:05:31 -0400 Subject: [PATCH 28/44] minor changes --- src/pcms/transfer/calculate_load_vector.hpp | 1 + .../transfer/conservative_projection_solver.cpp | 17 +++-------------- .../transfer/conservative_projection_solver.hpp | 1 + src/pcms/transfer/load_vector_integrator_mc.cpp | 6 +++--- 4 files changed, 8 insertions(+), 17 deletions(-) diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp index 0533b72e..7e22f910 100644 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ b/src/pcms/transfer/calculate_load_vector.hpp @@ -13,6 +13,7 @@ #define PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP #include #include +#include #include namespace pcms diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index 35d788e3..1c3b53b6 100644 --- a/src/pcms/transfer/conservative_projection_solver.cpp +++ b/src/pcms/transfer/conservative_projection_solver.cpp @@ -126,16 +126,8 @@ Omega_h::Reals solveGalerkinProjectionMI( Omega_h::Reals solveGalerkinProjectionMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, const int npoints_each_tri, SamplingMethod method, - const std::string sobol_filename); + const std::string sobol_filename) { - if ((PetscInt)source_values.size() != - source_mesh.coords().size() / source_mesh.dim()) { - std::cerr << "ERROR: source_values size (" << source_values.size() - << ") doesn't match expected size (" - << source_mesh.coords().size() / source_mesh.dim() << ")" - << std::endl; - throw std::runtime_error("source_values length mismatch"); - } Mat mass; PetscErrorCode ierr = calculateMassMatrix(target_mesh, &mass); @@ -144,8 +136,6 @@ Omega_h::Reals solveGalerkinProjectionMC( Vec vec; ierr = calculateLoadVectorMC(target_mesh, field_values_at_points, npoints_each_tri, method, &vec, sobol_filename); - ierr = calculateLoadVector(target_mesh, source_mesh, intersection, - source_values, &vec); CHKERRABORT(PETSC_COMM_WORLD, ierr); Vec x = solveLinearSystem(mass, vec); @@ -186,10 +176,9 @@ Omega_h::Reals computeRhsVectorMC(Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, const int npoints_each_tri, SamplingMethod method, - const std::string& sobol_filename = "") + const std::string& sobol_filename) { - { Vec vec; PetscErrorCode ierr; ierr = @@ -202,6 +191,6 @@ Omega_h::Reals computeRhsVectorMC(Omega_h::Mesh& target_mesh, CHKERRABORT(PETSC_COMM_WORLD, ierr); return rhsvector; - } +} } // namespace pcms diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index 2d3feee8..e04c3fa7 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -20,6 +20,7 @@ #include #include +#include namespace pcms { diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index c30329ea..28ea9bf6 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -177,7 +177,7 @@ KOKKOS_INLINE_FUNCTION Omega_h::Vector<3> compute_elm_load_vector_mc( const Kokkos::View& ref_barycentric_coords, const Omega_h::Reals& field_values_at_points, - const Omega_h::Real& elmVolume const int elm) + const Omega_h::Real& elmVolume, const int elm) { const int npoints_each_tri = ref_barycentric_coords.extent(0); @@ -191,7 +191,7 @@ Omega_h::Vector<3> compute_elm_load_vector_mc( sum1 += ref_barycentric_coords(i, 1) * f; sum2 += ref_barycentric_coords(i, 2) * f; } - return elmVolume / npoints_each_tri * {sum0, sum1, sum2}; + return elmVolume / npoints_each_tri * Omega_h::Vector<3>{sum0, sum1, sum2}; } Kokkos::View buildLoadVectorMC( @@ -203,7 +203,7 @@ Kokkos::View buildLoadVectorMC( Kokkos::View ref_barycentric_coords; if (method == SamplingMethod::SOBOL) { if (sobol_filename.empty()) { - throw std::runtime_error("Could not open sobol sample file : "; + throw std::runtime_error("Could not open sobol sample file"); } ref_barycentric_coords = read_sobol_barycentric_samples_from_file(sobol_filename); From 4e626c9aff87d8e5637eee1ce454602d5b8e5542 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 12:33:59 -0400 Subject: [PATCH 29/44] fix the function name --- src/pcms/transfer/calculate_load_vector.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/calculate_load_vector.cpp b/src/pcms/transfer/calculate_load_vector.cpp index 6da01206..3da94631 100644 --- a/src/pcms/transfer/calculate_load_vector.cpp +++ b/src/pcms/transfer/calculate_load_vector.cpp @@ -5,7 +5,7 @@ namespace pcms { -PetscErrorCode calculateLoadVector(Omega_h::Mesh& target_mesh, +PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values, From 5dca1356fd19ea22c29968fef1992a5c93bfdae3 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:40:48 -0400 Subject: [PATCH 30/44] update header --- test/test_load_vector.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp index b164b557..9406b592 100644 --- a/test/test_load_vector.cpp +++ b/test/test_load_vector.cpp @@ -88,7 +88,7 @@ TEST_CASE("Load vector computation on intersected regions", "[load_vector]") Omega_h::Write values(source_mesh.nverts(), 1.0); auto load_vector = - pcms::buildLoadVectorMIMI(target_mesh, source_mesh, intersection, values); + pcms::buildLoadVectorMI(target_mesh, source_mesh, intersection, values); auto load_vector_host = Kokkos::create_mirror(load_vector); Kokkos::deep_copy(load_vector_host, load_vector); From e0915d1eeb8178f75724f8ddb04690bddca06ee8 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:42:37 -0400 Subject: [PATCH 31/44] update header --- test/test_mesh_intersection_field_transfer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index 3dcf4a3b..d220a3ad 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -99,7 +99,7 @@ TEST_CASE("mesh intersection linear/constant conservation", source_linear[i] = x + y; }); - auto projected = pcms::solveGalerkinProjection( + auto projected = pcms::solveGalerkinProjectionMI( target_mesh, source_mesh, intersections, source_linear); auto projected_h = Omega_h::HostRead(projected); const auto tgt_coords_h = From 03103066af5be8fcc4c52032f0233267dac6d794 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:43:00 -0400 Subject: [PATCH 32/44] update function --- test/test_monte_carlo_projection_functions.cpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/test_monte_carlo_projection_functions.cpp b/test/test_monte_carlo_projection_functions.cpp index 72acf0ef..8b745952 100644 --- a/test/test_monte_carlo_projection_functions.cpp +++ b/test/test_monte_carlo_projection_functions.cpp @@ -15,7 +15,7 @@ constexpr double tol = 1e-12; -double linear_field(double x, double y) +OMEGA_H_INLINE double linear_field(double x, double y) { return 2.0 * x + 3.0 * y + 1.0; } @@ -111,6 +111,7 @@ TEST_CASE("read_sobol_barycentric_samples_from_file reads header and values", out << "1.0 0.0 0.0\n"; out << "0.0 1.0 0.0\n"; out << "0.2 0.3 0.5\n"; + out.close(); auto samples = pcms::read_sobol_barycentric_samples_from_file(filename); auto host = Kokkos::create_mirror_view(samples); @@ -249,7 +250,7 @@ TEST_CASE( auto results = pcms::localize_points_in_mesh(mesh, points); auto field_values = pcms::evaluate_field_from_point_localization( - mesh, nodal_field_values, points, results); + mesh, nodal_field_values, results); auto host_points = Kokkos::create_mirror_view(points); Kokkos::deep_copy(host_points, points); @@ -278,7 +279,7 @@ TEST_CASE("evaluate_field_from_point_localization works on multi-element mesh", auto results = pcms::localize_points_in_mesh(mesh, points); auto field_values = pcms::evaluate_field_from_point_localization( - mesh, nodal_field_values, points, results); + mesh, nodal_field_values, results); auto host_points = Kokkos::create_mirror_view(points); Kokkos::deep_copy(host_points, points); @@ -312,7 +313,7 @@ TEST_CASE("sampling mapping localization and field evaluation pipeline works", auto nodal_field_values = make_linear_nodal_field_values(source_mesh); auto sampled_field_values = pcms::evaluate_field_from_point_localization( - source_mesh, nodal_field_values, sampled_points, results); + source_mesh, nodal_field_values, results); auto host_points = Kokkos::create_mirror_view(sampled_points); Kokkos::deep_copy(host_points, sampled_points); From 4556bdc4d6851fc9da73a0dd7aca98e37241114d Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:44:07 -0400 Subject: [PATCH 33/44] minor changes --- src/pcms/transfer/conservative_projection_solver.hpp | 1 + src/pcms/transfer/load_vector_integrator.hpp | 1 - src/pcms/transfer/load_vector_integrator_mc.cpp | 3 ++- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index e04c3fa7..9ec9255d 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -21,6 +21,7 @@ #include #include +#include namespace pcms { diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index a634b0ee..501a3580 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -307,7 +307,6 @@ Kokkos::View localize_points_in_mesh( */ Omega_h::Reals evaluate_field_from_point_localization( Omega_h::Mesh& mesh, const Omega_h::Reals& nodal_field_values, - const Kokkos::View& points, const Kokkos::View& results); /** diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index 28ea9bf6..c96cf115 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -14,7 +14,8 @@ Kokkos::View read_sobol_barycentric_samples_from_file( } std::vector buffer; - size_t nrows = 0, ncols = 0; + size_t nrows = 0; + const size_t ncols = 3; std::string line; // skip header row From 24eccd5823b6332a390b1191d4fc257bedeff36c Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:48:36 -0400 Subject: [PATCH 34/44] add test case for monte carlo field transfer --- test/test_monte_carlo_field_transfer.cpp | 174 +++++++++++++++++++++++ 1 file changed, 174 insertions(+) create mode 100644 test/test_monte_carlo_field_transfer.cpp diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp new file mode 100644 index 00000000..eeb8c8d7 --- /dev/null +++ b/test/test_monte_carlo_field_transfer.cpp @@ -0,0 +1,174 @@ +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) +{ + const auto elem_areas = Omega_h::measure_elements_real(&mesh); + const auto elem_verts = mesh.ask_elem_verts(); + const auto nverts = mesh.nverts(); + + REQUIRE(static_cast(u.size()) == nverts); + + const auto elem_areas_h = Omega_h::HostRead(elem_areas); + const auto elem_verts_h = Omega_h::HostRead(elem_verts); + const auto u_h = Omega_h::HostRead(u); + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + const Omega_h::LO v0 = elem_verts_h[3 * e + 0]; + const Omega_h::LO v1 = elem_verts_h[3 * e + 1]; + const Omega_h::LO v2 = elem_verts_h[3 * e + 2]; + + const double area = elem_areas_h[e]; + const double avg = (u_h[v0] + u_h[v1] + u_h[v2]) / 3.0; + integral += area * avg; + } + return integral; +} + +void write_barycentric_samples_file( + const std::string& filename, + const std::vector>& bary) +{ + std::ofstream out(filename); + REQUIRE(out.is_open()); + out << "l0 l1 l2\n"; + for (const auto& s : bary) { + out << s[0] << " " << s[1] << " " << s[2] << "\n"; + } + out.close(); +} + +Omega_h::Reals evaluate_field_at_barycentric_samples( + Omega_h::Mesh& mesh, + const std::vector>& bary, + const std::function& f) +{ + const int npoints_each_tri = static_cast(bary.size()); + const auto coords_h = Omega_h::HostRead(mesh.coords()); + const auto ev2v_h = Omega_h::HostRead(mesh.ask_elem_verts()); + + Omega_h::Write vals(mesh.nelems() * npoints_each_tri); + + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + const Omega_h::LO v0 = ev2v_h[3 * e + 0]; + const Omega_h::LO v1 = ev2v_h[3 * e + 1]; + const Omega_h::LO v2 = ev2v_h[3 * e + 2]; + + const double x0 = coords_h[2 * v0 + 0]; + const double y0 = coords_h[2 * v0 + 1]; + const double x1 = coords_h[2 * v1 + 0]; + const double y1 = coords_h[2 * v1 + 1]; + const double x2 = coords_h[2 * v2 + 0]; + const double y2 = coords_h[2 * v2 + 1]; + + const int base = e * npoints_each_tri; + for (int i = 0; i < npoints_each_tri; ++i) { + const double l0 = bary[i][0]; + const double l1 = bary[i][1]; + const double l2 = bary[i][2]; + + const double x = l0 * x0 + l1 * x1 + l2 * x2; + const double y = l0 * y0 + l1 * y1 + l2 * y2; + + vals[base + i] = f(x, y); + } + } + + return Omega_h::Reals(vals); +} + + +TEST_CASE("monte carlo projection preserves constant and linear fields", + "[transfer][mc_projection]") +{ + Omega_h::Library lib; + + Omega_h::Reals coords({ + 0.0, 0.0, + 1.0, 0.0, + 1.0, 1.0, + 0.0, 1.0 + }); + + Omega_h::LOs ev2v_target({0, 1, 3, 1, 2, 3}); + Omega_h::Mesh target_mesh(&lib); + Omega_h::build_from_elems_and_coords( + &target_mesh, OMEGA_H_SIMPLEX, 2, ev2v_target, coords); + + // Deterministic barycentric samples for the reference triangle. + // Replace this with many more Sobol samples in real application + std::vector> bary; + bary.reserve(8); + bary.push_back({1.0/3.0, 1.0/3.0, 1.0/3.0}); + bary.push_back({0.6, 0.2, 0.2}); + bary.push_back({0.2, 0.6, 0.2}); + bary.push_back({0.2, 0.2, 0.6}); + bary.push_back({0.5, 0.4, 0.1}); + bary.push_back({0.5, 0.1, 0.4}); + bary.push_back({0.1, 0.5, 0.4}); + bary.push_back({0.1, 0.4, 0.5}); + + const std::string sobol_filename = "tmp_mc_barycentric_samples.txt"; + write_barycentric_samples_file(sobol_filename, bary); + + const int npoints_each_tri = static_cast(bary.size()); + + SECTION("constant field is preserved and conserved") + { + const double c = 2.0; + + auto field_values_at_points = evaluate_field_at_barycentric_samples( + target_mesh, bary, [c](double, double) { return c; }); + + Vec loadVec = nullptr; + auto projected = pcms::solveGalerkinProjectionMC( + target_mesh, field_values_at_points, npoints_each_tri, + SamplingMethod::SOBOL, &loadVec, sobol_filename); + + auto projected_h = Omega_h::HostRead(projected); + + REQUIRE(static_cast(projected.size()) == target_mesh.nverts()); + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + REQUIRE(projected_h[i] == Catch::Approx(c).margin(1e-10)); + } + + const double exact_integral = c * 1.0; // unit square area = 1 + const double projected_integral = + integrate_linear_field(target_mesh, projected); + REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(1e-10)); + } + + SECTION("linear field is reproduced approximately") + { + auto field_values_at_points = evaluate_field_at_barycentric_samples( + target_mesh, bary, [](double x, double y) { return x + y; }); + + Vec loadVec = nullptr; + auto projected = pcms::solveGalerkinProjectionMC( + target_mesh, field_values_at_points, npoints_each_tri, + SamplingMethod::SOBOL, &loadVec, sobol_filename); + + auto projected_h = Omega_h::HostRead(projected); + auto tgt_coords_h = Omega_h::HostRead(target_mesh.coords()); + + for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { + const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + REQUIRE(projected_h[i] == Catch::Approx(expected).margin(1e-6)); + } + + const double exact_integral = 1.0; // integral of x+y over unit square + const double projected_integral = + integrate_linear_field(target_mesh, projected); + REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(1e-6)); + } +} From 083144ccc3ec22a06d842e9fe5a5abcde0cf9ba0 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:50:37 -0400 Subject: [PATCH 35/44] add filename in CMakeLists.txt --- test/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 0312c77f..7063bd73 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -409,7 +409,8 @@ if(Catch2_FOUND) if(PCMS_ENABLE_PETSC) list(APPEND PCMS_UNIT_TEST_SOURCES - test_mesh_intersection_field_transfer.cpp) + test_mesh_intersection_field_transfer.cpp + test_monte_carlo_field_transfer.cpp) endif() add_executable(unit_tests ${PCMS_UNIT_TEST_SOURCES}) From 30785f12c8f919184128beaa15a39f69cdc3f876 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:51:48 -0400 Subject: [PATCH 36/44] update header --- test/test_monte_carlo_field_transfer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp index eeb8c8d7..002a9336 100644 --- a/test/test_monte_carlo_field_transfer.cpp +++ b/test/test_monte_carlo_field_transfer.cpp @@ -6,7 +6,7 @@ #include #include -#include +#include #include From 82f2fa84e6f1a31cab228e51cd36e1f2c0567583 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:55:09 -0400 Subject: [PATCH 37/44] update solveGalerkinProjectionMC parameters --- src/pcms/transfer/conservative_projection_solver.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index 9ec9255d..d13cdd12 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -90,7 +90,7 @@ Omega_h::Reals solveGalerkinProjectionMI( */ Omega_h::Reals solveGalerkinProjectionMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, - const int npoints_each_tri, SamplingMethod method, Vec* loadVec_out, + const int npoints_each_tri, SamplingMethod method, const std::string sobol_filename = ""); /** From 7a9785d1ee864a57e014dd2ecf7dedfcb0da0f88 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 13:58:18 -0400 Subject: [PATCH 38/44] update args in a solveGalerkinProjectionMC function --- test/test_monte_carlo_field_transfer.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp index 002a9336..bacbb76d 100644 --- a/test/test_monte_carlo_field_transfer.cpp +++ b/test/test_monte_carlo_field_transfer.cpp @@ -130,10 +130,9 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", auto field_values_at_points = evaluate_field_at_barycentric_samples( target_mesh, bary, [c](double, double) { return c; }); - Vec loadVec = nullptr; auto projected = pcms::solveGalerkinProjectionMC( target_mesh, field_values_at_points, npoints_each_tri, - SamplingMethod::SOBOL, &loadVec, sobol_filename); + pcms::SamplingMethod::SOBOL,sobol_filename); auto projected_h = Omega_h::HostRead(projected); @@ -153,10 +152,9 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", auto field_values_at_points = evaluate_field_at_barycentric_samples( target_mesh, bary, [](double x, double y) { return x + y; }); - Vec loadVec = nullptr; auto projected = pcms::solveGalerkinProjectionMC( target_mesh, field_values_at_points, npoints_each_tri, - SamplingMethod::SOBOL, &loadVec, sobol_filename); + pcms::SamplingMethod::SOBOL, sobol_filename); auto projected_h = Omega_h::HostRead(projected); auto tgt_coords_h = Omega_h::HostRead(target_mesh.coords()); From 5a3d44a1d2fb0b48a55de3f6da12f2e61d53943c Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 14:05:51 -0400 Subject: [PATCH 39/44] create a host view instead of device --- test/test_monte_carlo_field_transfer.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp index bacbb76d..6e273ce5 100644 --- a/test/test_monte_carlo_field_transfer.cpp +++ b/test/test_monte_carlo_field_transfer.cpp @@ -57,7 +57,7 @@ Omega_h::Reals evaluate_field_at_barycentric_samples( const auto coords_h = Omega_h::HostRead(mesh.coords()); const auto ev2v_h = Omega_h::HostRead(mesh.ask_elem_verts()); - Omega_h::Write vals(mesh.nelems() * npoints_each_tri); + Omega_h::HostWrite vals_h(mesh.nelems() * npoints_each_tri); for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { const Omega_h::LO v0 = ev2v_h[3 * e + 0]; @@ -80,11 +80,11 @@ Omega_h::Reals evaluate_field_at_barycentric_samples( const double x = l0 * x0 + l1 * x1 + l2 * x2; const double y = l0 * y0 + l1 * y1 + l2 * y2; - vals[base + i] = f(x, y); + vals_h[base + i] = f(x, y); } } - return Omega_h::Reals(vals); + return Omega_h::Reals(vals_h); } From 80f8aba3502b21a86dac6221d9c4a308d75907e9 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 14:47:20 -0400 Subject: [PATCH 40/44] update the API docs for MC projection --- .../conservative_projection_solver.hpp | 77 +++++++++++++++---- 1 file changed, 64 insertions(+), 13 deletions(-) diff --git a/src/pcms/transfer/conservative_projection_solver.hpp b/src/pcms/transfer/conservative_projection_solver.hpp index d13cdd12..0c0da649 100644 --- a/src/pcms/transfer/conservative_projection_solver.hpp +++ b/src/pcms/transfer/conservative_projection_solver.hpp @@ -5,7 +5,7 @@ * * Provides the main interface to perform Galerkin projection of scalar fields * from a source mesh to a target mesh using conservative transfer using a - * supermesh generated from mesh intersections. + * supermesh generated from mesh intersections or by using stochastic method. * * The solver computes the right-hand side (load vector), assembles the mass * matrix, and solves the resulting linear system to obtain projected nodal @@ -72,21 +72,72 @@ namespace pcms Omega_h::Reals solveGalerkinProjectionMI( Omega_h::Mesh& target_mesh, Omega_h::Mesh& source_mesh, const IntersectionResults& intersection, const Omega_h::Reals& source_values); - /** - * @brief Solves the conservative Galerkin projection using Monte Carlo - * integration of the load vector. + * @brief Solve the conservative Galerkin projection using Monte Carlo + * integration of the load vector on the target mesh. + * + * This routine computes the projected nodal field on the target mesh by + * approximating the Galerkin load vector with Monte Carlo integration over each + * target element. The sampled source-field values must already be provided at + * the physical sample points in each target element. + * + * The expected workflow is: + * 1. Generate or read barycentric sample coordinates on the reference triangle + * using either `generate_uniform_random_barycentric_coords()` or + * `read_sobol_barycentric_samples_from_file()`. + * 2. Map those reference barycentric coordinates to physical sample points in + * all target elements using + * `global_coords_from_ref_barycentric_coords()`. + * 3. Obtain source-field values at those physical sample points. + * - If a source mesh and source nodal field are available, this is done by + * localizing the sample points in the source mesh using + * `localize_points_in_mesh()` and then evaluating the source field with + * `evaluate_field_from_point_localization()`. + * - If no source mesh is available, this routine may still be used provided + * that the source-field values can be queried directly at physical points, + * for example from a black-box field-evaluation interface. + * 4. Pass the resulting sampled values as `field_values_at_points` to this + * routine. + * + * The array `field_values_at_points` is assumed to be stored element-by-element, + * with `npoints_each_tri` consecutive values for each target triangle. That is, + * for target element `e` and local sample index `i`, the sampled field value is + * stored at: + * @code + * field_values_at_points[e * npoints_each_tri + i] + * @endcode + * + * When `method == SamplingMethod::SOBOL`, the routine uses the barycentric + * samples read from `sobol_filename`. Otherwise, it generates uniform random + * barycentric samples internally according to `method`. + * + * @param[in] target_mesh Target 2D triangular mesh on which the projected field + * is represented. + * @param[in] field_values_at_points Source-field values evaluated at the + * sampled physical points in the target mesh, + * stored element-by-element. + * @param[in] npoints_each_tri Number of Monte Carlo sample points used in each + * target triangle. + * @param[in] method Sampling method used to define the reference-triangle + * barycentric sample pattern. + * @param[in] sobol_filename Path to a text file containing precomputed Sobol + * barycentric samples when Sobol sampling is selected. * - * @param target_mesh The target Omega_h mesh where the field is projected. - * @param field_values_at_points Source-field values evaluated at the sampled - * physical points in the target mesh. - * @param npoints_each_tri Number of sample points used in each target triangle. - * @param method Sampling method used to define the reference sample pattern. - * @param loadVec_out PETSc load vector associated with the sampled projection - * data, if required by the workflow. - * @param sobol_filename Path to the file containing precomputed Sobol samples - * when Sobol sampling is selected. * @return Projected nodal values on the target mesh. + * + * @note This routine does not require the source mesh directly. It only + * requires source-field values evaluated at the sampled physical points. + * This makes it applicable both to standard mesh-based transfer and to + * black-box settings where field values can be obtained through point + * queries. + * @note This routine assumes a 2D triangular target mesh. + * + * @see read_sobol_barycentric_samples_from_file + * @see generate_uniform_random_barycentric_coords + * @see global_coords_from_ref_barycentric_coords + * @see localize_points_in_mesh + * @see evaluate_field_from_point_localization + * @see computeRhsVectorMC */ Omega_h::Reals solveGalerkinProjectionMC( Omega_h::Mesh& target_mesh, const Omega_h::Reals& field_values_at_points, From e0c310d1e0e84f30880197a7800854335d5f48fa Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 17:19:03 -0400 Subject: [PATCH 41/44] add sobol samples --- test/test_monte_carlo_field_transfer.cpp | 158 ++++++++++++++++++----- 1 file changed, 126 insertions(+), 32 deletions(-) diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp index 6e273ce5..c2354dcb 100644 --- a/test/test_monte_carlo_field_transfer.cpp +++ b/test/test_monte_carlo_field_transfer.cpp @@ -10,6 +10,126 @@ #include +Kokkos::View make_sobol_barycentric_samples() +{ + static constexpr double data[100][3] = { + {0.77247435, 0.02761234, 0.19991331}, + {0.00066204, 0.93645417, 0.06288378}, + {0.24669859, 0.20114795, 0.55215346}, + {0.37601817, 0.44040299, 0.18357885}, + {0.45721860, 0.21639685, 0.32638456}, + {0.14873985, 0.50203922, 0.34922092}, + {0.08265261, 0.21955372, 0.69779366}, + {0.61711447, 0.30763833, 0.07524720}, + {0.55756878, 0.14595971, 0.29647151}, + {0.13181260, 0.55790360, 0.31028380}, + {0.17533046, 0.04910693, 0.77556261}, + {0.39746965, 0.60204336, 0.00048699}, + {0.32622724, 0.11952704, 0.55424572}, + {0.27416548, 0.62836549, 0.09746903}, + {0.04762209, 0.43946455, 0.51291336}, + {0.67610228, 0.17059590, 0.15330182}, + {0.70636627, 0.13868367, 0.15495006}, + {0.05423990, 0.50286400, 0.44289610}, + {0.26554739, 0.10223092, 0.63222169}, + {0.31251591, 0.56755339, 0.11993070}, + {0.44019781, 0.00429113, 0.55551107}, + {0.20250786, 0.75080234, 0.04668980}, + {0.10680470, 0.31658569, 0.57660961}, + {0.50457450, 0.32955129, 0.16587422}, + {0.59831281, 0.08088912, 0.32079807}, + {0.07370338, 0.70727787, 0.21901875}, + {0.15846319, 0.34444054, 0.49709627}, + {0.47101819, 0.31449160, 0.21449022}, + {0.34034786, 0.19216696, 0.46748519}, + {0.21550916, 0.57121237, 0.21327847}, + {0.02496562, 0.06798091, 0.90705347}, + {0.92252284, 0.06814746, 0.00932970}, + {0.85768371, 0.03637656, 0.10593973}, + {0.01619319, 0.68252604, 0.30128077}, + {0.22649636, 0.07900530, 0.69449834}, + {0.35119355, 0.59511537, 0.05369109}, + {0.48835031, 0.11281048, 0.39883921}, + {0.16631331, 0.65300768, 0.18067901}, + {0.06666197, 0.36178388, 0.57155415}, + {0.57652986, 0.24449693, 0.17897321}, + {0.52216375, 0.01874412, 0.45909213}, + {0.11478599, 0.86792052, 0.01729349}, + {0.19362412, 0.25593162, 0.55044426}, + {0.42505536, 0.36317736, 0.21176728}, + {0.30305820, 0.31329553, 0.38364627}, + {0.25341301, 0.38460776, 0.36197923}, + {0.06382536, 0.14763042, 0.78854422}, + {0.72968702, 0.22878500, 0.04152798}, + {0.65489003, 0.05202033, 0.29308965}, + {0.03867714, 0.80507099, 0.15625187}, + {0.28610086, 0.35098517, 0.36291398}, + {0.33688728, 0.36564553, 0.29746719}, + {0.41257474, 0.21920914, 0.36821612}, + {0.18330823, 0.55984932, 0.25684245}, + {0.12427171, 0.01628781, 0.85944048}, + {0.53763460, 0.44108409, 0.02128131}, + {0.64025214, 0.15451016, 0.20523770}, + {0.09045185, 0.55789331, 0.35165485}, + {0.14044321, 0.18374592, 0.67581086}, + {0.44167163, 0.43248795, 0.12584041}, + {0.36556721, 0.05186860, 0.58256418}, + {0.23481881, 0.68177609, 0.08340511}, + {0.00969458, 0.30812735, 0.68217807}, + {0.80334176, 0.14697275, 0.04968549}, + {0.79815936, 0.07689627, 0.12494437}, + {0.01472633, 0.55979307, 0.42548060}, + {0.23349852, 0.17606379, 0.59043770}, + {0.37341621, 0.49494930, 0.13163449}, + {0.45414551, 0.05244857, 0.49340592}, + {0.13694848, 0.78356864, 0.07948288}, + {0.09804517, 0.23897485, 0.66297999}, + {0.63200480, 0.25752600, 0.11046920}, + {0.54417024, 0.07635918, 0.37947058}, + {0.11645350, 0.75270655, 0.13083994}, + {0.18700282, 0.36095309, 0.45204410}, + {0.40096509, 0.30271633, 0.29631858}, + {0.32946503, 0.21998856, 0.45054641}, + {0.28751929, 0.45375342, 0.25872729}, + {0.03352118, 0.03288756, 0.93359126}, + {0.65789378, 0.33176733, 0.01033890}, + {0.74042350, 0.00658366, 0.25299284}, + {0.05660849, 0.90854533, 0.03484618}, + {0.25725536, 0.27042428, 0.47232036}, + {0.29337009, 0.47961516, 0.22701475}, + {0.41676634, 0.29019891, 0.29303475}, + {0.19477005, 0.45163734, 0.35359261}, + {0.10936446, 0.12581988, 0.76481566}, + {0.52412918, 0.39570142, 0.08016940}, + {0.57432353, 0.12821589, 0.29746058}, + {0.07183352, 0.68862912, 0.23953735}, + {0.16520587, 0.07277389, 0.76202024}, + {0.49782046, 0.45243601, 0.04974353}, + {0.36175953, 0.12969465, 0.50854582}, + {0.22280546, 0.59788723, 0.17930731}, + {0.02311096, 0.42476114, 0.55212791}, + {0.83893443, 0.10049592, 0.06056965}, + {0.91142950, 0.02213425, 0.06643625}, + {0.02989667, 0.78460762, 0.18549570}, + {0.21430992, 0.30905766, 0.47663242}, + {0.34768998, 0.37781691, 0.27449311}, + {0.48400348, 0.14291404, 0.37308247}, + {0.15508466, 0.60091608, 0.24399926} + }; + + Kokkos::View samples("sobol_bary", 100); + auto host = Kokkos::create_mirror_view(samples); + + for (int i = 0; i < 100; ++i) { + host(i, 0) = data[i][0]; + host(i, 1) = data[i][1]; + host(i, 2) = data[i][2]; + } + + Kokkos::deep_copy(samples, host); + return samples; +} + double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) { const auto elem_areas = Omega_h::measure_elements_real(&mesh); @@ -35,25 +155,15 @@ double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) return integral; } -void write_barycentric_samples_file( - const std::string& filename, - const std::vector>& bary) -{ - std::ofstream out(filename); - REQUIRE(out.is_open()); - out << "l0 l1 l2\n"; - for (const auto& s : bary) { - out << s[0] << " " << s[1] << " " << s[2] << "\n"; - } - out.close(); -} Omega_h::Reals evaluate_field_at_barycentric_samples( Omega_h::Mesh& mesh, - const std::vector>& bary, + const Kokkos::View& bary, const std::function& f) { - const int npoints_each_tri = static_cast(bary.size()); + const int npoints_each_tri = bary.extent(0); + const auto bary_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), bary); const auto coords_h = Omega_h::HostRead(mesh.coords()); const auto ev2v_h = Omega_h::HostRead(mesh.ask_elem_verts()); @@ -105,24 +215,8 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", Omega_h::build_from_elems_and_coords( &target_mesh, OMEGA_H_SIMPLEX, 2, ev2v_target, coords); - // Deterministic barycentric samples for the reference triangle. - // Replace this with many more Sobol samples in real application - std::vector> bary; - bary.reserve(8); - bary.push_back({1.0/3.0, 1.0/3.0, 1.0/3.0}); - bary.push_back({0.6, 0.2, 0.2}); - bary.push_back({0.2, 0.6, 0.2}); - bary.push_back({0.2, 0.2, 0.6}); - bary.push_back({0.5, 0.4, 0.1}); - bary.push_back({0.5, 0.1, 0.4}); - bary.push_back({0.1, 0.5, 0.4}); - bary.push_back({0.1, 0.4, 0.5}); - - const std::string sobol_filename = "tmp_mc_barycentric_samples.txt"; - write_barycentric_samples_file(sobol_filename, bary); - - const int npoints_each_tri = static_cast(bary.size()); - + auto bary = make_sobol_barycentric_samples(); + const int npoints_each_tri = bary.extent(0); SECTION("constant field is preserved and conserved") { const double c = 2.0; From 86285ff3bd9a419c342d4ea3c10ed4c74a7fcb18 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 23:30:52 -0400 Subject: [PATCH 42/44] added more sobol barycentric points --- test/test_monte_carlo_field_transfer.cpp | 355 +++++++++++++++++++++-- 1 file changed, 337 insertions(+), 18 deletions(-) diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp index c2354dcb..a7721756 100644 --- a/test/test_monte_carlo_field_transfer.cpp +++ b/test/test_monte_carlo_field_transfer.cpp @@ -12,7 +12,7 @@ Kokkos::View make_sobol_barycentric_samples() { - static constexpr double data[100][3] = { + static constexpr double data[][3] = { {0.77247435, 0.02761234, 0.19991331}, {0.00066204, 0.93645417, 0.06288378}, {0.24669859, 0.20114795, 0.55215346}, @@ -114,13 +114,313 @@ Kokkos::View make_sobol_barycentric_samples() {0.21430992, 0.30905766, 0.47663242}, {0.34768998, 0.37781691, 0.27449311}, {0.48400348, 0.14291404, 0.37308247}, - {0.15508466, 0.60091608, 0.24399926} + {0.15508466, 0.60091608, 0.24399926}, + {0.08106798, 0.10653392, 0.81239810}, + {0.59125760, 0.37895973, 0.02978268}, + {0.51037808, 0.22302710, 0.26659482}, + {0.09923018, 0.46590616, 0.43486366}, + {0.20608878, 0.14835960, 0.64555162}, + {0.42820563, 0.49839775, 0.07339662}, + {0.30554488, 0.03664965, 0.65780547}, + {0.26683109, 0.72591799, 0.00725092}, + {0.04918309, 0.32232089, 0.62849602}, + {0.70951817, 0.18863056, 0.10185128}, + {0.68554941, 0.10809439, 0.20635620}, + {0.04041478, 0.63331619, 0.32626903}, + {0.27830721, 0.00927806, 0.71241473}, + {0.31610340, 0.65117105, 0.03272556}, + {0.38928488, 0.07901356, 0.53170157}, + {0.17657948, 0.67520532, 0.14821519}, + {0.12612715, 0.41757651, 0.45629634}, + {0.55987203, 0.23834031, 0.20178765}, + {0.61447032, 0.02921133, 0.35631835}, + {0.08806645, 0.81064488, 0.10128867}, + {0.14753106, 0.24037411, 0.61209483}, + {0.46644535, 0.38526758, 0.14828707}, + {0.38713171, 0.25491614, 0.35795215}, + {0.24272908, 0.45719217, 0.30007875}, + {0.00758013, 0.19071151, 0.80170835}, + {0.75958884, 0.18196864, 0.05844252}, + {0.75642498, 0.07641374, 0.16716128}, + {0.00420685, 0.62520745, 0.37058570}, + {0.24168664, 0.03264327, 0.72567009}, + {0.38172359, 0.60846830, 0.00980811}, + {0.46561185, 0.08624877, 0.44813937}, + {0.14316113, 0.72845126, 0.12838760}, + {0.08758533, 0.40669566, 0.50571901}, + {0.60488600, 0.20208329, 0.19303072}, + {0.56507343, 0.04615150, 0.38877507}, + {0.12830905, 0.80315478, 0.06853617}, + {0.17934153, 0.20637694, 0.61428154}, + {0.39241969, 0.41909893, 0.18848138}, + {0.32026909, 0.26075533, 0.41897559}, + {0.28008103, 0.41270818, 0.30721079}, + {0.04340562, 0.21483531, 0.74175907}, + {0.68956069, 0.24440279, 0.06603652}, + {0.71782729, 0.06140852, 0.22076418}, + {0.05101451, 0.73865943, 0.21032606}, + {0.27004345, 0.31076297, 0.41919357}, + {0.30809653, 0.42173405, 0.27016942}, + {0.43305433, 0.17420283, 0.39274285}, + {0.20787862, 0.58903325, 0.20308814}, + {0.10230183, 0.07673665, 0.82096152}, + {0.51328268, 0.43555765, 0.05115967}, + {0.58899094, 0.20040958, 0.21059947}, + {0.07753719, 0.50501920, 0.41744360}, + {0.15398298, 0.13086653, 0.71515049}, + {0.47774603, 0.43947285, 0.08278112}, + {0.34722477, 0.01479171, 0.63798352}, + {0.20946864, 0.75724473, 0.03328664}, + {0.02961279, 0.35832657, 0.61206064}, + {0.87537282, 0.08492542, 0.03970177}, + {0.83027487, 0.00941023, 0.16031491}, + {0.02030022, 0.97498961, 0.00471017}, + {0.22099723, 0.26012123, 0.51888153}, + {0.35742635, 0.41555137, 0.22702228}, + {0.49434261, 0.23541477, 0.27024262}, + {0.16237386, 0.44450917, 0.39311698}, + {0.06994161, 0.16129971, 0.76875868}, + {0.56880831, 0.37150678, 0.05968491}, + {0.53224122, 0.12667801, 0.34108077}, + {0.10969413, 0.63166290, 0.25864297}, + {0.19954506, 0.09417638, 0.70627856}, + {0.41725926, 0.54386483, 0.03887591}, + {0.29800254, 0.16533382, 0.53666363}, + {0.25849374, 0.59305467, 0.14845158}, + {0.06004746, 0.37818531, 0.56176723}, + {0.74408352, 0.15190914, 0.10400734}, + {0.66896052, 0.13678923, 0.19425025}, + {0.03399412, 0.57824306, 0.38776282}, + {0.29280627, 0.13953571, 0.56765802}, + {0.33010489, 0.50887768, 0.16101743}, + {0.40657264, 0.03903293, 0.55439442}, + {0.18795978, 0.71097895, 0.10106128}, + {0.12024300, 0.25983911, 0.61991789}, + {0.54588448, 0.33244644, 0.12166908}, + {0.62857616, 0.05031450, 0.32110935}, + {0.09488745, 0.74370051, 0.16141204}, + {0.13548263, 0.41166905, 0.45284832}, + {0.44894084, 0.29505688, 0.25600228}, + {0.37040176, 0.22552560, 0.40407264}, + {0.23051849, 0.51481316, 0.25466836}, + {0.01277710, 0.00376135, 0.98346155}, + {0.78725897, 0.19948623, 0.01325480}, + {0.82195808, 0.03258204, 0.14545988}, + {0.00990231, 0.85926894, 0.13082875}, + {0.23940874, 0.34932966, 0.41126160}, + {0.36586971, 0.33044399, 0.30368631}, + {0.44726339, 0.18952789, 0.36320872}, + {0.14105568, 0.56099359, 0.29795073}, + {0.09385519, 0.04442906, 0.86171575}, + {0.64176458, 0.35329500, 0.00494043}, + {0.53540888, 0.18463415, 0.27995698}, + {0.12130680, 0.51241518, 0.36627803}, + {0.18206767, 0.20116483, 0.61676751}, + {0.40813395, 0.47629744, 0.11556862}, + {0.33439267, 0.07444933, 0.59115800}, + {0.28325405, 0.66170341, 0.05504254}, + {0.03694131, 0.27051844, 0.69254025}, + {0.64907563, 0.25101550, 0.09990886}, + {0.72535642, 0.07845303, 0.19619055}, + {0.06114412, 0.68153434, 0.25732154}, + {0.25184627, 0.05381980, 0.69433393}, + {0.29944780, 0.62015638, 0.08039582}, + {0.42245179, 0.10883681, 0.46871140}, + {0.19098435, 0.60927365, 0.19974200}, + {0.11306857, 0.37229210, 0.51463933}, + {0.51776683, 0.29293065, 0.18930251}, + {0.58506052, 0.00368122, 0.41125826}, + {0.06672258, 0.88474617, 0.04853125}, + {0.17063862, 0.28852920, 0.54083218}, + {0.48842080, 0.33968153, 0.17189767}, + {0.35585174, 0.31039706, 0.33375120}, + {0.22737862, 0.42156721, 0.35105417}, + {0.01924874, 0.12284985, 0.85790141}, + {0.86263256, 0.11208975, 0.02527769}, + {0.95817945, 0.01840331, 0.02341725}, + {0.02650898, 0.48796507, 0.48552595}, + {0.21820855, 0.13415044, 0.64764100}, + {0.34265084, 0.56259210, 0.09475706}, + {0.47578167, 0.01991735, 0.50430098}, + {0.15987053, 0.81813824, 0.02199123}, + {0.07643449, 0.29925505, 0.62431045}, + {0.60122538, 0.25238860, 0.14638602}, + {0.50320834, 0.11597001, 0.38082164}, + {0.10343688, 0.71171168, 0.18485144}, + {0.20163755, 0.30104856, 0.49731389}, + {0.43486931, 0.31896405, 0.24616665}, + {0.31243986, 0.17957018, 0.50798996}, + {0.26069115, 0.51451830, 0.22479055}, + {0.05419924, 0.09454110, 0.85125966}, + {0.69437709, 0.27862917, 0.02699374}, + {0.67447426, 0.02971719, 0.29580855}, + {0.04437040, 0.86481066, 0.09081893}, + {0.27342374, 0.21592583, 0.51065043}, + {0.32165352, 0.50049288, 0.17785360}, + {0.39712442, 0.25964003, 0.34323555}, + {0.17111712, 0.51386702, 0.31501586}, + {0.13159658, 0.17992854, 0.68847488}, + {0.54971857, 0.34826318, 0.10201825}, + {0.62247930, 0.13601247, 0.24150822}, + {0.08445530, 0.61801344, 0.29753125}, + {0.15111408, 0.02468658, 0.82419933}, + {0.46031064, 0.52180822, 0.01788114}, + {0.38020573, 0.08993594, 0.52985833}, + {0.24806226, 0.62803177, 0.12390598}, + {0.00328046, 0.49226494, 0.50445459}, + {0.77699014, 0.12422558, 0.09878428}, + {0.78037086, 0.10871693, 0.11091221}, + {0.00256173, 0.55443851, 0.44299976}, + {0.24905797, 0.11131422, 0.63962781}, + {0.37905057, 0.51672443, 0.10422500}, + {0.46073568, 0.01611880, 0.52314552}, + {0.15080906, 0.82042003, 0.02877092}, + {0.08470579, 0.33223038, 0.58306383}, + {0.62179392, 0.25425402, 0.12395206}, + {0.55356305, 0.09393561, 0.35250133}, + {0.12964776, 0.67030717, 0.20004507}, + {0.17319929, 0.35717574, 0.46962497}, + {0.39432064, 0.37472435, 0.23095501}, + {0.32343227, 0.20286983, 0.47369790}, + {0.27172408, 0.53543229, 0.19284363}, + {0.04563221, 0.08775343, 0.86661436}, + {0.67069807, 0.29776921, 0.03153272}, + {0.69996875, 0.02960151, 0.27042974}, + {0.05237750, 0.86509673, 0.08252577}, + {0.26298510, 0.19013757, 0.54687733}, + {0.30993604, 0.48248765, 0.20757630}, + {0.43705964, 0.21192089, 0.35101946}, + {0.20012840, 0.45200540, 0.34786619}, + {0.10481590, 0.20660818, 0.68857592}, + {0.50078670, 0.39757326, 0.10164003}, + {0.60302891, 0.12737390, 0.26959718}, + {0.07562574, 0.58790810, 0.33646616}, + {0.16072509, 0.03092437, 0.80835054}, + {0.47435812, 0.51251802, 0.01312386}, + {0.34304639, 0.11093623, 0.54601738}, + {0.21791652, 0.67152784, 0.11055563}, + {0.02677604, 0.42745403, 0.54576994}, + {0.95303562, 0.02357760, 0.02338678}, + {0.87170440, 0.01642484, 0.11187077}, + {0.01798891, 0.79865357, 0.18335752}, + {0.22893986, 0.37212694, 0.39893320}, + {0.35393519, 0.35199286, 0.29407195}, + {0.49180135, 0.17835918, 0.32983947}, + {0.16859828, 0.54943574, 0.28196599}, + {0.06857135, 0.00935645, 0.92207219}, + {0.58099726, 0.39674815, 0.02225459}, + {0.51824256, 0.20256397, 0.27919347}, + {0.11277663, 0.53835958, 0.34886380}, + {0.19126783, 0.15451537, 0.65421680}, + {0.42200355, 0.43374657, 0.14424988}, + {0.30051665, 0.05120278, 0.64828057}, + {0.25088902, 0.66214977, 0.08696121}, + {0.06194140, 0.27103811, 0.66702048}, + {0.72275939, 0.20036182, 0.07687879}, + {0.64981714, 0.09746043, 0.25272243}, + {0.03670424, 0.69158962, 0.27170614}, + {0.28361681, 0.07963785, 0.63674534}, + {0.33404969, 0.61522431, 0.05072600}, + {0.40934759, 0.14334693, 0.44730548}, + {0.18115458, 0.66156829, 0.15727714}, + {0.12212382, 0.34779570, 0.53008048}, + {0.53380322, 0.27246895, 0.19372783}, + {0.64522998, 0.01711943, 0.33765059}, + {0.09252528, 0.89571913, 0.01175559}, + {0.14249520, 0.29163504, 0.56586976}, + {0.44508584, 0.36394138, 0.19097278}, + {0.36854563, 0.28925850, 0.34219587}, + {0.23714459, 0.39839435, 0.36446106}, + {0.01161403, 0.17781180, 0.81057417}, + {0.81252136, 0.16329852, 0.02418012}, + {0.78848435, 0.00015018, 0.21136546}, + {0.01254582, 0.92905218, 0.05840200}, + {0.23085637, 0.27458481, 0.49455882}, + {0.37003918, 0.42218906, 0.20777176}, + {0.45024456, 0.26024872, 0.28950672}, + {0.13461870, 0.46572118, 0.39966012}, + {0.09568060, 0.12180175, 0.78251765}, + {0.62656969, 0.30714723, 0.06628307}, + {0.54861319, 0.13276215, 0.31862466}, + {0.11887325, 0.64619030, 0.23493644}, + {0.18948259, 0.05067674, 0.75984066}, + {0.40454388, 0.52324937, 0.07220675}, + {0.33263740, 0.13121882, 0.53614379}, + {0.29037173, 0.53950416, 0.17012411}, + {0.03574861, 0.39595055, 0.56830085}, + {0.66379106, 0.20214016, 0.13406878}, + {0.74883632, 0.10187882, 0.14928486}, + {0.05873302, 0.55569368, 0.38557330}, + {0.26012063, 0.17519029, 0.56468907}, + {0.29624352, 0.56193058, 0.14182590}, + {0.42022474, 0.06972676, 0.51004850}, + {0.19743123, 0.74688264, 0.05568613}, + {0.11163233, 0.24122366, 0.64714401}, + {0.52863305, 0.33412098, 0.13724597}, + {0.56934041, 0.07519494, 0.35546465}, + {0.06966320, 0.80052515, 0.12981165}, + {0.16264765, 0.39240943, 0.44494292}, + {0.49383070, 0.26702800, 0.23914130}, + {0.35859181, 0.21465364, 0.42675454}, + {0.22007784, 0.50374629, 0.27617587}, + {0.02106424, 0.05698114, 0.92195462}, + {0.82610383, 0.17259021, 0.00130596}, + {0.88979418, 0.04060211, 0.06960370}, + {0.02783713, 0.66322450, 0.30893837}, + {0.21161353, 0.01571390, 0.77267257}, + {0.34458805, 0.62964267, 0.02576928}, + {0.48011698, 0.07982212, 0.44006090}, + {0.15255870, 0.71413914, 0.13330216}, + {0.07887744, 0.44623393, 0.47488863}, + {0.58606714, 0.22789712, 0.18603574}, + {0.51475924, 0.04019846, 0.44504230}, + {0.10146980, 0.80640087, 0.09212934}, + {0.20878503, 0.24261614, 0.54859883}, + {0.43173781, 0.42295975, 0.14530244}, + {0.30847232, 0.29219966, 0.39932802}, + {0.26973070, 0.44749009, 0.28277921}, + {0.05128846, 0.20523387, 0.74347767}, + {0.71701920, 0.22061320, 0.06236759}, + {0.69198032, 0.06937895, 0.23864074}, + {0.04265677, 0.75300864, 0.20433459}, + {0.28112109, 0.27769243, 0.44118648}, + {0.31921562, 0.38849754, 0.29228684}, + {0.39279721, 0.15350470, 0.45369810}, + {0.17902601, 0.56525353, 0.25572046}, + {0.12857214, 0.09527275, 0.77615511}, + {0.56447838, 0.39985163, 0.03566999}, + {0.60927291, 0.17521490, 0.21551219}, + {0.08573032, 0.46508628, 0.44918340}, + {0.14517521, 0.13865896, 0.71616583}, + {0.46245073, 0.45660097, 0.08094830}, + {0.38367568, 0.02845870, 0.58786562}, + {0.24005796, 0.74555625, 0.01438578}, + {0.00541770, 0.31312798, 0.68145433}, + {0.75140083, 0.15578997, 0.09280920}, + {0.76478783, 0.04584958, 0.18936259}, + {0.00636664, 0.74940485, 0.24422850}, + {0.24436353, 0.31486656, 0.44076991}, + {0.38516864, 0.37069689, 0.24413447}, + {0.46963028, 0.15117586, 0.37919386}, + {0.14551144, 0.61433574, 0.24015282}, + {0.08992623, 0.07002544, 0.84004834}, + {0.61002495, 0.34622587, 0.04374919}, + {0.56046086, 0.21034739, 0.22919175}, + {0.12586480, 0.47279218, 0.40134302}, + {0.17689407, 0.10863813, 0.71446780}, + {0.38890953, 0.49946420, 0.11162627}, + {0.31715206, 0.00964574, 0.67320221}, + {0.27727121, 0.68718320, 0.03554560}, + {0.04116187, 0.33275098, 0.62608715}, + {0.68317891, 0.20807902, 0.10874207} }; - Kokkos::View samples("sobol_bary", 100); + constexpr int npts = sizeof(data) / sizeof(data[0]); + + Kokkos::View samples("sobol_bary", npts); auto host = Kokkos::create_mirror_view(samples); - for (int i = 0; i < 100; ++i) { + for (int i = 0; i < npts; ++i) { host(i, 0) = data[i][0]; host(i, 1) = data[i][1]; host(i, 2) = data[i][2]; @@ -130,6 +430,24 @@ Kokkos::View make_sobol_barycentric_samples() return samples; } +void write_barycentric_samples_file( + const std::string& filename, + const Kokkos::View& bary) +{ + const auto bary_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), bary); + + std::ofstream out(filename); + REQUIRE(out.is_open()); + out << "l0 l1 l2\n"; + for (std::size_t i = 0; i < bary.extent(0); ++i) { + out << bary_h(i, 0) << " " + << bary_h(i, 1) << " " + << bary_h(i, 2) << "\n"; + } + out.close(); +} + double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) { const auto elem_areas = Omega_h::measure_elements_real(&mesh); @@ -155,13 +473,12 @@ double integrate_linear_field(Omega_h::Mesh& mesh, const Omega_h::Reals& u) return integral; } - Omega_h::Reals evaluate_field_at_barycentric_samples( Omega_h::Mesh& mesh, - const Kokkos::View& bary, + const Kokkos::View& bary, const std::function& f) { - const int npoints_each_tri = bary.extent(0); + const int npoints_each_tri = bary.extent(0); const auto bary_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), bary); const auto coords_h = Omega_h::HostRead(mesh.coords()); @@ -183,9 +500,9 @@ Omega_h::Reals evaluate_field_at_barycentric_samples( const int base = e * npoints_each_tri; for (int i = 0; i < npoints_each_tri; ++i) { - const double l0 = bary[i][0]; - const double l1 = bary[i][1]; - const double l2 = bary[i][2]; + const double l0 = bary_h(i, 0); + const double l1 = bary_h(i, 1); + const double l2 = bary_h(i, 2); const double x = l0 * x0 + l1 * x1 + l2 * x2; const double y = l0 * y0 + l1 * y1 + l2 * y2; @@ -197,7 +514,6 @@ Omega_h::Reals evaluate_field_at_barycentric_samples( return Omega_h::Reals(vals_h); } - TEST_CASE("monte carlo projection preserves constant and linear fields", "[transfer][mc_projection]") { @@ -217,6 +533,9 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", auto bary = make_sobol_barycentric_samples(); const int npoints_each_tri = bary.extent(0); + const std::string sobol_filename = "tmp_mc_barycentric_samples.txt"; + write_barycentric_samples_file(sobol_filename, bary); + SECTION("constant field is preserved and conserved") { const double c = 2.0; @@ -226,19 +545,19 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", auto projected = pcms::solveGalerkinProjectionMC( target_mesh, field_values_at_points, npoints_each_tri, - pcms::SamplingMethod::SOBOL,sobol_filename); + pcms::SamplingMethod::SOBOL, sobol_filename); auto projected_h = Omega_h::HostRead(projected); REQUIRE(static_cast(projected.size()) == target_mesh.nverts()); for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - REQUIRE(projected_h[i] == Catch::Approx(c).margin(1e-10)); + REQUIRE(projected_h[i] == Catch::Approx(c).margin(4e-2)); } - const double exact_integral = c * 1.0; // unit square area = 1 + const double exact_integral = c * 1.0; const double projected_integral = integrate_linear_field(target_mesh, projected); - REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(1e-10)); + REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(1e-2)); } SECTION("linear field is reproduced approximately") @@ -255,12 +574,12 @@ TEST_CASE("monte carlo projection preserves constant and linear fields", for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { const double expected = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; - REQUIRE(projected_h[i] == Catch::Approx(expected).margin(1e-6)); + REQUIRE(projected_h[i] == Catch::Approx(expected).margin(4e-2)); } - const double exact_integral = 1.0; // integral of x+y over unit square + const double exact_integral = 1.0; const double projected_integral = integrate_linear_field(target_mesh, projected); - REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(1e-6)); + REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(5e-4)); } } From a92a3d1bd07868677f5cb90a227dc40da74685d1 Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 23:34:59 -0400 Subject: [PATCH 43/44] change triangular element to just element --- src/pcms/transfer/calculate_load_vector.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp index 7e22f910..f2831831 100644 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ b/src/pcms/transfer/calculate_load_vector.hpp @@ -62,7 +62,7 @@ PetscErrorCode calculateLoadVectorMI(Omega_h::Mesh& target_mesh, * @brief Assembles the global load vector using Monte Carlo integration. * * This function computes the unassembled local load vector contributions for - * each triangular element in the target mesh using Monte Carlo integration and + * each element in the target mesh using Monte Carlo integration and * then assembles them into a global PETSc vector in COO format. * * The element-local contributions are computed from source-field values already From 228d58a49c8c3b8162f19e34269ff0f8209f1fba Mon Sep 17 00:00:00 2001 From: abhiyanpaudel Date: Thu, 19 Mar 2026 23:36:16 -0400 Subject: [PATCH 44/44] minor changes --- src/pcms/transfer/load_vector_integrator_mc.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp index c96cf115..f8435496 100644 --- a/src/pcms/transfer/load_vector_integrator_mc.cpp +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -4,7 +4,7 @@ namespace pcms { // TODO:: create a function that can sample sobol sequences instead of -// generating in python and using +// reading from file Kokkos::View read_sobol_barycentric_samples_from_file( std::string file_path) {