diff --git a/src/ipc/broad_phase/cuda/lbvh.cu b/src/ipc/broad_phase/cuda/lbvh.cu index b5b62f340..c63136c09 100644 --- a/src/ipc/broad_phase/cuda/lbvh.cu +++ b/src/ipc/broad_phase/cuda/lbvh.cu @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -354,25 +355,32 @@ namespace { } const size_t num_nodes = size_t(2) * n - 1; - bvh.nodes.resize(num_nodes); - bvh.rightmost_leaves.resize(num_nodes); + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("resize_bvh"); + bvh.nodes.resize(num_nodes); + bvh.rightmost_leaves.resize(num_nodes); - for (auto& codes : impl.morton_codes) { - codes.resize(n); + for (auto& codes : impl.morton_codes) { + codes.resize(n); + } + for (auto& ids : impl.box_ids) { + ids.resize(n); + } + // Only the visitation counts need zeroing; a memset is the + // cheapest way to do it. + impl.construction_infos.resize(num_nodes); + impl.construction_infos.zero(); + IPC_TOOLKIT_CUDA_CHECK(cudaMemsetAsync(d_root, 0xFF, sizeof(int))); } - for (auto& ids : impl.box_ids) { - ids.resize(n); + + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("compute_morton_codes"); + compute_morton_codes_kernel<<< + kernel_grid_size(n), KERNEL_BLOCK_SIZE>>>( + d_box_min, d_box_max, n, impl.domain.data(), dim, + impl.morton_codes[0].data(), impl.box_ids[0].data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); } - // Only the visitation counts need zeroing; a memset is the cheapest - // way to do it. - impl.construction_infos.resize(num_nodes); - impl.construction_infos.zero(); - IPC_TOOLKIT_CUDA_CHECK(cudaMemsetAsync(d_root, 0xFF, sizeof(int))); - - compute_morton_codes_kernel<<>>( - d_box_min, d_box_max, n, impl.domain.data(), dim, - impl.morton_codes[0].data(), impl.box_ids[0].data()); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); // Radix sort the (code, id) pairs by code. // @@ -392,29 +400,39 @@ namespace { // Keeping that scratch in a persistent buffer is what makes this // allocation-free per build (thrust::sort_by_key would cudaMalloc and // cudaFree it every call). - size_t temp_bytes = 0; - IPC_TOOLKIT_CUDA_CHECK( - cub::DeviceRadixSort::SortPairs( - nullptr, temp_bytes, keys, values, n)); - impl.sort_temp.resize(temp_bytes); - IPC_TOOLKIT_CUDA_CHECK( - cub::DeviceRadixSort::SortPairs( - impl.sort_temp.data(), temp_bytes, keys, values, n)); - - // Current() is the sorted half of each ping-pong pair. - build_hierarchy_kernel<<>>( - d_box_min, d_box_max, keys.Current(), values.Current(), n, - bvh.nodes.data(), bvh.rightmost_leaves.data(), - impl.construction_infos.data(), d_root); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("sort_morton_codes"); + size_t temp_bytes = 0; + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceRadixSort::SortPairs( + nullptr, temp_bytes, keys, values, n)); + impl.sort_temp.resize(temp_bytes); + IPC_TOOLKIT_CUDA_CHECK( + cub::DeviceRadixSort::SortPairs( + impl.sort_temp.data(), temp_bytes, keys, values, n)); + } - swap_root_kernel<<<1, 1>>>( - bvh.nodes.data(), bvh.rightmost_leaves.data(), d_root); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("build_hierarchy"); + // Current() is the sorted half of each ping-pong pair. + build_hierarchy_kernel<<>>( + d_box_min, d_box_max, keys.Current(), values.Current(), n, + bvh.nodes.data(), bvh.rightmost_leaves.data(), + impl.construction_infos.data(), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } - patch_left_kernel<<>>( - bvh.nodes.data(), static_cast(num_nodes), d_root); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("patch_root"); + swap_root_kernel<<<1, 1>>>( + bvh.nodes.data(), bvh.rightmost_leaves.data(), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + + patch_left_kernel<<< + kernel_grid_size(num_nodes), KERNEL_BLOCK_SIZE>>>( + bvh.nodes.data(), static_cast(num_nodes), d_root); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } } /// @brief Reduce the Morton-normalization domain (min of mins, max of @@ -425,6 +443,8 @@ namespace { /// reads it directly, with no host round-trip. void compute_domain(LBVH::Impl& impl, const int n_vertices) { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("compute_domain"); + Domain init; for (int k = 0; k < 3; ++k) { init.min[k] = std::numeric_limits::max(); @@ -453,6 +473,8 @@ namespace { void upload_vertices( Eigen::ConstRef vertices, DeviceBuffer& d) { + IPC_TOOLKIT_PROFILE_BLOCK("upload_vertices"); + const size_t n = static_cast(vertices.size()); if (vertices.innerStride() == 1 && vertices.outerStride() == vertices.rows()) { @@ -471,6 +493,8 @@ namespace { std::vector& h, DeviceBuffer& d) { + IPC_TOOLKIT_PROFILE_BLOCK("upload_connectivity"); + const size_t n = M.rows(); h.resize(Cols * n); for (size_t i = 0; i < n; ++i) { @@ -518,24 +542,29 @@ namespace { upload_connectivity<3>(faces, impl.h_faces, impl.faces); // Build edge/face boxes on the device from the vertex boxes. - impl.ebox_min.resize(3 * size_t(n_edges)); - impl.ebox_max.resize(3 * size_t(n_edges)); - if (n_edges > 0) { - build_edge_boxes_kernel<<< - kernel_grid_size(n_edges), KERNEL_BLOCK_SIZE>>>( - impl.vbox_min.data(), impl.vbox_max.data(), impl.edges.data(), - n_edges, impl.ebox_min.data(), impl.ebox_max.data()); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); - } + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("build_edge_face_boxes"); + impl.ebox_min.resize(3 * size_t(n_edges)); + impl.ebox_max.resize(3 * size_t(n_edges)); + if (n_edges > 0) { + build_edge_boxes_kernel<<< + kernel_grid_size(n_edges), KERNEL_BLOCK_SIZE>>>( + impl.vbox_min.data(), impl.vbox_max.data(), + impl.edges.data(), n_edges, impl.ebox_min.data(), + impl.ebox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } - impl.fbox_min.resize(3 * size_t(n_faces)); - impl.fbox_max.resize(3 * size_t(n_faces)); - if (n_faces > 0) { - build_face_boxes_kernel<<< - kernel_grid_size(n_faces), KERNEL_BLOCK_SIZE>>>( - impl.vbox_min.data(), impl.vbox_max.data(), impl.faces.data(), - n_faces, impl.fbox_min.data(), impl.fbox_max.data()); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + impl.fbox_min.resize(3 * size_t(n_faces)); + impl.fbox_max.resize(3 * size_t(n_faces)); + if (n_faces > 0) { + build_face_boxes_kernel<<< + kernel_grid_size(n_faces), KERNEL_BLOCK_SIZE>>>( + impl.vbox_min.data(), impl.vbox_max.data(), + impl.faces.data(), n_faces, impl.fbox_min.data(), + impl.fbox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } } // The CPU normalizes all three BVHs by the vertex box domain. @@ -554,7 +583,10 @@ namespace { // The one synchronization of the build: surfaces any kernel fault and // lets the roots be read back -- all three at once. - IPC_TOOLKIT_CUDA_CHECK(cudaDeviceSynchronize()); + { + IPC_TOOLKIT_PROFILE_BLOCK("synchronize_and_read_roots"); + IPC_TOOLKIT_CUDA_CHECK(cudaDeviceSynchronize()); + } // A hierarchy build that never reaches its root leaves -1 behind. The // traversal always starts at node 0, which would then be an arbitrary @@ -875,6 +907,8 @@ namespace { /// @return The number of candidate pairs emitted. template size_t run_traversal(LBVH::Impl& impl) { + IPC_TOOLKIT_PROFILE_BLOCK("traverse"); + using T = Traversal; const LBVH::Impl::DeviceBVH& source = T::source(impl); const LBVH::Impl::DeviceBVH& target = T::target(impl); @@ -961,20 +995,33 @@ namespace { if (count == 0) { return; } - std::vector h_a(count), h_b(count); - buf.a.download(h_a.data()); - buf.b.download(h_b.data()); - - out.reserve(count); - if (filter.accepts_all()) { - for (size_t k = 0; k < count; ++k) { - out.emplace_back(h_a[k], h_b[k]); - } - } else { - for (size_t k = 0; k < count; ++k) { - if (can_collide(h_a[k], h_b[k])) { + std::vector h_a, h_b; + { + // Staging for the pairs: two int32 arrays the size of the + // candidate set (42 MB on Cloth-Ball), allocated and zeroed. + IPC_TOOLKIT_PROFILE_BLOCK("allocate_host_buffers"); + h_a.resize(count); + h_b.resize(count); + } + { + IPC_TOOLKIT_PROFILE_BLOCK("download_pairs"); + buf.a.download(h_a.data()); + buf.b.download(h_b.data()); + } + + { + IPC_TOOLKIT_PROFILE_BLOCK("construct_candidates"); + out.reserve(count); + if (filter.accepts_all()) { + for (size_t k = 0; k < count; ++k) { out.emplace_back(h_a[k], h_b[k]); } + } else { + for (size_t k = 0; k < count; ++k) { + if (can_collide(h_a[k], h_b[k])) { + out.emplace_back(h_a[k], h_b[k]); + } + } } } } @@ -996,6 +1043,7 @@ namespace { // materialize too, so another call cannot overwrite the buffer while // it is being copied to the host. const std::lock_guard lock(impl.mutex); + IPC_TOOLKIT_PROFILE_BLOCK("cuda::LBVH::detect_candidates"); run_traversal(impl); materialize( Traversal::buffer(impl), filter, can_collide, out); @@ -1080,6 +1128,8 @@ void LBVH::build( Eigen::ConstRef faces, const double inflation_radius) { + IPC_TOOLKIT_PROFILE_BLOCK("cuda::LBVH::build"); + assert(vertices_t0.rows() == vertices_t1.rows()); assert(vertices_t0.cols() == vertices_t1.cols()); assert(vertices_t0.rows() <= std::numeric_limits::max()); @@ -1108,13 +1158,16 @@ void LBVH::build( same_vertices ? device.vertices_t0.data() : device.vertices_t1.data(); // Build vertex boxes on the device (always 3-wide storage). - device.vbox_min.resize(3 * size_t(n_vertices)); - device.vbox_max.resize(3 * size_t(n_vertices)); - build_vertex_boxes_kernel<<< - kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( - device.vertices_t0.data(), d_vertices_t1, n_vertices, dim, - inflation_radius, device.vbox_min.data(), device.vbox_max.data()); - IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + { + IPC_TOOLKIT_PROFILE_BLOCK_CUDA("build_vertex_boxes"); + device.vbox_min.resize(3 * size_t(n_vertices)); + device.vbox_max.resize(3 * size_t(n_vertices)); + build_vertex_boxes_kernel<<< + kernel_grid_size(n_vertices), KERNEL_BLOCK_SIZE>>>( + device.vertices_t0.data(), d_vertices_t1, n_vertices, dim, + inflation_radius, device.vbox_min.data(), device.vbox_max.data()); + IPC_TOOLKIT_CUDA_CHECK(cudaGetLastError()); + } build_from_vertex_boxes(device, dim, n_vertices, edges, faces); } diff --git a/src/ipc/broad_phase/lbvh.cpp b/src/ipc/broad_phase/lbvh.cpp index e39fee3a4..3182be40d 100644 --- a/src/ipc/broad_phase/lbvh.cpp +++ b/src/ipc/broad_phase/lbvh.cpp @@ -401,10 +401,16 @@ void LBVH::detect_candidates( tbb::enumerable_thread_specific> storage; - independent_traversal( - source, target, rightmost_leaves, can_collide, storage); + { + IPC_TOOLKIT_PROFILE_BLOCK("traverse"); + independent_traversal( + source, target, rightmost_leaves, can_collide, storage); + } - merge_thread_local_vectors(storage, candidates); + { + IPC_TOOLKIT_PROFILE_BLOCK("merge_thread_local_candidates"); + merge_thread_local_vectors(storage, candidates); + } } void LBVH::detect_vertex_vertex_candidates( diff --git a/src/ipc/utils/cuda/cuda_profiler.cuh b/src/ipc/utils/cuda/cuda_profiler.cuh new file mode 100644 index 000000000..474751d57 --- /dev/null +++ b/src/ipc/utils/cuda/cuda_profiler.cuh @@ -0,0 +1,70 @@ +#pragma once + +#include + +#ifdef IPC_TOOLKIT_WITH_PROFILER + +#include + +#include + +namespace ipc::cuda { + +/// @brief A ProfilePoint timer that measures device time with CUDA events. +/// +/// Kernel launches are asynchronous, so a host clock around one measures the +/// launch and not the work. Recording events on the stream instead measures +/// when the device ran the enclosed stage. +/// +/// Reading the elapsed time waits on the stop event, so a profiling build +/// serializes stages that would otherwise overlap and its totals are larger +/// than the same build without the profiler. Stage *shares* are what this is +/// for; take absolute build and detect times from an uninstrumented build. +class CudaEventTimer { +public: + CudaEventTimer() + { + IPC_TOOLKIT_CUDA_CHECK(cudaEventCreate(&m_start)); + IPC_TOOLKIT_CUDA_CHECK(cudaEventCreate(&m_stop)); + } + + ~CudaEventTimer() + { + cudaEventDestroy(m_start); + cudaEventDestroy(m_stop); + } + + CudaEventTimer(const CudaEventTimer&) = delete; + CudaEventTimer& operator=(const CudaEventTimer&) = delete; + + void start() { IPC_TOOLKIT_CUDA_CHECK(cudaEventRecord(m_start)); } + + void stop() { IPC_TOOLKIT_CUDA_CHECK(cudaEventRecord(m_stop)); } + + // NOLINTNEXTLINE(readability-identifier-naming) + double getElapsedTimeInMilliSec() const + { + IPC_TOOLKIT_CUDA_CHECK(cudaEventSynchronize(m_stop)); + float ms = 0; + IPC_TOOLKIT_CUDA_CHECK(cudaEventElapsedTime(&ms, m_start, m_stop)); + return static_cast(ms); + } + +private: + cudaEvent_t m_start {}, m_stop {}; +}; + +} // namespace ipc::cuda + +/// @brief Profile a device-side stage, timed with CUDA events. +#define IPC_TOOLKIT_PROFILE_BLOCK_CUDA(...) \ + ipc::ProfilePoint \ + IPC_TOOLKIT_PROFILE_BLOCK_CONCAT( \ + __ipc_cuda_profile_point_, __COUNTER__)(__VA_ARGS__); \ + ZoneScopedN(__VA_ARGS__) + +#else + +#define IPC_TOOLKIT_PROFILE_BLOCK_CUDA(...) ZoneScopedN(__VA_ARGS__) + +#endif diff --git a/tests/src/tests/broad_phase/CMakeLists.txt b/tests/src/tests/broad_phase/CMakeLists.txt index 572cf94ac..5f44702e8 100644 --- a/tests/src/tests/broad_phase/CMakeLists.txt +++ b/tests/src/tests/broad_phase/CMakeLists.txt @@ -19,6 +19,7 @@ set(SOURCES if(IPC_TOOLKIT_WITH_CUDA) list(APPEND SOURCES test_gpu_lbvh.cu + benchmark_lbvh_stages.cu ) endif() diff --git a/tests/src/tests/broad_phase/benchmark_lbvh_stages.cu b/tests/src/tests/broad_phase/benchmark_lbvh_stages.cu new file mode 100644 index 000000000..5c16087da --- /dev/null +++ b/tests/src/tests/broad_phase/benchmark_lbvh_stages.cu @@ -0,0 +1,147 @@ +// Per-stage breakdown of ipc::LBVH and ipc::cuda::LBVH, for the stacked-bar +// benchmark figure. Both broad phases are instrumented with +// IPC_TOOLKIT_PROFILE_BLOCK, so this only has to drive them and hand the +// profiler tree out; the stage names come from the library, not from here. +// +// Device stages are timed with CUDA events (see +// ipc/utils/cuda/cuda_profiler.cuh), which waits on each stage's stop event. +// That serializes stages the uninstrumented build overlaps, so the totals +// here run larger than the fused build and detect times. Use this for stage +// *shares*; take absolute times from a build without the profiler. +// +// Run with: +// IPC_TOOLKIT_BENCH_OUTPUT=stages.json ./ipc_toolkit_tests "[lbvh_stages]" +// +// Environment: +// IPC_TOOLKIT_BENCH_SAMPLES timed calls per scene and phase (default 10) +// IPC_TOOLKIT_BENCH_OUTPUT write the results as JSON to this path + +#include + +#if defined(IPC_TOOLKIT_WITH_CUDA) && defined(IPC_TOOLKIT_WITH_PROFILER) + +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include +#include +#include + +using namespace ipc; + +namespace { + +struct Scene { + std::string name, mesh_t0, mesh_t1; +}; + +const std::vector SCENES = { + { "Cloth-Funnel", "cloth-funnel/227.ply", "cloth-funnel/228.ply" }, + { "Armadillo-Rollers", "armadillo-rollers/326.ply", + "armadillo-rollers/327.ply" }, + { "Rod-Twist", "rod-twist/3036.ply", "rod-twist/3037.ply" }, + { "Cloth-Ball", "cloth_ball92.ply", "cloth_ball93.ply" }, + { "N-Body-Simulation", "n-body-simulation/balls16_18.ply", + "n-body-simulation/balls16_19.ply" }, + { "Puffer-Ball", "puffer-ball/20.ply", "puffer-ball/21.ply" }, +}; + +int env_int(const char* name, const int fallback) +{ + const char* v = std::getenv(name); + return (v != nullptr && *v != '\0') ? std::atoi(v) : fallback; +} + +/// @brief Run `f` `samples` times with the profiler cleared first, and return +/// the resulting scope tree. Each scope carries its own accumulated time_ms +/// and count, so the consumer divides rather than this dividing for it. +template nlohmann::json profile_calls(const int samples, F&& f) +{ + f(); // warm up: first-touch allocation, learned buffer capacities + profiler().clear(); + for (int i = 0; i < samples; ++i) { + f(); + } + return profiler().data(); +} + +} // namespace + +TEST_CASE("Benchmark LBVH stages", "[!benchmark][broad_phase][lbvh_stages]") +{ + tests::skip_if_no_cuda_device(); + + constexpr double inflation_radius = 0; + const int samples = env_int("IPC_TOOLKIT_BENCH_SAMPLES", 10); + + nlohmann::json report; + report["samples"] = samples; + report["scenes"] = nlohmann::json::array(); + + for (const auto& [name, mesh_t0, mesh_t1] : SCENES) { + Eigen::MatrixXd vertices_t0, vertices_t1; + Eigen::MatrixXi edges, faces; + REQUIRE(tests::load_mesh(mesh_t0, vertices_t0, edges, faces)); + REQUIRE(tests::load_mesh(mesh_t1, vertices_t1, edges, faces)); + + nlohmann::json scene; + scene["scene"] = name; + scene["num_faces"] = faces.rows(); + scene["num_edges"] = edges.rows(); + scene["num_vertices"] = vertices_t0.rows(); + + ipc::LBVH cpu; + ipc::cuda::LBVH gpu; + + scene["cpu_build"] = profile_calls(samples, [&]() { + cpu.build(vertices_t0, vertices_t1, edges, faces, inflation_radius); + }); + scene["gpu_build"] = profile_calls(samples, [&]() { + gpu.build(vertices_t0, vertices_t1, edges, faces, inflation_radius); + }); + + // Both trees are now built; time detection against them. The output + // vector is fresh per call, as in the Catch2 detect benchmarks, so + // the candidate storage is allocated here and not reused across + // calls -- on Cloth-Ball that is 42 MB per call. + size_t n_cpu = 0, n_gpu = 0; + scene["cpu_detect"] = profile_calls(samples, [&]() { + std::vector candidates; + cpu.detect_edge_edge_candidates(candidates); + n_cpu = candidates.size(); + }); + scene["gpu_detect"] = profile_calls(samples, [&]() { + std::vector candidates; + gpu.detect_edge_edge_candidates(candidates); + n_gpu = candidates.size(); + }); + // The two broad phases must agree, or the stage split is comparing + // different amounts of work. + REQUIRE(n_gpu == n_cpu); + scene["num_ee_candidates"] = n_cpu; + + report["scenes"].push_back(scene); + } + + profiler().clear(); + + const char* out = std::getenv("IPC_TOOLKIT_BENCH_OUTPUT"); + if (out != nullptr && *out != '\0') { + std::ofstream file(out); + REQUIRE(file.is_open()); + file << report.dump(2) << std::endl; + } else { + WARN("IPC_TOOLKIT_BENCH_OUTPUT unset; stage results not written"); + } +} + +#endif