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) diff --git a/src/pcms/transfer/calculate_load_vector.cpp b/src/pcms/transfer/calculate_load_vector.cpp index f351a4ef..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, @@ -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); diff --git a/src/pcms/transfer/calculate_load_vector.hpp b/src/pcms/transfer/calculate_load_vector.hpp index b1968338..f2831831 100644 --- a/src/pcms/transfer/calculate_load_vector.hpp +++ b/src/pcms/transfer/calculate_load_vector.hpp @@ -13,14 +13,18 @@ #define PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP #include #include +#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 +44,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 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, Vec* loadVec_out, + const std::string sobol_filename = ""); } // namespace pcms #endif // PCMS_TRANSFER_CALCULATE_LOAD_VECTOR_HPP diff --git a/src/pcms/transfer/conservative_projection_solver.cpp b/src/pcms/transfer/conservative_projection_solver.cpp index 148454f8..1c3b53b6 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,8 +104,8 @@ Omega_h::Reals solveGalerkinProjection(Omega_h::Mesh& target_mesh, CHKERRABORT(PETSC_COMM_WORLD, ierr); Vec vec; - 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); Vec x = solveLinearSystem(mass, vec); @@ -123,15 +122,46 @@ 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 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) +{ + + 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); + 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 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 +171,26 @@ 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 d3e46a37..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 @@ -20,53 +20,224 @@ #include #include +#include +#include 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. + * 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. + * + * 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. + * 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 Solve the conservative Galerkin projection using Monte Carlo + * integration of the load vector on the target mesh. * - * @return A vector of nodal values on the target mesh after projection - * (Omega_h::Reals). + * 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. * + * @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, + const int npoints_each_tri, SamplingMethod method, + const std::string sobol_filename = ""); -Omega_h::Reals solveGalerkinProjection(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); -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 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 diff --git a/src/pcms/transfer/load_vector_integrator.hpp b/src/pcms/transfer/load_vector_integrator.hpp index 4e1d2fad..501a3580 100644 --- a/src/pcms/transfer/load_vector_integrator.hpp +++ b/src/pcms/transfer/load_vector_integrator.hpp @@ -23,10 +23,23 @@ #include #include #include +#include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include namespace pcms { + +enum class SamplingMethod { SOBOL, UNIFORM }; /** * @brief Computes the load vector for each target element in the conservative * field transfer. @@ -136,9 +149,10 @@ 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); + /// Holds projection and conservation error metrics returned by /// evaluate_pro_and_cons_errors(). struct Errors @@ -215,6 +229,130 @@ 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& 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 diff --git a/src/pcms/transfer/load_vector_integrator_mc.cpp b/src/pcms/transfer/load_vector_integrator_mc.cpp new file mode 100644 index 00000000..f8435496 --- /dev/null +++ b/src/pcms/transfer/load_vector_integrator_mc.cpp @@ -0,0 +1,243 @@ +#include "pcms/transfer/load_vector_integrator.hpp" +#include "pcms/utility/assert.h" + +namespace pcms +{ +// TODO:: create a function that can sample sobol sequences instead of +// reading from file +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; + const size_t ncols = 3; + 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); + ++nrows; + } + + Kokkos::View device_samples("device_data", nrows); + + auto host_samples = Kokkos::create_mirror_view(device_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].element_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 +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 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) { + 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; + } + return elmVolume / npoints_each_tri * Omega_h::Vector<3>{sum0, sum1, sum2}; +} + +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) +{ + + 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) { + 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]; + } + }); + + return elmLoadVector; +} +} // namespace pcms diff --git a/src/pcms/transfer/load_vector_integrator.cpp b/src/pcms/transfer/load_vector_integrator_mi.cpp similarity index 99% rename from src/pcms/transfer/load_vector_integrator.cpp rename to src/pcms/transfer/load_vector_integrator_mi.cpp index 477f744b..efb25969 100644 --- a/src/pcms/transfer/load_vector_integrator.cpp +++ b/src/pcms/transfer/load_vector_integrator_mi.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, diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c27a6de2..7063bd73 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() @@ -408,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}) diff --git a/test/test_load_vector.cpp b/test/test_load_vector.cpp index 1a38df50..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::buildLoadVector(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); @@ -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); diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index 339e8c93..d220a3ad 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()); @@ -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 = diff --git a/test/test_monte_carlo_field_transfer.cpp b/test/test_monte_carlo_field_transfer.cpp new file mode 100644 index 00000000..a7721756 --- /dev/null +++ b/test/test_monte_carlo_field_transfer.cpp @@ -0,0 +1,585 @@ +#include +#include +#include +#include +#include +#include + +#include +#include + +#include + +Kokkos::View make_sobol_barycentric_samples() +{ + static constexpr double data[][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}, + {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} + }; + + 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 < npts; ++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; +} + +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); + 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; +} + +Omega_h::Reals evaluate_field_at_barycentric_samples( + Omega_h::Mesh& mesh, + const Kokkos::View& bary, + const std::function& f) +{ + 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()); + + 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]; + 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_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; + + vals_h[base + i] = f(x, y); + } + } + + return Omega_h::Reals(vals_h); +} + +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); + + 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; + + auto field_values_at_points = evaluate_field_at_barycentric_samples( + target_mesh, bary, [c](double, double) { return c; }); + + auto projected = pcms::solveGalerkinProjectionMC( + target_mesh, field_values_at_points, npoints_each_tri, + 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(4e-2)); + } + + 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-2)); + } + + 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; }); + + auto projected = pcms::solveGalerkinProjectionMC( + target_mesh, field_values_at_points, npoints_each_tri, + pcms::SamplingMethod::SOBOL, 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(4e-2)); + } + + const double exact_integral = 1.0; + const double projected_integral = + integrate_linear_field(target_mesh, projected); + REQUIRE(projected_integral == Catch::Approx(exact_integral).margin(5e-4)); + } +} diff --git a/test/test_monte_carlo_projection_functions.cpp b/test/test_monte_carlo_projection_functions.cpp new file mode 100644 index 00000000..8b745952 --- /dev/null +++ b/test/test_monte_carlo_projection_functions.cpp @@ -0,0 +1,331 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +constexpr double tol = 1e-12; + +OMEGA_H_INLINE 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"; + out.close(); + + 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).element_id >= 0); + CHECK(host_results(1).element_id >= 0); + CHECK(host_results(2).element_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, 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, 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, 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)); + } +}