From 4cb57af9908babeaa54781f3777ccc241a2b0667 Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Wed, 23 Oct 2024 17:00:53 -0400 Subject: [PATCH 1/2] add test for find_exit_face in pumipic_adjacency.tpp --- test/CMakeLists.txt | 1 + test/test_find_exit_face.cpp | 233 +++++++++++++++++++++++++++++++++++ test/testing.cmake | 5 + 3 files changed, 239 insertions(+) create mode 100644 test/test_find_exit_face.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 94a88763..a918838f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -22,6 +22,7 @@ make_test(search2d search2d.cpp) make_test(pseudoXGCm pseudoXGCm.cpp) make_test(pseudoXGCm_scatter pseudoXGCm_scatter.cpp) make_test(loadSerialMesh loadSerialMesh.cpp) +make_test(test_find_exit_face test_find_exit_face.cpp) include(testing.cmake) bob_end_subdir() diff --git a/test/test_find_exit_face.cpp b/test/test_find_exit_face.cpp new file mode 100644 index 00000000..2f82651f --- /dev/null +++ b/test/test_find_exit_face.cpp @@ -0,0 +1,233 @@ +/* + * Test for find_exit_face() function in pumipic_adjacency.tpp + * This test is for 2D mesh + * It shows that when the exis face is computed using `useBcc` flag on, it + fails to find the correct exit face + + * The test mesh is a square as below with the following edges: + + 1(-5,5) e3 2(5,5) + *-----------------------------------* + | \ o / | + | \ / | + | \ e4 e6 / | + | \ / | + | \ / | +e0 | * 4(0,0) | e5 + | / \ | + | e2 / \ e7 | + | / \ | + | / \ | + | / y x \ | + *-----------------------------------* + 0(-5,-5) e1 3(5,-5) + + * Two rays are cast from the point x->o and y->o. + * The expected exit faces are 7 and 2 respectively + + x(2.520540, -3.193840) y(-3.03604, -3.193840) + o(-0.294799, 3.072990) + + TODO It fails for the x->o ray and returns edge 2 instead of 7 +*/ + +#include "pumipic_adjacency.hpp" +#include "pumipic_kktypes.hpp" +#include "pumipic_library.hpp" +#include "pumipic_mesh.hpp" +#include +#include +#include +#include +#include +#include + +using particle_structs::lid_t; +using particle_structs::MemberTypes; +using particle_structs::SellCSigma; +using pumipic::fp_t; +using pumipic::Vector3d; + +namespace o = Omega_h; +namespace p = pumipic; +namespace ps = particle_structs; + +// To demonstrate push and adjacency search we store: +//-two fp_t[3] arrays, 'Vector3d', for the current and +// computed (pre adjacency search) positions, and +//-an integer to store the particles id +typedef MemberTypes Particle; +typedef ps::ParticleStructure PS; +typedef Kokkos::DefaultExecutionSpace ExeSpace; + +void print_2D_mesh_edges(o::Mesh &mesh) { + printf("------------------ Mesh Info ------------------\n"); + const auto edge2node = mesh.ask_down(o::EDGE, o::VERT).ab2b; + const auto coords = mesh.coords(); + // print edge nodes and their coordinates given an edge id + auto print_edge_nodes = OMEGA_H_LAMBDA(o::LO edge_id) { + + printf("Edge %d nodes (%d, %d): ", edge_id, edge2node[2 * edge_id], + edge2node[2 * edge_id + 1]); + for (int i = 0; i < 2; i++) { + printf("%d (%f, %f) ", edge2node[2 * edge_id + i], + coords[edge2node[2 * edge_id + i] * 2 + 0], + coords[edge2node[2 * edge_id + i] * 2 + 1]); + } + printf("\n"); + }; + o::parallel_for(mesh.nedges(), print_edge_nodes); + // add a barrier to make sure all the prints are done + printf("------------------ Mesh Info ------------------\n"); +} + +bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { + printf("Test: 2D intersection...\n"); + o::Library &olib = lib->omega_h_lib(); + // load mesh + o::Mesh mesh = Omega_h::gmsh::read(mesh_fname, olib.self()); + printf("Mesh loaded successfully with %d elements\n", mesh.nelems()); + // o::vtk::write_parallel("square4el.vtk", &mesh); + + print_2D_mesh_edges(mesh); + + Omega_h::Write owners(mesh.nelems(), 0); + p::Mesh picparts(mesh, owners); + o::Mesh *p_mesh = picparts.mesh(); + p_mesh->ask_elem_verts(); + + // create particles + Kokkos::TeamPolicy policy; +#ifdef PP_USE_GPU + printf("Using GPU for simulation..."); + policy = + Kokkos::TeamPolicy(10, Kokkos::AUTO()); +#else + printf("Using CPU for simulation..."); + policy = Kokkos::TeamPolicy(100, 1); +#endif + + o::Int ne = p_mesh->nelems(); + PS::kkLidView ptcls_per_elem("ptcls_per_elem", ne); + PS::kkGidView element_gids("element_gids", ne); + o::GOs mesh_element_gids = picparts.globalIds(picparts.dim()); + o::parallel_for( + ne, OMEGA_H_LAMBDA(const int &i) { + element_gids(i) = mesh_element_gids[i]; + if (i == 1) { + ptcls_per_elem(i) = 2; + } else { + ptcls_per_elem(i) = 0; + } + }); + +#ifdef PP_ENABLE_CAB + PS *ptcls = new p::DPS(policy, ne, 2, ptcls_per_elem, element_gids); + printf("DPS Particle structure created successfully\n"); +#else + PS *ptcls = PS *ptcls = new SellCSigma( + policy, 10, 10, ne, 2, ptcls_per_elem, element_gids); + printf("SellCSigma Particle structure created successfully\n"); +#endif + + // set particle position + auto particle_init_position = ptcls->get<0>(); + auto particle_final_position = ptcls->get<1>(); + auto pid_d = ptcls->get<2>(); + auto setIDs = PS_LAMBDA(const int &eid, const int &pid, const bool &mask) { + pid_d(pid) = pid; + }; + ps::parallel_for(ptcls, setIDs); + + // const o::Vector<3> start_point {0.0, 2.0, -3.0}; + // const o::Vector<3> end_point {0.0, -4.0, -3.0}; + const o::Vector<3> start_point{2.52054, -3.19384, 0}; + const o::Vector<3> start_point2{-3.03604, -3.193840, 0}; + const o::Vector<3> end_point{-0.294799, 3.07299, 0}; + + printf("Start point y: %f %f %f\n", start_point[0], start_point[1], + start_point[2]); + printf("Start point x: %f %f %f\n", start_point2[0], start_point2[1], + start_point2[2]); + printf("End point: %f %f %f\n", end_point[0], end_point[1], end_point[2]); + + auto set_initial_and_final_position = + PS_LAMBDA(const int &e, const int &pid, const int &mask) { + // see test calculation notebook + if (pid == 0) { + particle_init_position(pid, 0) = start_point[0]; + particle_init_position(pid, 1) = start_point[1]; + particle_init_position(pid, 2) = start_point[2]; + + particle_final_position(pid, 0) = end_point[0]; + particle_final_position(pid, 1) = end_point[1]; + particle_final_position(pid, 2) = end_point[2]; + } + if (pid == 1) { + particle_init_position(pid, 0) = start_point2[0]; + particle_init_position(pid, 1) = start_point2[1]; + particle_init_position(pid, 2) = start_point2[2]; + + particle_final_position(pid, 0) = end_point[0]; + particle_final_position(pid, 1) = end_point[1]; + particle_final_position(pid, 2) = end_point[2]; + } + }; + + ps::parallel_for(ptcls, set_initial_and_final_position, + "set_initial_and_final_position"); + + // search for intersection + const auto psCapacity = ptcls->capacity(); + o::Write elem_ids; + o::Write xpoints_d(3 * psCapacity, 0.0, "intersection points"); + o::Write xface_id(psCapacity, 0.0, "intersection faces"); + + const auto elmArea = measure_elements_real(&mesh); + bool useBcc = true; + o::Real tol = p::compute_tolerance_from_area(elmArea); + auto inter_points = + o::Write(2 * psCapacity, 0.0, "intersection points"); + o::Write lastExit(psCapacity, -1, "search_last_exit"); + o::Write ptcl_done(psCapacity, 0, "search_ptcl_done"); + + o::Write elem_ids_2d(psCapacity, 1, "elem_ids for find_exit_face"); + + p::find_exit_face(*p_mesh, ptcls, particle_init_position, + particle_final_position, elem_ids_2d, ptcl_done, elmArea, + useBcc, lastExit, inter_points, tol); + + printf("Last Exit Face: "); + + auto lastExit_h = o::HostRead(lastExit); + bool passed_right = (lastExit_h[0] == 7); + bool passed_left = (lastExit_h[1] == 2); + + if (!passed_right) { + printf("Expected exit face: 7, got %d\n", lastExit_h[0]); + } + if (!passed_left) { + printf("Expected exit face: 2, got %d\n", lastExit_h[1]); + } + + // clean up + delete ptcls; + + return passed_right && passed_left; +} + +int main(int argc, char **argv) { + p::Library lib(&argc, &argv); + if (argc != 2) { + printf("Usage: %s \n", argv[0]); + exit(1); + } + std::string mesh_fname = argv[1]; + bool passed = test_2D_intersection(mesh_fname, &lib); + if (!passed) { + printf("Test failed. Run in verbose mode to view details.\n"); + return 1; + } + + return 0; +} \ No newline at end of file diff --git a/test/testing.cmake b/test/testing.cmake index 48e813e5..3b6580a2 100644 --- a/test/testing.cmake +++ b/test/testing.cmake @@ -1,3 +1,8 @@ +# unit tests +mpi_test(test_find_exit_face 1 + ./test_find_exit_face + ${TEST_DATA_DIR}/square2d4elem.msh) + #simple tests mpi_test(barycentric_3 1 ./barycentric test1) From afc8c745820a60ec4360f56229e34c2ae72993e8 Mon Sep 17 00:00:00 2001 From: Fuad Hasibul Hasan Date: Thu, 24 Oct 2024 14:49:18 -0400 Subject: [PATCH 2/2] test find_exit_face without bcc based search --- test/test_find_exit_face.cpp | 98 +++++++++++++++++++++++++++--------- 1 file changed, 74 insertions(+), 24 deletions(-) diff --git a/test/test_find_exit_face.cpp b/test/test_find_exit_face.cpp index 2f82651f..17385514 100644 --- a/test/test_find_exit_face.cpp +++ b/test/test_find_exit_face.cpp @@ -60,6 +60,37 @@ typedef MemberTypes Particle; typedef ps::ParticleStructure PS; typedef Kokkos::DefaultExecutionSpace ExeSpace; +OMEGA_H_DEVICE bool all_positive(const o::Vector<3>& a, Omega_h::Real tol = EPSILON) { + auto isPos = 1; + for (o::LO i = 0; i < a.size(); ++i) { + const auto gtez = o::are_close(a[i], 0.0, tol, tol) || a[i] > 0; + isPos = isPos && gtez; + } + return isPos; +} + +bool is_inside(o::Mesh &mesh, o::LO elem_id, const o::Vector<2> point) { + const auto &coords = mesh.coords(); + const auto &faces2nodes = mesh.ask_verts_of(o::FACE); + + o::Write inside(2, 0); + + auto is_inside_lambda = OMEGA_H_LAMBDA(o::LO id){ + const auto current_el_verts = o::gather_verts<3>(faces2nodes, elem_id); + const Omega_h::Few, 3> current_el_vert_coords = + o::gather_vectors<3, 2>(coords, current_el_verts); + // check if the particle is in the element + o::Vector<3> bcc = o::barycentric_from_global<2, 2>( + point, current_el_vert_coords); + // if all bcc are positive, the point is inside the element + inside[0] = all_positive(bcc,0.0); + }; + o::parallel_for(1, is_inside_lambda); + auto host_inside = o::HostWrite(inside); + + return bool(host_inside[0]); +} + void print_2D_mesh_edges(o::Mesh &mesh) { printf("------------------ Mesh Info ------------------\n"); const auto edge2node = mesh.ask_down(o::EDGE, o::VERT).ab2b; @@ -81,7 +112,7 @@ void print_2D_mesh_edges(o::Mesh &mesh) { printf("------------------ Mesh Info ------------------\n"); } -bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { +bool test_2D_intersection(const std::string mesh_fname, p::Library *lib, bool useBcc) { printf("Test: 2D intersection...\n"); o::Library &olib = lib->omega_h_lib(); // load mesh @@ -141,15 +172,27 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { // const o::Vector<3> start_point {0.0, 2.0, -3.0}; // const o::Vector<3> end_point {0.0, -4.0, -3.0}; - const o::Vector<3> start_point{2.52054, -3.19384, 0}; - const o::Vector<3> start_point2{-3.03604, -3.193840, 0}; - const o::Vector<3> end_point{-0.294799, 3.07299, 0}; + const o::Vector<2> start_point{2.52054, -3.19384}; + const o::Vector<2> start_point2{-3.03604, -3.193840}; + const o::Vector<2> end_point{-0.294799, 3.07299}; + if (!is_inside(mesh, 1, start_point)) { + printf("Error: Start point is not inside the expected element.\n"); + exit(1); + } + if (!is_inside(mesh, 1, start_point2)) { + printf("Error: Start point 1 is not inside the expected element.\n"); + exit(1); + } + if (is_inside(mesh, 1, end_point)) { + printf("Error: End point is incorrectly inside the expected element.\n"); + exit(1); + } - printf("Start point y: %f %f %f\n", start_point[0], start_point[1], - start_point[2]); - printf("Start point x: %f %f %f\n", start_point2[0], start_point2[1], - start_point2[2]); - printf("End point: %f %f %f\n", end_point[0], end_point[1], end_point[2]); + printf("Start point y: %f %f\n", start_point[0], + start_point[1]); + printf("Start point x: %f %f\n", start_point2[0], + start_point2[1]); + printf("End point: %f %f\n", end_point[0], end_point[1]); auto set_initial_and_final_position = PS_LAMBDA(const int &e, const int &pid, const int &mask) { @@ -157,20 +200,20 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { if (pid == 0) { particle_init_position(pid, 0) = start_point[0]; particle_init_position(pid, 1) = start_point[1]; - particle_init_position(pid, 2) = start_point[2]; + particle_init_position(pid, 2) = 0.0; particle_final_position(pid, 0) = end_point[0]; particle_final_position(pid, 1) = end_point[1]; - particle_final_position(pid, 2) = end_point[2]; + particle_final_position(pid, 2) = 0.0; } if (pid == 1) { particle_init_position(pid, 0) = start_point2[0]; particle_init_position(pid, 1) = start_point2[1]; - particle_init_position(pid, 2) = start_point2[2]; + particle_init_position(pid, 2) = 0.0; particle_final_position(pid, 0) = end_point[0]; particle_final_position(pid, 1) = end_point[1]; - particle_final_position(pid, 2) = end_point[2]; + particle_final_position(pid, 2) = 0.0; } }; @@ -179,12 +222,10 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { // search for intersection const auto psCapacity = ptcls->capacity(); - o::Write elem_ids; o::Write xpoints_d(3 * psCapacity, 0.0, "intersection points"); - o::Write xface_id(psCapacity, 0.0, "intersection faces"); + o::Write xface_id(psCapacity, -1, "intersection faces"); const auto elmArea = measure_elements_real(&mesh); - bool useBcc = true; o::Real tol = p::compute_tolerance_from_area(elmArea); auto inter_points = o::Write(2 * psCapacity, 0.0, "intersection points"); @@ -197,17 +238,15 @@ bool test_2D_intersection(const std::string mesh_fname, p::Library *lib) { particle_final_position, elem_ids_2d, ptcl_done, elmArea, useBcc, lastExit, inter_points, tol); - printf("Last Exit Face: "); - auto lastExit_h = o::HostRead(lastExit); bool passed_right = (lastExit_h[0] == 7); bool passed_left = (lastExit_h[1] == 2); if (!passed_right) { - printf("Expected exit face: 7, got %d\n", lastExit_h[0]); + printf("ERROR: Expected exit face: 7, got %d\n", lastExit_h[0]); } if (!passed_left) { - printf("Expected exit face: 2, got %d\n", lastExit_h[1]); + printf("ERROR: Expected exit face: 2, got %d\n", lastExit_h[1]); } // clean up @@ -223,11 +262,22 @@ int main(int argc, char **argv) { exit(1); } std::string mesh_fname = argv[1]; - bool passed = test_2D_intersection(mesh_fname, &lib); - if (!passed) { - printf("Test failed. Run in verbose mode to view details.\n"); + + bool passed_bcc = test_2D_intersection(mesh_fname, &lib, true); + if (!passed_bcc) { + printf("ERROR: Test failed with BCC. Run in verbose mode to view details.\n"); + } + + bool passed_woBcc = test_2D_intersection(mesh_fname, &lib, false); + if (!passed_woBcc) { + printf("ERROR: Test failed without BCC. Run in verbose mode to view details.\n"); + } + + if (passed_bcc && passed_woBcc) { + printf("All tests passed!\n"); + return 0; + } else { return 1; } - return 0; } \ No newline at end of file