From ec0ca918267f97ab4c0ce8b5a326af4b2be5d5fd Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Thu, 6 Aug 2026 09:36:25 -0400 Subject: [PATCH] test: share common setup helpers via field_test_utils.h Factor the repeated mesh construction, function-space creation, field sampling and evaluation boilerplate out of the individual test translation units and into field_test_utils.h: - BuildUnitSquare / MakeP1Space / MakeP0Space for the unit-square simplex meshes used across the transfer and evaluation tests - IntegrateP0Field / IntegrateP1Field for conservation checks - EvaluateAndAssemble for the sample-then-assemble pattern shared by the linear form integrator tests - StandardOutsideCoords2D for out-of-bounds policy tests SetField now dispatches on the arity of the supplied callable so the same helper covers 2D and 3D layouts, and it samples the layout's DOF coordinates directly instead of copying them into a fixed-rank view. The evaluation helpers use rank-2 output views through MakeRank2View rather than hand-built Rank2Views over rank-1 storage. No behavioral change to the tests themselves. --- test/field_test_utils.h | 202 +++++++++++++---- test/test_eqdsk.cpp | 37 +--- test/test_field_evaluation.cpp | 36 +-- test/test_field_exchange_planner.cpp | 4 +- test/test_field_interpolation.cpp | 30 +-- test/test_interpolation_class.cpp | 118 ++-------- test/test_interpolation_on_ltx_mesh.cpp | 48 ++-- test/test_localization_factory.cpp | 13 +- .../test_mesh_intersection_field_transfer.cpp | 209 +++++------------- ...st_omega_h_intersection_rhs_integrator.cpp | 138 +++--------- test/test_omega_h_lagrange_field.cpp | 38 +--- test/test_omega_h_mass_integrator.cpp | 42 +--- test/test_omega_h_mc_rhs_integrator.cpp | 108 ++------- test/test_point_evaluator.cpp | 96 +++----- ...lynomial_reconstruction_mls_evaluation.cpp | 112 +++------- test/test_spr_meshfields.cpp | 18 +- test/test_svd_serial.cpp | 4 +- test/test_uniform_grid_field.cpp | 82 +++---- test/test_xgc_reverse_classification.cpp | 2 +- 19 files changed, 440 insertions(+), 897 deletions(-) diff --git a/test/field_test_utils.h b/test/field_test_utils.h index 5aaac78e..2ddd2f45 100644 --- a/test/field_test_utils.h +++ b/test/field_test_utils.h @@ -3,6 +3,7 @@ #include #include +#include "pcms/configuration.h" #include "pcms/field/field.h" #include "pcms/field/field_data.h" #include "pcms/field/field_evaluator_factory.h" @@ -11,10 +12,22 @@ #include "pcms/field/out_of_bounds_policy.h" #include "pcms/coupler/field_serializer.h" #include "pcms/field/coordinate_system.h" -#include "pcms/localization/adj_search.hpp" #include "pcms/utility/arrays.h" #include "pcms/utility/memory_spaces.h" +#ifdef PCMS_ENABLE_OMEGA_H +#include +#include +#include +#include "pcms/field/function_space/lagrange.h" +#include "pcms/localization/adj_search.hpp" +#endif +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +#include "pcms/transfer/linear_form_integrator.hpp" +#endif #include +#include +#include +#include #include #include @@ -36,6 +49,62 @@ inline std::vector StandardEvalCoords2D() return {0.1, 0.2, 0.5, 0.5, 0.7, 0.3, 0.9, 0.1, 0.2, 0.8}; } +// Points strictly outside a unit [0,1]^2 box mesh. +inline std::vector StandardOutsideCoords2D() +{ + return {-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; +} + +#ifdef PCMS_ENABLE_OMEGA_H +// Builds a unit-square 2D simplex mesh from the given element connectivity and +// adds the geometric classification tags required to build an Omega_h-backed +// Lagrange function space. +inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, + const Omega_h::LOs& ev2v) +{ + const Omega_h::Reals coords({ + 0.0, 0.0, // v0 + 1.0, 0.0, // v1 + 1.0, 1.0, // v2 + 0.0, 1.0 // v3 + }); + Omega_h::Mesh mesh(&lib); + Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); + for (Omega_h::Int dim = 0; dim <= 2; ++dim) { + mesh.add_tag( + dim, "class_dim", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::I8(dim))); + mesh.add_tag( + dim, "class_id", 1, + Omega_h::Read(mesh.nents(dim), Omega_h::ClassId(0))); + } + return mesh; +} + +// Convenience overload: diagonal=0 splits the square along vertices (1,3), +// diagonal=1 along (0,2). +inline Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) +{ + return BuildUnitSquare(lib, (diagonal == 0) + ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) + : Omega_h::LOs({0, 1, 2, 0, 2, 3})); +} + +inline std::shared_ptr MakeP1Space( + Omega_h::Mesh& mesh, const std::string& global_id_name = "global") +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 1, 1, CoordinateSystem::Cartesian, global_id_name, + LagrangeFunctionSpace::Backend::OmegaH); +} + +inline std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) +{ + return LagrangeFunctionSpace::FromMesh( + mesh, 0, 1, CoordinateSystem::Cartesian, "global", + LagrangeFunctionSpace::Backend::OmegaH); +} + inline bool AreArraysEqualUnordered( const Omega_h::HostRead& array1, const Omega_h::HostRead& array2, int start, int end) @@ -76,6 +145,40 @@ inline std::vector CopyOmegaHRealsToVector(const Omega_h::Reals& coords) coords_read.data() + coords_read.size()); } +inline double IntegrateP0Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + integral += measures_h[e] * values[e]; + } + return integral; +} + +inline double IntegrateP1Field(Omega_h::Mesh& mesh, const Field& field) +{ + const auto values = FlattenToRank1View(field.GetDOFHolderDataHost()); + const auto measures = Omega_h::measure_elements_real(&mesh); + const auto measures_h = Omega_h::HostRead(measures); + const auto elem_verts_h = + Omega_h::HostRead(mesh.ask_elem_verts()); + const int verts_per_elem = mesh.dim() + 1; + + double integral = 0.0; + for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { + double avg = 0.0; + for (int k = 0; k < verts_per_elem; ++k) { + avg += values[elem_verts_h[verts_per_elem * e + k]]; + } + integral += measures_h[e] * (avg / verts_per_elem); + } + return integral; +} +#endif // PCMS_ENABLE_OMEGA_H + // Copy coordinates from device memory to a host view. // This handles potential layout mismatches between host and device memory // spaces. @@ -86,14 +189,8 @@ inline Kokkos::View CopyCoordinatesToHost( auto coords_view = Kokkos::View("coords_view", nents, dim); auto coords_view_device = - Kokkos::create_mirror(DeviceMemorySpace(), coords_view); - Kokkos::parallel_for( - "copy_coords_to_host_view", Kokkos::RangePolicy<>(0, nents), - KOKKOS_LAMBDA(int i) { - for (int d = 0; d < dim; ++d) { - coords_view_device(i, d) = coords_device(i, d); - } - }); + Kokkos::create_mirror_view(DeviceMemorySpace(), coords_view); + ConvertMismatchLayoutView2D(coords_view_device, coords_device); Kokkos::deep_copy(coords_view, coords_view_device); return coords_view; } @@ -122,44 +219,49 @@ inline std::vector EvaluateReferenceFunction(const std::vector& pts, expected_host.data() + expected_host.extent(0)); } -template +template struct SetFieldFunctor { - Kokkos::View data; - Kokkos::View coords; + DataView data; + CoordsView coords; Func f; - SetFieldFunctor(Kokkos::View data_, - Kokkos::View coords_, Func f_) + SetFieldFunctor(DataView data_, CoordsView coords_, Func f_) : data(data_), coords(coords_), f(f_) { } KOKKOS_INLINE_FUNCTION - void operator()(int i) const { data(i) = f(coords(i, 0), coords(i, 1)); } + void operator()(int i) const + { + if constexpr (std::is_invocable_v) { + data(i) = f(coords(i, 0), coords(i, 1), coords(i, 2)); + } else { + data(i) = f(coords(i, 0), coords(i, 1)); + } + } }; -// Set scalar DOF data by sampling func at each DOF-holder coordinate. +// Set scalar DOF data by sampling func at each DOF-holder coordinate. The +// arity of func selects the spatial dimension: func(x, y) for 2D layouts, +// func(x, y, z) for 3D. template inline void SetField(const FieldLayout& layout, FieldData& field, Func func) { using MemorySpace = typename ExecutionSpace::memory_space; + static_assert(std::is_invocable_v || + std::is_invocable_v, + "SetField requires func(x, y) or func(x, y, z)"); + auto dof_coords = layout.GetDOFHolderCoordinates().GetValues(); int n = static_cast(dof_coords.extent(0)); - Kokkos::View coords_device("coords_device", n); - Kokkos::parallel_for( - "field_test_utils_copy_coords", Kokkos::RangePolicy(0, n), - KOKKOS_LAMBDA(int i) { - coords_device(i, 0) = dof_coords(i, 0); - coords_device(i, 1) = dof_coords(i, 1); - }); Kokkos::View data_device("data_device", n); Kokkos::parallel_for("field_test_utils_set_field", Kokkos::RangePolicy(0, n), - SetFieldFunctor{data_device, coords_device, func}); + SetFieldFunctor{data_device, dof_coords, func}); auto data_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), data_device); @@ -253,12 +355,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, { int n = static_cast(pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - // Rank2View out(eval.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -267,8 +365,8 @@ void CheckEvaluation(const PointEvaluator& evaluator, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i) << " expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + << " got=" << out_host(i, 0) << " expected=" << expected[i]); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } } @@ -295,16 +393,13 @@ inline void CheckFillMode(const PointEvaluator& evaluator, { int n = static_cast(outside_pts.size()) / 2; - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator.Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator.Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } @@ -339,11 +434,8 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, auto evaluator = factory->template CreatePointEvaluator( EvaluationRequest::FromCoordinates(device_coords.coordinate_view, policy)); - Kokkos::View out_device("out_device", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out(out_device.data(), n, 1); - evaluator->Evaluate(field, out); + Kokkos::View out_device("out_device", n, 1); + evaluator->Evaluate(field, MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); @@ -352,16 +444,34 @@ void CheckEvaluationWithFill(const Factory& factory, const Field& field, for (int i = 0; i < n; ++i) { INFO("Point " << i << " (" << pts[2 * static_cast(i)] << ", " << pts[2 * static_cast(i) + 1] << ")" - << " got=" << out_host(i)); + << " got=" << out_host(i, 0)); if (is_inside[i]) { INFO(" expected=" << expected[i]); - REQUIRE(out_host(i) == Catch::Approx(expected[i]).margin(abs_tol)); + REQUIRE(out_host(i, 0) == Catch::Approx(expected[i]).margin(abs_tol)); } else { - REQUIRE(out_host(i) == fill_value); + REQUIRE(out_host(i, 0) == fill_value); } } } +#if defined(PCMS_ENABLE_PETSC) && defined(PCMS_ENABLE_MESHFIELDS) +// Evaluates source_field at the integrator's sample points and assembles the +// load vector. +inline void EvaluateAndAssemble( + LinearFormIntegrator& integrator, + const std::shared_ptr& source_space, + const Field& source_field) +{ + const auto& pts = integrator.GetIntegrationPoints(); + const std::size_t npts = pts.GetValues().extent(0); + auto evaluator = source_space->CreatePointEvaluator( + EvaluationRequest::FromCoordinates(pts)); + Kokkos::View sampled("sampled", npts, 1); + evaluator->Evaluate(source_field, MakeRank2View(sampled)); + integrator.Assemble(MakeConstRank2View(sampled)); +} +#endif // PCMS_ENABLE_PETSC && PCMS_ENABLE_MESHFIELDS + } // namespace pcms::test #endif // PCMS_TEST_FIELD_TEST_UTILS_H diff --git a/test/test_eqdsk.cpp b/test/test_eqdsk.cpp index d546fbbe..c51f46b1 100644 --- a/test/test_eqdsk.cpp +++ b/test/test_eqdsk.cpp @@ -7,6 +7,7 @@ #include #include #include +#include "field_test_utils.h" #include #include @@ -212,42 +213,24 @@ TEST_CASE("EQDSKData with SplineFunctionSpace") const int num_eval_points = 3; - auto eval_coords_host = Kokkos::View( - "eval_coords_host", num_eval_points, 2); - for (size_t i = 0; i < num_eval_points; ++i) { - eval_coords_host(i, 0) = eval_coords[2 * i]; // R - eval_coords_host(i, 1) = eval_coords[2 * i + 1]; // Z - } - - auto eval_coords_device = - Kokkos::View("eval_coords_device", - num_eval_points, 2); - pcms::DeepCopyMismatchLayouts(eval_coords_device, eval_coords_host); - - auto coords_view = pcms::MakeRank2View(eval_coords_device); - auto coord_view = pcms::CoordinateView{ - CoordinateSystem::Cartesian, coords_view}; - auto eval_request = pcms::EvaluationRequest::FromCoordinates(coord_view); + auto device_coords = pcms::test::CreateDeviceCoordinateView( + eval_coords, CoordinateSystem::Cartesian); + auto eval_request = + pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view); auto evaluator = spline_space->CreatePointEvaluator(eval_request); - auto eval_results_1d = Kokkos::View( - "eval_results", num_eval_points); - - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View - eval_results(eval_results_1d.data(), num_eval_points, 1); - - evaluator->Evaluate(psi_field, eval_results); + Kokkos::View eval_results( + "eval_results", num_eval_points, 1); + evaluator->Evaluate(psi_field, pcms::MakeRank2View(eval_results)); auto results_host = Kokkos::create_mirror_view_and_copy( - pcms::HostMemorySpace(), eval_results_1d); + pcms::HostMemorySpace(), eval_results); // Verify that all results are finite for (int i = 0; i < num_eval_points; ++i) { - REQUIRE(std::isfinite(results_host(i))); + REQUIRE(std::isfinite(results_host(i, 0))); } } } diff --git a/test/test_field_evaluation.cpp b/test/test_field_evaluation.cpp index ff238e11..85af38b4 100644 --- a/test/test_field_evaluation.cpp +++ b/test/test_field_evaluation.cpp @@ -34,10 +34,10 @@ TEST_CASE("evaluate linear 2d omega_h_field") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } #ifdef PCMS_ENABLE_MESHFIELDS @@ -50,35 +50,11 @@ TEST_CASE("evaluate quadratic 2d meshfields_field") mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - // Quadratic DOF holders span vertices and edge midpoints; set them inline. - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - test_f[i] = sin_f(mesh_coords[2 * static_cast(i)], - mesh_coords[2 * static_cast(i) + 1]); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto ep = Omega_h::gather_verts<2>(edge_verts, i); - Real cx = (mesh_coords[2 * static_cast(ep[0])] + - mesh_coords[2 * static_cast(ep[1])]) / - 2; - Real cy = (mesh_coords[2 * static_cast(ep[0]) + 1] + - mesh_coords[2 * static_cast(ep[1]) + 1]) / - 2; - test_f[nverts + i] = sin_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); + // Quadratic DOF holders span vertices and edge midpoints; the layout's DOF + // coordinates cover both, so SetField samples sin_f at every holder. auto field = factory->CreateFunction(); - field.GetData().SetDOFHolderDataHost( - pcms::Rank2View(test_f_host.data(), - test_f_host.size(), 1)); + pcms::test::SetField( + field, OMEGA_H_LAMBDA(Real x, Real y) { return sin_f(x, y); }); pcms::test::CheckEvaluation( factory, field, kEvalCoords, diff --git a/test/test_field_exchange_planner.cpp b/test/test_field_exchange_planner.cpp index 61b5bc54..99e6181b 100644 --- a/test/test_field_exchange_planner.cpp +++ b/test/test_field_exchange_planner.cpp @@ -51,9 +51,7 @@ TEST_CASE("GID messages insert headers around compact field payloads", "gid_message", plan.msg_size + 2 * static_cast(pcms::ent_offsets_len)); pcms::GenericFieldExchangePlanner planner; - planner.FillGidMessage(layout, plan, - pcms::Rank1View( - message.data(), message.size())); + planner.FillGidMessage(layout, plan, pcms::MakeRank1View(message)); const pcms::GO expected[] = { 0, 2, 2, 2, 2, 0, 2, 0, 2, 2, 2, 2, 1, 3, diff --git a/test/test_field_interpolation.cpp b/test/test_field_interpolation.cpp index bc88ab97..f8e2432d 100644 --- a/test/test_field_interpolation.cpp +++ b/test/test_field_interpolation.cpp @@ -55,35 +55,11 @@ TEST_CASE("interpolate quadratic 2d meshfields_field") auto factory2 = pcms::LagrangeFunctionSpace::FromMesh( mesh, 2, 1, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::MeshFields); - auto layout = factory2->GetLayout(); - const auto nverts = mesh.nents(0); - const auto nedges = mesh.nents(1); - auto mesh_coords = mesh.coords(); - auto edge_verts = mesh.ask_verts_of(1); - Omega_h::Write test_f(nverts + nedges); - Omega_h::parallel_for( - nverts, OMEGA_H_LAMBDA(int i) { - Real x = mesh_coords[2 * i + 0]; - Real y = mesh_coords[2 * i + 1]; - test_f[i] = interpolation_linear_f(x, y); - }); - Omega_h::parallel_for( - nedges, OMEGA_H_LAMBDA(int i) { - auto endpoints = Omega_h::gather_verts<2>(edge_verts, i); - Real x0 = mesh_coords[2 * endpoints[0] + 0]; - Real y0 = mesh_coords[2 * endpoints[0] + 1]; - Real x1 = mesh_coords[2 * endpoints[1] + 0]; - Real y1 = mesh_coords[2 * endpoints[1] + 1]; - Real cx = (x0 + x1) / 2; - Real cy = (y0 + y1) / 2; - test_f[nverts + i] = interpolation_linear_f(cx, cy); - }); - - Omega_h::HostWrite test_f_host(test_f); auto field = factory2->CreateFunction(); auto interpolated = factory2->CreateFunction(); - field.SetDOFHolderDataHost(pcms::Rank2View( - test_f_host.data(), test_f_host.size(), 1)); + pcms::test::SetField( + field, + OMEGA_H_LAMBDA(Real x, Real y) { return interpolation_linear_f(x, y); }); pcms::Interpolator interp(*factory2, *factory2); interp.Apply(field, interpolated); diff --git a/test/test_interpolation_class.cpp b/test/test_interpolation_class.cpp index 7a34af4d..878f6288 100644 --- a/test/test_interpolation_class.cpp +++ b/test/test_interpolation_class.cpp @@ -10,36 +10,10 @@ #include #include #include +#include "field_test_utils.h" #include #include -#include - -bool areArraysEqualUnordered(const Omega_h::HostRead& array1, - const Omega_h::HostRead& array2, - int start, int end) -{ - // Ensure the indices are valid - assert(start >= 0 && end <= array1.size() && start <= end); - assert(start >= 0 && end <= array2.size() && start <= end); - - // Use frequency maps to count occurrences of each value - std::unordered_map freq1, freq2; - - for (int i = start; i < end; ++i) { - freq1[array1[i]]++; - freq2[array2[i]]++; - } - - // Compare the frequency maps - if (freq1 != freq2) { - pcms::printError("[ERROR] Arrays differ in the range [%d, %d)\n", start, - end); - return false; - } - - return true; -} void translate_mesh(Omega_h::Mesh* mesh, Omega_h::Vector<2> translation_vector) { @@ -135,28 +109,13 @@ TEST_CASE("Test MLSMeshInterpolation") auto mls_single = pcms::MLSMeshInterpolation(source_mesh, 0.12, 15, 3, true, 0.0, 5.0); - auto source_points_reals = - pcms::get_entity_centroids(source_mesh, Omega_h::FACE); - auto source_points_host = - Omega_h::HostRead(source_points_reals); - auto source_points_host_write = - Omega_h::HostWrite(source_points_host.size()); - for (int i = 0; i < source_points_host.size(); i++) { - source_points_host_write[i] = source_points_host[i]; - } - auto source_points_view = pcms::Rank1View( - source_points_host_write.data(), source_points_host_write.size()); - - auto target_points_reals = source_mesh.coords(); - auto target_points_host = - Omega_h::HostRead(target_points_reals); - auto target_points_host_write = - Omega_h::HostWrite(target_points_host.size()); - for (int i = 0; i < target_points_host.size(); i++) { - target_points_host_write[i] = target_points_host[i]; - } - auto target_points_view = pcms::Rank1View( - target_points_host_write.data(), target_points_host_write.size()); + auto source_points_vec = pcms::test::CopyOmegaHRealsToVector( + pcms::get_entity_centroids(source_mesh, Omega_h::FACE)); + auto source_points_view = pcms::make_array_view(source_points_vec); + + auto target_points_vec = + pcms::test::CopyOmegaHRealsToVector(source_mesh.coords()); + auto target_points_view = pcms::make_array_view(target_points_vec); REQUIRE(source_mesh.dim() == 2); pcms::printInfo("Point cloud based search...\n"); auto point_mls = pcms::MLSPointCloudInterpolation( @@ -173,14 +132,11 @@ TEST_CASE("Test MLSMeshInterpolation") source_mesh.nverts()); Omega_h::HostWrite exact_values_at_nodes(source_sinxcosy_node); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); - pcms::Rank1View - point_cloud_interpolatedArrayView( - point_cloud_interpolated_data_hwrite.data(), - point_cloud_interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); + auto point_cloud_interpolatedArrayView = + pcms::make_array_view(point_cloud_interpolated_data_hwrite); OMEGA_H_CHECK_PRINTF(sourceArrayView.size() == mls_single.getSourceSize(), "Source size mismatch: %zu vs %zu\n", @@ -222,33 +178,8 @@ TEST_CASE("Test MLSMeshInterpolation") ///*****************************// auto mesh_based_supports = mls_single.getSupports(); auto point_cloud_based_supports = point_mls.getSupports(); - auto mesh_based_support_ptr_host = - Omega_h::HostRead(mesh_based_supports.supports_ptr); - auto mesh_based_support_idx_host = - Omega_h::HostRead(mesh_based_supports.supports_idx); - auto point_cloud_based_support_ptr_host = - Omega_h::HostRead(point_cloud_based_supports.supports_ptr); - auto point_cloud_based_support_idx_host = - Omega_h::HostRead(point_cloud_based_supports.supports_idx); - - REQUIRE(point_cloud_based_support_idx_host.size() == - mesh_based_support_idx_host.size()); - REQUIRE(point_cloud_based_support_ptr_host.size() == - mesh_based_support_ptr_host.size()); - for (int i = 0; i < mesh_based_support_ptr_host.size(); i++) { - REQUIRE(point_cloud_based_support_ptr_host[i] == - mesh_based_support_ptr_host[i]); - } - - for (int i = 0; i < mesh_based_support_ptr_host.size() - 1; i++) { - auto start = mesh_based_support_ptr_host[i]; - auto end = mesh_based_support_ptr_host[i + 1]; - - bool isEqual = - areArraysEqualUnordered(mesh_based_support_idx_host, - point_cloud_based_support_idx_host, start, end); - REQUIRE(isEqual); - } + pcms::test::CheckSupportResultsEquivalent(point_cloud_based_supports, + mesh_based_supports); // Check if the point cloud interpolation is same as the MLS interpolation pcms::printDebugInfo("Interpolated data size: %d\n", @@ -289,10 +220,9 @@ TEST_CASE("Test MLSMeshInterpolation") Omega_h::HostWrite interpolated_data_hwrite( mls_double.getTargetSize()); - pcms::Rank1View sourceArrayView( - source_data_host_write.data(), source_data_host_write.size()); - pcms::Rank1View interpolatedArrayView( - interpolated_data_hwrite.data(), interpolated_data_hwrite.size()); + auto sourceArrayView = pcms::make_array_view(source_data_host_write); + auto interpolatedArrayView = + pcms::make_array_view(interpolated_data_hwrite); mls_double.eval(sourceArrayView, interpolatedArrayView); @@ -326,10 +256,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto target_points = source_points; - auto source_points_view = pcms::Rank1View( - source_points.data(), source_points.size()); - auto target_points_view = pcms::Rank1View( - target_points.data(), target_points.size()); + auto source_points_view = pcms::make_array_view(source_points); + auto target_points_view = pcms::make_array_view(target_points); // Degree-1 polynomial that depends on z to catch accidental 2D behavior. auto source_values = Omega_h::HostWrite(27); @@ -341,10 +269,8 @@ TEST_CASE("MLSPointCloudInterpolation honors provided dimension in eval") } auto output_values = Omega_h::HostWrite(27, "output_values"); - auto source_values_view = pcms::Rank1View( - source_values.data(), source_values.size()); - auto output_values_view = pcms::Rank1View( - output_values.data(), output_values.size()); + auto source_values_view = pcms::make_array_view(source_values); + auto output_values_view = pcms::make_array_view(output_values); auto mls = pcms::MLSPointCloudInterpolation( source_points_view, target_points_view, 3, 2.5, 10, 1, true, 0.0, 5.0); diff --git a/test/test_interpolation_on_ltx_mesh.cpp b/test/test_interpolation_on_ltx_mesh.cpp index a348acae..141ebfa8 100644 --- a/test/test_interpolation_on_ltx_mesh.cpp +++ b/test/test_interpolation_on_ltx_mesh.cpp @@ -83,16 +83,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") printf("[INFO] Degas2 Mesh loaded from %s with %d elements\n", degas2_mesh_filename.c_str(), degas2_num_elems); const auto degas2_mesh_centroids_view = - pcms::Rank1View( - degas2_mesh_centroids_host.data(), degas2_mesh_centroids_host.size()); + pcms::make_array_view(degas2_mesh_centroids_host); auto xgc_mesh_points = read_xgc_mesh_nodes(ltx_mesh_base_filename + ".node"); const int xgc_num_nodes = xgc_mesh_points.size() / 2; printf("[INFO] XGC Mesh loaded from %s with %d points\n", ltx_mesh_base_filename.c_str(), xgc_num_nodes); - const auto xgc_mesh_points_view = - pcms::Rank1View( - xgc_mesh_points.data(), xgc_mesh_points.size()); + const auto xgc_mesh_points_view = pcms::make_array_view(xgc_mesh_points); auto xgc_to_degas2_interpolator = pcms::MLSPointCloudInterpolation( xgc_mesh_points_view, degas2_mesh_centroids_view, 2, 0.000001, 10, 1, true, @@ -120,18 +117,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") // ------------------ First Interpolation ------------------ // const auto density_at_xgc_nodes_view = - pcms::Rank1View(density_at_xgc_nodes.data(), - density_at_xgc_nodes.size()); - const auto temp_at_xgc_nodes_view = - pcms::Rank1View(temp_at_xgc_nodes.data(), - temp_at_xgc_nodes.size()); + pcms::make_array_view(density_at_xgc_nodes); + const auto temp_at_xgc_nodes_view = pcms::make_array_view(temp_at_xgc_nodes); const auto density_at_degas2_centroids_view = - pcms::Rank1View( - density_at_degas2_centroids.data(), density_at_degas2_centroids.size()); + pcms::make_array_view(density_at_degas2_centroids); const auto temp_at_degas2_centroids_view = - pcms::Rank1View( - temp_at_degas2_centroids.data(), temp_at_degas2_centroids.size()); + pcms::make_array_view(temp_at_degas2_centroids); Omega_h::HostWrite interpolated_xgc_density(degas2_num_elems); Omega_h::HostWrite interpolated_xgc_temp(degas2_num_elems); @@ -139,17 +131,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") Omega_h::HostWrite interpolated_degas2_temp(xgc_num_nodes); const auto interpolated_xgc_density_view = - pcms::Rank1View( - interpolated_xgc_density.data(), interpolated_xgc_density.size()); + pcms::make_array_view(interpolated_xgc_density); const auto interpolated_xgc_temp_view = - pcms::Rank1View( - interpolated_xgc_temp.data(), interpolated_xgc_temp.size()); + pcms::make_array_view(interpolated_xgc_temp); const auto interpolated_degas2_density_view = - pcms::Rank1View( - interpolated_degas2_density.data(), interpolated_degas2_density.size()); + pcms::make_array_view(interpolated_degas2_density); const auto interpolated_degas2_temp_view = - pcms::Rank1View( - interpolated_degas2_temp.data(), interpolated_degas2_temp.size()); + pcms::make_array_view(interpolated_degas2_temp); xgc_to_degas2_interpolator.eval(density_at_xgc_nodes_view, interpolated_xgc_density_view); @@ -174,21 +162,13 @@ TEST_CASE("Test Interpolation on LTX Mesh", "[interpolation]") degas2_num_elems); const auto interpolated_back_density_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_density_at_xgc_nodes.data(), - interpolated_back_density_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_density_at_xgc_nodes); const auto interpolated_back_temp_at_xgc_nodes_view = - pcms::Rank1View( - interpolated_back_temp_at_xgc_nodes.data(), - interpolated_back_temp_at_xgc_nodes.size()); + pcms::make_array_view(interpolated_back_temp_at_xgc_nodes); const auto interpolated_back_density_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_density_at_degas2_centroids.data(), - interpolated_back_density_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_density_at_degas2_centroids); const auto interpolated_back_temp_at_degas2_centroids_view = - pcms::Rank1View( - interpolated_back_temp_at_degas2_centroids.data(), - interpolated_back_temp_at_degas2_centroids.size()); + pcms::make_array_view(interpolated_back_temp_at_degas2_centroids); degas2_to_xgc_interpolator.eval(interpolated_xgc_density_view, interpolated_back_density_at_xgc_nodes_view); diff --git a/test/test_localization_factory.cpp b/test/test_localization_factory.cpp index 571ac175..d65e8fc0 100644 --- a/test/test_localization_factory.cpp +++ b/test/test_localization_factory.cpp @@ -14,6 +14,7 @@ #include "pcms/localization/point_cloud_localization.h" #include "pcms/discretization/discretization/omega_h.hpp" #include "pcms/utility/arrays.h" +#include "pcms/utility/omega_h_array_utils.h" #include "field_test_utils.h" namespace @@ -46,17 +47,7 @@ TEST_CASE( auto target_device = pcms::test::CreateDeviceCoordinateView( target_coords, pcms::CoordinateSystem::Cartesian, dim); - auto coords_read = Omega_h::HostRead(source_coords); - - // auto coords_dev = Kokkos::create_mirror_view_and_copy( - // Kokkos::DefaultExecutionSpace{}, coords_host); - auto coords_dev = Kokkos::View( - "point_cloud_coords", mesh.nverts(), dim); - auto coords_host = Kokkos::create_mirror(pcms::HostMemorySpace(), coords_dev); - for (int i = 0; i < mesh.nverts(); ++i) - for (int d = 0; d < dim; ++d) - coords_host(i, d) = coords_read[i * dim + d]; - Kokkos::deep_copy(coords_dev, coords_host); + auto coords_dev = pcms::ConvertCoordsTo2D(source_coords, mesh.nverts(), dim); auto layout = std::make_shared( dim, coords_dev, pcms::CoordinateSystem::Cartesian); diff --git a/test/test_mesh_intersection_field_transfer.cpp b/test/test_mesh_intersection_field_transfer.cpp index fbb34ddb..a627ab81 100644 --- a/test/test_mesh_intersection_field_transfer.cpp +++ b/test/test_mesh_intersection_field_transfer.cpp @@ -7,53 +7,14 @@ #include #include +#include #include "field_test_utils.h" -#include #include namespace { -// Builds a unit-square 2D simplex mesh from the given element connectivity and -// adds the geometric classification tags required to build an Omega_h-backed -// Lagrange function space. -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, const Omega_h::LOs& ev2v) -{ - const Omega_h::Reals coords({ - 0.0, 0.0, // v0 - 1.0, 0.0, // v1 - 1.0, 1.0, // v2 - 0.0, 1.0 // v3 - }); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - for (Omega_h::Int dim = 0; dim <= 2; ++dim) { - mesh.add_tag( - dim, "class_dim", 1, - Omega_h::Read(mesh.nents(dim), Omega_h::I8(dim))); - mesh.add_tag( - dim, "class_id", 1, - Omega_h::Read(mesh.nents(dim), Omega_h::ClassId(0))); - } - return mesh; -} - -std::shared_ptr MakeP1Space( - Omega_h::Mesh& mesh, const std::string& global_id_name = "global") -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, global_id_name, - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - -std::shared_ptr MakeP0Space(Omega_h::Mesh& mesh) -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 0, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - void AddReorderedVertexGlobalIds(Omega_h::Mesh& mesh) { mesh.add_tag( @@ -68,64 +29,6 @@ void AddSparseVertexGlobalIds(Omega_h::Mesh& mesh) Omega_h::Read({102, 7, 41, 19}, "sparse_global")); } -// Integral over the mesh of an order-0 (piecewise-constant) field whose values -// are stored in element order. Exact by construction. -double IntegrateP0Field(Omega_h::Mesh& mesh, - const pcms::Field& field) -{ - const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); - const auto elem_areas = Omega_h::measure_elements_real(&mesh); - const auto elem_areas_h = Omega_h::HostRead(elem_areas); - - double integral = 0.0; - for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { - integral += elem_areas_h[e] * values[e]; - } - return integral; -} - -// Per-element centroid coordinates in element order. -std::vector> ElementCentroids(Omega_h::Mesh& mesh) -{ - const auto coords_h = Omega_h::HostRead(mesh.coords()); - const auto elem_verts_h = - Omega_h::HostRead(mesh.ask_elem_verts()); - std::vector> centroids(mesh.nelems()); - for (Omega_h::LO e = 0; e < mesh.nelems(); ++e) { - double cx = 0.0, cy = 0.0; - for (int k = 0; k < 3; ++k) { - const Omega_h::LO v = elem_verts_h[3 * e + k]; - cx += coords_h[2 * v + 0]; - cy += coords_h[2 * v + 1]; - } - centroids[e] = {cx / 3.0, cy / 3.0}; - } - return centroids; -} - -// Integral over the mesh of an order-1 nodal field whose values are stored in -// vertex order (as produced by Field::GetDOFHolderDataHost for the -// Omega_h backend), computed exactly via per-element vertex averaging. -double IntegrateP1Field(Omega_h::Mesh& mesh, - const pcms::Field& field) -{ - const auto values = pcms::FlattenToRank1View(field.GetDOFHolderDataHost()); - const auto elem_areas = Omega_h::measure_elements_real(&mesh); - const auto elem_verts = mesh.ask_elem_verts(); - const auto elem_areas_h = Omega_h::HostRead(elem_areas); - const auto elem_verts_h = Omega_h::HostRead(elem_verts); - - 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 avg = (values[v0] + values[v1] + values[v2]) / 3.0; - integral += elem_areas_h[e] * avg; - } - return integral; -} - } // namespace // Fields that already live in the target order-1 space (constants and affine @@ -139,12 +42,12 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", // Source and target triangulate the same square along opposite diagonals. Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -166,8 +69,9 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field is reproduced on target vertices and conserved") @@ -179,34 +83,31 @@ TEST_CASE("OmegaHConservativeProjection reproduces constant and linear fields", const auto target_values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), + 2); 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]; + const double expected = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); REQUIRE(target_values[i] == Catch::Approx(expected).margin(1e-9)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } -// For a field that does not live in the target space the projection is not an -// interpolation, but conservation of the integral is the defining property of -// the conservative (Galerkin) projection and must still hold exactly. This -// also exercises a second Apply with a different source field to confirm the -// cached integrators/evaluator/factorization are reused without stale state. TEST_CASE("OmegaHConservativeProjection conserves the integral", "[transfer][mesh_intersection]") { Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -218,8 +119,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", return x * x + x * y + 0.5 * y * y; }); projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); // Second Apply with a different source field — verifies cached state is // reused correctly and is not stale. @@ -227,8 +129,9 @@ TEST_CASE("OmegaHConservativeProjection conserves the integral", source, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return 2.0 * x - y + 0.5; }); projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-12)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-12)); } TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " @@ -238,13 +141,13 @@ TEST_CASE("OmegaHConservativeProjection writes reordered target GIDs in local " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); AddReorderedVertexGlobalIds(target_mesh); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh, "reordered_global"); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "reordered_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -277,13 +180,13 @@ TEST_CASE("OmegaHConservativeProjection maps sparse target GIDs to active " Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); AddSparseVertexGlobalIds(target_mesh); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh, "sparse_global"); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh, "sparse_global"); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -318,12 +221,12 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP0Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_space = pcms::test::MakeP0Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -345,8 +248,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { REQUIRE(target_values[i] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } SECTION("non-constant P0 source conserves the integral") @@ -357,8 +261,9 @@ TEST_CASE("OmegaHConservativeProjection P0 source to P1 target", projection.Apply(source, target); - REQUIRE(IntegrateP1Field(target_mesh, target) == - Catch::Approx(IntegrateP0Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP1Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP0Field(source_mesh, source)) + .margin(1e-10)); } } @@ -372,12 +277,12 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", Omega_h::Library lib; Omega_h::Mesh source_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 2, 0, 2, 3})); Omega_h::Mesh target_mesh = - BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); + pcms::test::BuildUnitSquare(lib, Omega_h::LOs({0, 1, 3, 1, 2, 3})); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP0Space(target_mesh); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP0Space(target_mesh); auto source = source_space->CreateFunction(); auto target = target_space->CreateFunction(); @@ -399,8 +304,9 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { REQUIRE(target_values[e] == Catch::Approx(c).margin(1e-10)); } - REQUIRE(IntegrateP0Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-10)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-10)); } SECTION("linear field projects to cell average and conserves the integral") @@ -414,13 +320,16 @@ TEST_CASE("OmegaHConservativeProjection P1 source to P0 target", const auto target_values = pcms::FlattenToRank1View(target.GetDOFHolderDataHost()); - const auto centroids = ElementCentroids(target_mesh); - REQUIRE(static_cast(target_values.size()) == centroids.size()); - for (std::size_t e = 0; e < centroids.size(); ++e) { - const double expected = f(centroids[e][0], centroids[e][1]); + const auto centroids_h = Omega_h::HostRead( + pcms::get_entity_centroids(target_mesh, target_mesh.dim())); + REQUIRE(static_cast(target_values.size()) == + target_mesh.nelems()); + for (Omega_h::LO e = 0; e < target_mesh.nelems(); ++e) { + const double expected = f(centroids_h[2 * e + 0], centroids_h[2 * e + 1]); REQUIRE(target_values[e] == Catch::Approx(expected).margin(1e-9)); } - REQUIRE(IntegrateP0Field(target_mesh, target) == - Catch::Approx(IntegrateP1Field(source_mesh, source)).margin(1e-9)); + REQUIRE(pcms::test::IntegrateP0Field(target_mesh, target) == + Catch::Approx(pcms::test::IntegrateP1Field(source_mesh, source)) + .margin(1e-9)); } } diff --git a/test/test_omega_h_intersection_rhs_integrator.cpp b/test/test_omega_h_intersection_rhs_integrator.cpp index 18ea11bb..f5d667d9 100644 --- a/test/test_omega_h_intersection_rhs_integrator.cpp +++ b/test/test_omega_h_intersection_rhs_integrator.cpp @@ -21,37 +21,6 @@ namespace { -// Build a unit-square 2D simplex mesh. -// diagonal=0: T0=(0,1,3), T1=(1,2,3) -// diagonal=1: T0=(0,1,2), T1=(0,2,3) -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) -{ - const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) - : Omega_h::LOs({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - // Independent reference for the assembled conservative load vector // b_j = \int phi_j^target f dx // when f is exactly representable in the target order-1 space (constant or @@ -100,17 +69,7 @@ void CheckAssembledLoadMatches( const pcms::Field& source_field, const std::unordered_map& expected) { - const auto pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator.Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); Vec vec = integrator.GetVector(); const PetscScalar* vec_array = nullptr; @@ -135,26 +94,19 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: integration points lie inside " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto integrator = pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); const auto raw_coords = integrator->GetIntegrationPoints().GetValues(); - auto raw_coords_view = Kokkos::View>( - raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); - auto raw_coords_host = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, raw_coords_view); + auto raw_coords_host = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); const std::size_t n = raw_coords_host.extent(0); REQUIRE(n > 0); @@ -172,32 +124,18 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: zero source field gives zero " "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); // Default-constructed field has zero data. auto integrator = pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space); - const auto& pts = integrator->GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - - integrator->Assemble(pcms::MakeConstRank2View(sampled)); + pcms::test::EvaluateAndAssemble(*integrator, source_space, source_field); Vec vec = integrator->GetVector(); PetscReal norm = 0.0; @@ -214,15 +152,11 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: constant field matches " // For a constant source field the assembled load vector must equal the // target consistent mass matrix applied to the constant nodal vector. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); const double c = 2.0; auto source_field = source_space->CreateFunction(); @@ -251,15 +185,11 @@ TEST_CASE( // vector must equal the target consistent mass matrix applied to the nodal // values of x + y. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -268,11 +198,11 @@ TEST_CASE( const auto target_layout = std::dynamic_pointer_cast( target_space->GetLayout()); - const auto tgt_coords_h = - Omega_h::HostRead(target_mesh.coords()); + const auto tgt_coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); std::vector g(target_mesh.nverts()); for (int i = 0; i < target_mesh.nverts(); ++i) { - g[i] = tgt_coords_h[2 * i + 0] + tgt_coords_h[2 * i + 1]; + g[i] = tgt_coords_h(i, 0) + tgt_coords_h(i, 1); } const auto expected = ExpectedLoadByGid(target_mesh, *target_layout, g); @@ -285,26 +215,22 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", "[rhs_integrator]") { Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); SECTION("multi-component source space throws") { auto source_space = pcms::LagrangeFunctionSpace::FromMesh( source_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::test::MakeP1Space(target_mesh); REQUIRE_THROWS( pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space)); } SECTION("multi-component target space throws") { - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 2, pcms::CoordinateSystem::Cartesian, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); @@ -317,18 +243,14 @@ TEST_CASE("OmegaHIntersectionRHSIntegrator: rejects invalid layouts", auto source_space = pcms::LagrangeFunctionSpace::FromMesh( source_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); - auto target_space = pcms::LagrangeFunctionSpace::FromMesh( - target_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto target_space = pcms::test::MakeP1Space(target_mesh); REQUIRE_THROWS( pcms::BuildOmegaHConservativeRHSIntegrator(*source_space, *target_space)); } SECTION("non-Cartesian target coordinate system throws") { - auto source_space = pcms::LagrangeFunctionSpace::FromMesh( - source_mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto source_space = pcms::test::MakeP1Space(source_mesh); auto target_space = pcms::LagrangeFunctionSpace::FromMesh( target_mesh, 1, 1, pcms::CoordinateSystem::Cylindrical, "global", pcms::LagrangeFunctionSpace::Backend::OmegaH); diff --git a/test/test_omega_h_lagrange_field.cpp b/test/test_omega_h_lagrange_field.cpp index c17d45b1..0dd5ffc7 100644 --- a/test/test_omega_h_lagrange_field.cpp +++ b/test/test_omega_h_lagrange_field.cpp @@ -155,10 +155,10 @@ TEST_CASE("OmegaHLagrangeField order-1: linear function evaluation") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // The same linear evaluation test run on the MeshFields-backed order-1 field @@ -173,10 +173,10 @@ TEST_CASE("MeshFieldsAdapter order-1: linear function evaluation (shared util)") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckEvaluation( factory, field, pcms::test::StandardEvalCoords2D(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") @@ -189,10 +189,10 @@ TEST_CASE("OmegaHLagrangeField order-1: out-of-bounds FILL mode") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); Real fill_value = -999.0; - std::vector outside{-0.5, 0.5, 1.5, 0.5, 0.5, -0.5, 0.5, 1.5}; + auto outside = pcms::test::StandardOutsideCoords2D(); pcms::test::CheckFillMode(factory, field, fill_value, outside); } @@ -206,7 +206,7 @@ TEST_CASE("OmegaHLagrangeField order-1: serialize / deserialize round-trip") pcms::test::SetField( field.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::CheckSerializeDeserialize(*factory->GetLayout(), field.GetData()); } @@ -270,24 +270,9 @@ TEST_CASE("OmegaHLagrangeField order-0: constant field evaluation") 1); field.GetData().SetDOFHolderDataHost(view); - auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; - auto device_coords = - pcms::test::CreateDeviceCoordinateView(pts, factory->GetCoordinateSystem()); - auto evaluator = factory->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - - Kokkos::View eval_device("eval", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out = pcms::Rank2View( - eval_device.data(), n, 1); - evaluator->Evaluate(field, out); - auto eval_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), eval_device); - - for (int i = 0; i < n; ++i) - REQUIRE(eval_host(i) == Catch::Approx(kValue)); + pcms::test::CheckEvaluation( + factory, field, pcms::test::StandardEvalCoords2D(), + OMEGA_H_LAMBDA(Real, Real) { return kValue; }); } TEST_CASE("OmegaHLagrangeField order-0: out-of-bounds FILL mode") @@ -357,7 +342,8 @@ TEST_CASE("OmegaHLagrangeField: field valid after layout destruction") } // factory goes out of scope; field keeps layout alive pcms::test::SetField( - *field, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + *field, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Just verify data was set correctly (no evaluator needed for this lifetime // test) auto data = field->GetDOFHolderDataHost(); diff --git a/test/test_omega_h_mass_integrator.cpp b/test/test_omega_h_mass_integrator.cpp index c9013f9d..c49a2e7d 100644 --- a/test/test_omega_h_mass_integrator.cpp +++ b/test/test_omega_h_mass_integrator.cpp @@ -8,6 +8,7 @@ #include #include #include +#include "field_test_utils.h" #include #include #include @@ -23,33 +24,6 @@ namespace { -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib) -{ - const 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); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - std::map, pcms::Real> BuildReferenceMassMap( Omega_h::Mesh& mesh, const pcms::FunctionSpace& space) { @@ -98,11 +72,9 @@ TEST_CASE( "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto ref = BuildReferenceMassMap(mesh, *space); @@ -128,11 +100,9 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", // each node: row_sum(i) = integral(N_i dx). For a regular mesh of unit // area with uniform nodal distribution the sum over all nodes equals 1. Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); - auto space = pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); + auto space = pcms::test::MakeP1Space(mesh); const auto layout = std::dynamic_pointer_cast( @@ -169,7 +139,7 @@ TEST_CASE("OmegaHMassIntegrator: row sums match lumped mass", TEST_CASE("OmegaHMassIntegrator: rejects invalid layouts", "[mass_integrator]") { Omega_h::Library lib; - auto mesh = BuildUnitSquare(lib); + auto mesh = pcms::test::BuildUnitSquare(lib, 0); SECTION("multi-component space throws") { diff --git a/test/test_omega_h_mc_rhs_integrator.cpp b/test/test_omega_h_mc_rhs_integrator.cpp index 43172fce..4cf8a651 100644 --- a/test/test_omega_h_mc_rhs_integrator.cpp +++ b/test/test_omega_h_mc_rhs_integrator.cpp @@ -15,69 +15,6 @@ #include #include -// --------------------------------------------------------------------------- -// Helpers -// --------------------------------------------------------------------------- - -namespace -{ - -// Build a unit-square 2D simplex mesh. -// diagonal=0: T0=(0,1,3), T1=(1,2,3) -// diagonal=1: T0=(0,1,2), T1=(0,2,3) -Omega_h::Mesh BuildUnitSquare(Omega_h::Library& lib, int diagonal) -{ - const Omega_h::Reals coords({0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0}); - Omega_h::LOs ev2v = (diagonal == 0) ? Omega_h::LOs({0, 1, 3, 1, 2, 3}) - : Omega_h::LOs({0, 1, 2, 0, 2, 3}); - Omega_h::Mesh mesh(&lib); - Omega_h::build_from_elems_and_coords(&mesh, OMEGA_H_SIMPLEX, 2, ev2v, coords); - mesh.add_tag( - 0, "class_dim", 1, - Omega_h::Read(mesh.nverts(), Omega_h::I8(0))); - mesh.add_tag( - 0, "class_id", 1, - Omega_h::Read(mesh.nverts(), Omega_h::ClassId(0))); - mesh.add_tag( - 1, "class_dim", 1, - Omega_h::Read(mesh.nedges(), Omega_h::I8(1))); - mesh.add_tag( - 1, "class_id", 1, - Omega_h::Read(mesh.nedges(), Omega_h::ClassId(0))); - mesh.add_tag( - 2, "class_dim", 1, - Omega_h::Read(mesh.nelems(), Omega_h::I8(2))); - mesh.add_tag( - 2, "class_id", 1, - Omega_h::Read(mesh.nelems(), Omega_h::ClassId(0))); - return mesh; -} - -std::shared_ptr MakeP1Space(Omega_h::Mesh& mesh) -{ - return pcms::LagrangeFunctionSpace::FromMesh( - mesh, 1, 1, pcms::CoordinateSystem::Cartesian, "global", - pcms::LagrangeFunctionSpace::Backend::OmegaH); -} - -// Evaluates source_field at the integrator's sample points and assembles. -void EvaluateAndAssemble( - pcms::LinearFormIntegrator& integrator, - const std::shared_ptr& source_space, - const pcms::Field& source_field) -{ - const auto& pts = integrator.GetIntegrationPoints(); - const std::size_t npts = pts.GetValues().extent(0); - auto evaluator = source_space->CreatePointEvaluator( - pcms::EvaluationRequest::FromCoordinates(pts)); - Kokkos::View sampled("sampled", npts, - 1); - evaluator->Evaluate(source_field, pcms::MakeRank2View(sampled)); - integrator.Assemble(pcms::MakeConstRank2View(sampled)); -} - -} // namespace - // --------------------------------------------------------------------------- // Monte Carlo RHS integrator // --------------------------------------------------------------------------- @@ -86,20 +23,17 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: sample points lie inside the domain", "[mc_rhs_integrator]") { Omega_h::Library lib; - auto target_mesh = BuildUnitSquare(lib, 0); - auto target_space = MakeP1Space(target_mesh); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto target_space = pcms::test::MakeP1Space(target_mesh); const int samples_per_element = 16; for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { pcms::OmegaHMonteCarloRHSIntegrator integrator( *target_space, samples_per_element, sampling); const auto raw_coords = integrator.GetIntegrationPoints().GetValues(); - auto coords_view = Kokkos::View>( - raw_coords.data_handle(), raw_coords.extent(0), raw_coords.extent(1)); - auto coords_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace{}, coords_view); + auto coords_h = pcms::test::CopyCoordinatesToHost( + raw_coords, static_cast(raw_coords.extent(0)), + static_cast(raw_coords.extent(1))); REQUIRE( coords_h.extent(0) == @@ -120,10 +54,10 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", // For f = c the estimator sums to c * |domain| for any sample placement, // because the P1 basis functions partition unity at every sample point. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( @@ -131,7 +65,7 @@ TEST_CASE("OmegaHMonteCarloRHSIntegrator: constant field integrates exactly", for (const auto sampling : {pcms::MonteCarloSampling::UniformRandom}) { pcms::OmegaHMonteCarloRHSIntegrator integrator(*target_space, 8, sampling); - EvaluateAndAssemble(integrator, source_space, source_field); + pcms::test::EvaluateAndAssemble(integrator, source_space, source_field); PetscScalar sum = 0.0; VecSum(integrator.GetVector(), &sum); @@ -152,10 +86,10 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " // of the source field) equals the source field, so the sampled residual is // identically zero and the projection is exact regardless of sample count. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); auto target_field = target_space->CreateFunction(); @@ -170,11 +104,11 @@ TEST_CASE("OmegaHControlVariateProjection: exact for fields in the target " const auto values = pcms::FlattenToRank1View(target_field.GetDOFHolderDataHost()); - const auto coords_h = Omega_h::HostRead(target_mesh.coords()); + const auto coords_h = pcms::test::CopyCoordinatesToHost( + pcms::MakeConstRank2View(target_mesh.coords(), 2), target_mesh.nverts(), 2); REQUIRE(static_cast(values.size()) == target_mesh.nverts()); for (Omega_h::LO i = 0; i < target_mesh.nverts(); ++i) { - const double expected = - 3.0 * coords_h[2 * i + 0] - coords_h[2 * i + 1] + 0.25; + const double expected = 3.0 * coords_h(i, 0) - coords_h(i, 1) + 0.25; CAPTURE(i, expected, values[i]); CHECK(values[i] == Catch::Approx(expected).margin(1e-9)); } @@ -188,10 +122,10 @@ TEST_CASE("OmegaHControlVariateProjection: reduces error vs plain Monte Carlo", // seeds the control-variate projection must be closer to the exact // (intersection-quadrature) projection than the plain Monte Carlo one. Omega_h::Library lib; - auto source_mesh = BuildUnitSquare(lib, 1); - auto target_mesh = BuildUnitSquare(lib, 0); - auto source_space = MakeP1Space(source_mesh); - auto target_space = MakeP1Space(target_mesh); + auto source_mesh = pcms::test::BuildUnitSquare(lib, 1); + auto target_mesh = pcms::test::BuildUnitSquare(lib, 0); + auto source_space = pcms::test::MakeP1Space(source_mesh); + auto target_space = pcms::test::MakeP1Space(target_mesh); auto source_field = source_space->CreateFunction(); pcms::test::SetField( diff --git a/test/test_point_evaluator.cpp b/test/test_point_evaluator.cpp index dcc57087..e9e235a9 100644 --- a/test/test_point_evaluator.cpp +++ b/test/test_point_evaluator.cpp @@ -34,7 +34,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -44,7 +44,7 @@ TEST_CASE("PointEvaluator: OmegaH order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } // ============================================================================ @@ -67,45 +67,24 @@ TEST_CASE("PointEvaluator: same evaluator reused for two FieldData objects") // field_a: linear_f; field_b: constant 42 pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); - // Create the PointEvaluator once + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - // Evaluate field_a then field_b with the same evaluator - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } // ============================================================================ @@ -125,11 +104,10 @@ TEST_CASE("PointEvaluator: OmegaH order-1 out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); // Points clearly outside [0,1]^2 - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -158,7 +136,7 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -167,7 +145,8 @@ TEST_CASE("PointEvaluator: UniformGrid order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") @@ -184,7 +163,7 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = @@ -193,7 +172,8 @@ TEST_CASE("PointEvaluator: SplineFunctionSpace uniform-grid evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 1e-8); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 1e-8); } // ============================================================================ @@ -308,7 +288,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; @@ -318,7 +298,7 @@ TEST_CASE("PointEvaluator: MeshFields order-1 linear evaluation") pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); pcms::test::CheckEvaluation( *evaluator, field_data, pts, - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); } TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") @@ -334,10 +314,9 @@ TEST_CASE("PointEvaluator: MeshFields out-of-bounds fill") auto field_data = factory->CreateFunction(); pcms::test::SetField( field_data.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); - const std::vector outside_pts = {-0.5, 0.5, 1.5, 0.5, - 0.5, -0.5, 0.5, 1.5}; + const auto outside_pts = pcms::test::StandardOutsideCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView( outside_pts, CoordinateSystem::Cartesian); pcms::OutOfBoundsPolicy policy{pcms::OutOfBoundsMode::FILL, -999.0}; @@ -362,43 +341,24 @@ TEST_CASE( auto field_b = factory->CreateFunction(); pcms::test::SetField( field_a.GetData(), *factory->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); pcms::test::SetField( field_b.GetData(), *factory->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); auto pts = pcms::test::StandardEvalCoords2D(); - int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = factory->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto view_a = pcms::Rank2View( - out_a_device.data(), n, 1); - auto view_b = pcms::Rank2View( - out_b_device.data(), n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(pcms::HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(1e-10)); - REQUIRE(out_b_host(i) == Catch::Approx(42.0).margin(1e-10)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return Real(42); }); } TEST_CASE("LagrangeFunctionSpace: MeshFields rejects multi-component fields") diff --git a/test/test_polynomial_reconstruction_mls_evaluation.cpp b/test/test_polynomial_reconstruction_mls_evaluation.cpp index 888a83b1..415d4968 100644 --- a/test/test_polynomial_reconstruction_mls_evaluation.cpp +++ b/test/test_polynomial_reconstruction_mls_evaluation.cpp @@ -77,12 +77,6 @@ pcms::MLSOptions DefaultTestOptions3D() return opts; } -// Interior query points for a unit box — same as StandardEvalCoords2D. -std::vector QueryPoints() -{ - return pcms::test::StandardEvalCoords2D(); -} - pcms::MLSOptions SweepTestOptions(unsigned degree, pcms::RadialBasisFunction basis) { @@ -122,7 +116,7 @@ void CheckPolynomialReproduction(unsigned degree, auto field = fs->CreateFunction(); pcms::test::SetField(field.GetData(), *fs->GetLayout(), func); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( @@ -149,7 +143,9 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: reproduces " CheckPolynomialReproduction( 0, basis, OMEGA_H_LAMBDA(Real, Real) { return 3.14; }, 5e-3); CheckPolynomialReproduction( - 1, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }, 5e-3); + 1, basis, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); CheckPolynomialReproduction( 2, basis, OMEGA_H_LAMBDA(Real x, Real y) { return x * x + x * y + 2.0 * y * y; }, @@ -175,43 +171,26 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: same PointEvaluator " pcms::test::SetField( field_a.GetData(), *fs->GetLayout(), - OMEGA_H_LAMBDA(Real x, Real y) { return x + 2.0 * y; }); + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }); const Real cval = 7.0; pcms::test::SetField( field_b.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return cval; }); - auto pts = QueryPoints(); - int n = static_cast(pts.size()) / 2; + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); + // Create the PointEvaluator once and reuse it for both fields. auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_a_device("out_a", n); - Kokkos::View out_b_device("out_b", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View view_a(out_a_device.data(), - n, 1); - Rank2View view_b(out_b_device.data(), - n, 1); - - evaluator->Evaluate(field_a, view_a); - evaluator->Evaluate(field_b, view_b); - - auto out_a_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_a_device); - auto out_b_host = - Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_b_device); - - for (int i = 0; i < n; ++i) { - Real x = pts[2 * static_cast(i)], - y = pts[2 * static_cast(i) + 1]; - REQUIRE(out_a_host(i) == - Catch::Approx(pcms::test::linear_f(x, y)).margin(5e-3)); - REQUIRE(out_b_host(i) == Catch::Approx(cval).margin(5e-3)); - } + pcms::test::CheckEvaluation( + *evaluator, field_a, pts, + OMEGA_H_LAMBDA(Real x, Real y) { return pcms::test::linear_f(x, y); }, + 5e-3); + pcms::test::CheckEvaluation( + *evaluator, field_b, pts, OMEGA_H_LAMBDA(Real, Real) { return cval; }, + 5e-3); } // ============================================================================ @@ -228,7 +207,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " coords_view, CoordinateSystem::Cartesian); auto field = fs->CreateFunction(); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); @@ -236,13 +215,8 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: Evaluate throws for " pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); // Two-component output — must throw - Kokkos::View out_device("out", - static_cast(n) * 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 2); - REQUIRE_THROWS(evaluator->Evaluate(field, out_view)); + Kokkos::View out_device("out", n, 2); + REQUIRE_THROWS(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); } // ============================================================================ @@ -263,24 +237,20 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: default MLSOptions — " field.GetData(), *fs->GetLayout(), OMEGA_H_LAMBDA(Real, Real) { return Real(1.0); }); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); int n = static_cast(pts.size()) / 2; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", n); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - auto out_view = - Rank2View(out_device.data(), n, 1); + Kokkos::View out_device("out", n, 1); // Just verify it runs without error and returns finite values - REQUIRE_NOTHROW(evaluator->Evaluate(field, out_view)); + REQUIRE_NOTHROW(evaluator->Evaluate(field, pcms::MakeRank2View(out_device))); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); for (int i = 0; i < n; ++i) - REQUIRE(std::isfinite(out_host(i))); + REQUIRE(std::isfinite(out_host(i, 0))); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " @@ -293,7 +263,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -310,7 +280,7 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: CreatePointEvaluator " auto fs = pcms::PolynomialReconstructionFunctionSpace::Create( coords_view, CoordinateSystem::Cylindrical, DefaultTestOptions()); - auto pts = QueryPoints(); + auto pts = pcms::test::StandardEvalCoords2D(); auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cylindrical); REQUIRE_THROWS(fs->CreatePointEvaluator( @@ -342,16 +312,12 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: radius option is " pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 1); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 1, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 1, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); - REQUIRE(out_host(0) == Catch::Approx(1.0).margin(1e-8)); + REQUIRE(out_host(0, 0) == Catch::Approx(1.0).margin(1e-8)); } TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " @@ -375,34 +341,24 @@ TEST_CASE("PolynomialReconstructionFunctionSpace MLS: 3D point clouds preserve " } auto field = fs->CreateFunction(); - std::vector dof_values(27); - for (int i = 0; i < 27; ++i) { - const Real x = src[3 * i + 0]; - const Real y = src[3 * i + 1]; - const Real z = src[3 * i + 2]; - dof_values[i] = x + 2.0 * y + 3.0 * z; - } - field.GetData().SetDOFHolderDataHost(Rank2View( - dof_values.data(), static_cast(dof_values.size()), 1)); + pcms::test::SetField( + field, + KOKKOS_LAMBDA(Real x, Real y, Real z) { return x + 2.0 * y + 3.0 * z; }); std::vector pts{0.5, 0.5, 0.25, 0.5, 0.5, 0.75}; auto device_coords = pcms::test::CreateDeviceCoordinateView(pts, CoordinateSystem::Cartesian, 3); auto evaluator = fs->CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View out_device("out", 2); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - Rank2View out_view(out_device.data(), - 2, 1); - evaluator->Evaluate(field, out_view); + Kokkos::View out_device("out", 2, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(out_device)); auto out_host = Kokkos::create_mirror_view_and_copy(HostMemorySpace(), out_device); // This 3D case is a small, low-resolution support cloud with MLS weights // built from a finite-radius neighborhood rather than exact nodal lookup, so // it is checked with a slightly looser tolerance than the denser 2D tests. - REQUIRE(out_host(0) == Catch::Approx(2.25).margin(1e-2)); - REQUIRE(out_host(1) == Catch::Approx(3.75).margin(1e-2)); - REQUIRE(out_host(1) - out_host(0) == Catch::Approx(1.5).margin(1e-2)); + REQUIRE(out_host(0, 0) == Catch::Approx(2.25).margin(1e-2)); + REQUIRE(out_host(1, 0) == Catch::Approx(3.75).margin(1e-2)); + REQUIRE(out_host(1, 0) - out_host(0, 0) == Catch::Approx(1.5).margin(1e-2)); } diff --git a/test/test_spr_meshfields.cpp b/test/test_spr_meshfields.cpp index 42a6acf0..7ad4470a 100644 --- a/test/test_spr_meshfields.cpp +++ b/test/test_spr_meshfields.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -127,22 +128,7 @@ TEST_CASE("meshfields_spr_test") const auto& ntargets = mesh.nverts(); - Write source_coordinates( - dim * nfaces, 0, "stores coordinates of cell centroid of each tri element"); - - const auto& faces2nodes = mesh.ask_down(FACE, VERT).ab2b; - - Kokkos::parallel_for( - "calculate the centroid in each tri element", nfaces, - OMEGA_H_LAMBDA(const LO id) { - const auto current_el_verts = gather_verts<3>(faces2nodes, id); - const Omega_h::Few, 3> current_el_vert_coords = - gather_vectors<3, 2>(target_coordinates, current_el_verts); - auto centroid = average(current_el_vert_coords); - int index = 2 * id; - source_coordinates[index] = centroid[0]; - source_coordinates[index + 1] = centroid[1]; - }); + const auto source_coordinates = pcms::get_entity_centroids(mesh, FACE); pcms::Points source_points; diff --git a/test/test_svd_serial.cpp b/test/test_svd_serial.cpp index bdc23290..1774e29b 100644 --- a/test/test_svd_serial.cpp +++ b/test/test_svd_serial.cpp @@ -38,7 +38,7 @@ TEST_CASE("test_serial_svd") Kokkos::deep_copy(A_data, host_A_data); - Kokkos::View rhs_data("Device rhs data", row, column); + Kokkos::View rhs_data("Device rhs data", row); auto host_rhs_data = Kokkos::create_mirror_view(rhs_data); host_rhs_data(0) = 1.28571; @@ -54,7 +54,7 @@ TEST_CASE("test_serial_svd") { Kokkos::View result("result", row, row); - Kokkos::View transpose_expected("result", row, row); + Kokkos::View transpose_expected("transpose_expected", row, row); Kokkos::deep_copy(result, 0.0); Kokkos::deep_copy(transpose_expected, 0.0); team_policy tp(1, Kokkos::AUTO); diff --git a/test/test_uniform_grid_field.cpp b/test/test_uniform_grid_field.cpp index 632bd4db..699ec26a 100644 --- a/test/test_uniform_grid_field.cpp +++ b/test/test_uniform_grid_field.cpp @@ -72,7 +72,7 @@ void VerifyUniformGridFieldValues( int vertex_id = j * (grid.divisions[0] + 1) + i; pcms::Real x = ug_coords.GetValues()(vertex_id, 0); pcms::Real y = ug_coords.GetValues()(vertex_id, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); pcms::Real actual = ug_field_data[vertex_id]; REQUIRE(std::abs(expected - actual) <= 1e-10); } @@ -151,21 +151,16 @@ TEST_CASE("UniformGrid order-0 field creation and evaluation") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); - - REQUIRE(results_host(0) == Catch::Approx(1.0)); - REQUIRE(results_host(1) == Catch::Approx(2.0)); - REQUIRE(results_host(2) == Catch::Approx(3.0)); - REQUIRE(results_host(3) == Catch::Approx(4.0)); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); + + REQUIRE(results_host(0, 0) == Catch::Approx(1.0)); + REQUIRE(results_host(1, 0) == Catch::Approx(2.0)); + REQUIRE(results_host(2, 0) == Catch::Approx(3.0)); + REQUIRE(results_host(3, 0) == Catch::Approx(4.0)); } TEST_CASE("UniformGrid field data operations", "[uniform_grid_field]") @@ -236,26 +231,21 @@ TEST_CASE("UniformGrid field evaluation - piecewise constant") auto evaluator = eval_factory.CreatePointEvaluator( pcms::EvaluationRequest::FromCoordinates(device_coords.coordinate_view)); - Kokkos::View results_host("results_host", - 4); - Kokkos::View results_device( - "results_device", 4); - using LayoutPolicy = - pcms::detail::default_layout_for_memory_space_t; - pcms::Rank2View out( - results_device.data(), 4, 1); - evaluator->Evaluate(field, out); - Kokkos::deep_copy(results_host, results_device); + Kokkos::View results_device( + "results_device", 4, 1); + evaluator->Evaluate(field, pcms::MakeRank2View(results_device)); + auto results_host = Kokkos::create_mirror_view_and_copy( + pcms::HostMemorySpace(), results_device); // Check results - interpolated from vertices // Cell 0 center (2.5, 2.5): avg of v0,v1,v3,v4 = (1.0+1.5+2.0+2.5)/4 = 1.75 // Cell 1 center (7.5, 2.5): avg of v1,v2,v4,v5 = (1.5+2.0+2.5+3.0)/4 = 2.25 // Cell 2 center (2.5, 7.5): avg of v3,v4,v6,v7 = (2.0+2.5+3.0+3.5)/4 = 2.75 // Cell 3 center (7.5, 7.5): avg of v4,v5,v7,v8 = (2.5+3.0+3.5+4.0)/4 = 3.25 - REQUIRE(std::abs(results_host(0) - 1.75) < 1e-10); - REQUIRE(std::abs(results_host(1) - 2.25) < 1e-10); - REQUIRE(std::abs(results_host(2) - 2.75) < 1e-10); - REQUIRE(std::abs(results_host(3) - 3.25) < 1e-10); + REQUIRE(std::abs(results_host(0, 0) - 1.75) < 1e-10); + REQUIRE(std::abs(results_host(1, 0) - 2.25) < 1e-10); + REQUIRE(std::abs(results_host(2, 0) - 2.75) < 1e-10); + REQUIRE(std::abs(results_host(3, 0) - 3.25) < 1e-10); } TEST_CASE("UniformGrid field serialization") @@ -328,8 +318,9 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); pcms::UniformGrid<2> grid; grid.edge_length = {1.0, 1.0}; @@ -354,7 +345,7 @@ TEST_CASE("Transfer from OmegaH field to UniformGrid field") for (int i = 0; i < num_ug_nodes; ++i) { pcms::Real x = ug_coords_host(i, 0); pcms::Real y = ug_coords_host(i, 1); - pcms::Real expected = x + 2.0 * y; + pcms::Real expected = pcms::test::linear_f(x, y); REQUIRE(std::abs(transferred_data[i] - expected) < 1e-6); } } @@ -641,8 +632,9 @@ TEST_CASE("UniformGrid workflow") mesh, 1, 1, pcms::CoordinateSystem::Cartesian); auto omega_h_field = omega_h_factory->CreateFunction(); pcms::test::SetField( - omega_h_field, - OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { return x + 2.0 * y; }); + omega_h_field, OMEGA_H_LAMBDA(pcms::Real x, pcms::Real y) { + return pcms::test::linear_f(x, y); + }); auto ug_factory = pcms::LagrangeFunctionSpace::FromUniformGrid( grid, 1, pcms::CoordinateSystem::Cartesian); @@ -658,23 +650,11 @@ TEST_CASE("UniformGrid workflow") auto ug_coords_host_view = pcms::test::CopyCoordinatesToHost(ug_coords_device_view, 25, 2); - auto ug_field_data_device = ug_field.GetDOFHolderData(); - Kokkos::View ug_field_data_device_view( - "", 25); - Kokkos::parallel_for( - "CopyFieldDataToView", 25, KOKKOS_LAMBDA(int i) { - ug_field_data_device_view(i) = ug_field_data_device(i, 0); - }); - auto ug_field_data_host_view = - Kokkos::View("", 25); - Kokkos::deep_copy(ug_field_data_host_view, ug_field_data_device_view); - - pcms::Rank2View ug_coords( - ug_coords_host_view.data(), 25, 2); pcms::CoordinateView ug_coords_view( - pcms::CoordinateSystem::Cartesian, ug_coords); - pcms::Rank1View ug_field_data( - ug_field_data_host_view.data(), 25); + pcms::CoordinateSystem::Cartesian, + pcms::MakeConstRank2View(ug_coords_host_view)); + const auto ug_field_data = + pcms::FlattenToRank1View(ug_field.GetDOFHolderDataHost()); VerifyUniformGridFieldValues(grid, ug_coords_view, ug_field_data); VerifyMaskFieldValues(grid, mask_field); diff --git a/test/test_xgc_reverse_classification.cpp b/test/test_xgc_reverse_classification.cpp index 4b87fcb9..5603b52b 100644 --- a/test/test_xgc_reverse_classification.cpp +++ b/test/test_xgc_reverse_classification.cpp @@ -39,7 +39,7 @@ TEST_CASE("reverse classification") } auto vec = rc.Serialize(); pcms::ReverseClassificationVertex rc_deserialized; - pcms::Rank1View av{vec.data(), vec.size()}; + auto av = pcms::make_array_view(vec); rc_deserialized.Deserialize(av); REQUIRE(rc_deserialized == rc); }