From 759b362b2943da1e3d36b62714a51c82595de410 Mon Sep 17 00:00:00 2001 From: Francesco Biscani Date: Mon, 30 Jul 2018 16:58:43 +0200 Subject: [PATCH] More testing work. --- test/CMakeLists.txt | 1 + test/accel_ordering.cpp | 186 ++++++++++++++++++++++++++++++++++++++++ test/accuracy.cpp | 3 +- 3 files changed, 189 insertions(+), 1 deletion(-) create mode 100644 test/accel_ordering.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4cb87d8..2e29726 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -11,3 +11,4 @@ endfunction() ADD_RAKAU_TESTCASE(accuracy) ADD_RAKAU_TESTCASE(update) ADD_RAKAU_TESTCASE(simd) +ADD_RAKAU_TESTCASE(accel_ordering) diff --git a/test/accel_ordering.cpp b/test/accel_ordering.cpp new file mode 100644 index 0000000..6fca7b9 --- /dev/null +++ b/test/accel_ordering.cpp @@ -0,0 +1,186 @@ +// Copyright 2018 Francesco Biscani (bluescarni@gmail.com) +// +// This file is part of the rakau library. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#include + +#define CATCH_CONFIG_MAIN +#include "catch.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include "test_utils.hpp" + +using namespace rakau; +using namespace rakau_test; + +using fp_types = std::tuple; + +static std::mt19937 rng; + +TEST_CASE("accelerations ordering") +{ + tuple_for_each(fp_types{}, [](auto x) { + using fp_type = decltype(x); + constexpr auto bsize = static_cast(10), theta = static_cast(.01); + constexpr auto s = 10000u; + // Get random particles in a box 1/10th the size of the domain. + auto parts = get_uniform_particles<3>(s, bsize / fp_type(10), rng); + octree t(bsize, {parts.begin() + s, parts.begin() + 2u * s, parts.begin() + 3u * s, parts.begin()}, s, + 16, 256); + // Select randomly some particle indices to track. + using size_type = typename decltype(t)::size_type; + std::vector track_idx(100); + std::uniform_int_distribution idist(0, s - 1u); + std::generate(track_idx.begin(), track_idx.end(), [&idist]() { return idist(rng); }); + // Get the exact accelerations on the particles. + std::vector> exact_accs; + for (auto idx : track_idx) { + exact_accs.emplace_back(t.exact_acc_o(idx)); + } + // Now let's rotate the particles around the z axis by a random amount. + std::uniform_real_distribution urd(fp_type(0), fp_type(2) * boost::math::constants::pi()); + auto rot_angle = urd(rng); + t.update_particles_u([rot_angle](const auto &r) { + for (auto i = 0u; i < s; ++i) { + const auto x0 = r[0].first[i]; + const auto y0 = r[1].first[i]; + const auto z0 = r[2].first[i]; + const auto r0 = std::hypot(x0, y0, z0); + const auto th0 = std::acos(z0 / r0); + const auto phi0 = std::atan2(y0, x0); + const auto phi1 = phi0 + rot_angle; + r[0].first[i] = r0 * std::sin(th0) * std::cos(phi1); + r[1].first[i] = r0 * std::sin(th0) * std::sin(phi1); + r[2].first[i] = r0 * std::cos(th0); + } + }); + // Now let's get the tree-calculated accelerations. + std::array, 3> t_accs; + t.accs_o(t_accs, theta); + // Let's compare the absolute values of the accelerations. + for (size_type i = 0; i < track_idx.size(); ++i) { + const auto eacc = std::sqrt(exact_accs[i][0] * exact_accs[i][0] + exact_accs[i][1] * exact_accs[i][1] + + exact_accs[i][2] * exact_accs[i][2]); + const auto tacc = std::sqrt(t_accs[0][track_idx[i]] * t_accs[0][track_idx[i]] + + t_accs[1][track_idx[i]] * t_accs[1][track_idx[i]] + + t_accs[2][track_idx[i]] * t_accs[2][track_idx[i]]); + const auto rdiff = std::abs(eacc - tacc) / eacc; + std::cout << "rdiff=" << rdiff << '\n'; + if (std::numeric_limits::is_iec559) { + // The usual arbitrary values. + if (std::is_same_v) { + REQUIRE(rdiff <= fp_type(2E-11)); + } else { + REQUIRE(rdiff <= fp_type(2E-3)); + } + } + } + // Do another rotation and re-check. + rot_angle = urd(rng); + t.update_particles_u([rot_angle](const auto &r) { + for (auto i = 0u; i < s; ++i) { + const auto x0 = r[0].first[i]; + const auto y0 = r[1].first[i]; + const auto z0 = r[2].first[i]; + const auto r0 = std::hypot(x0, y0, z0); + const auto th0 = std::acos(z0 / r0); + const auto phi0 = std::atan2(y0, x0); + const auto phi1 = phi0 + rot_angle; + r[0].first[i] = r0 * std::sin(th0) * std::cos(phi1); + r[1].first[i] = r0 * std::sin(th0) * std::sin(phi1); + r[2].first[i] = r0 * std::cos(th0); + } + }); + t.accs_o(t_accs, theta); + for (size_type i = 0; i < track_idx.size(); ++i) { + const auto eacc = std::sqrt(exact_accs[i][0] * exact_accs[i][0] + exact_accs[i][1] * exact_accs[i][1] + + exact_accs[i][2] * exact_accs[i][2]); + const auto tacc = std::sqrt(t_accs[0][track_idx[i]] * t_accs[0][track_idx[i]] + + t_accs[1][track_idx[i]] * t_accs[1][track_idx[i]] + + t_accs[2][track_idx[i]] * t_accs[2][track_idx[i]]); + const auto rdiff = std::abs(eacc - tacc) / eacc; + std::cout << "rdiff=" << rdiff << '\n'; + if (std::numeric_limits::is_iec559) { + if (std::is_same_v) { + REQUIRE(rdiff <= fp_type(2E-11)); + } else { + REQUIRE(rdiff <= fp_type(2E-3)); + } + } + } + // Add a small delta to all coordinates and re-check. + std::uniform_real_distribution urd2(fp_type(0), fp_type(3)); + auto delta = urd2(rng); + t.update_particles_u([delta](const auto &r) { + for (auto i = 0u; i < s; ++i) { + r[0].first[i] += delta; + r[1].first[i] += delta; + r[2].first[i] += delta; + } + }); + t.accs_o(t_accs, theta); + for (size_type i = 0; i < track_idx.size(); ++i) { + const auto eacc = std::sqrt(exact_accs[i][0] * exact_accs[i][0] + exact_accs[i][1] * exact_accs[i][1] + + exact_accs[i][2] * exact_accs[i][2]); + const auto tacc = std::sqrt(t_accs[0][track_idx[i]] * t_accs[0][track_idx[i]] + + t_accs[1][track_idx[i]] * t_accs[1][track_idx[i]] + + t_accs[2][track_idx[i]] * t_accs[2][track_idx[i]]); + const auto rdiff = std::abs(eacc - tacc) / eacc; + std::cout << "rdiff=" << rdiff << '\n'; + if (std::numeric_limits::is_iec559) { + if (std::is_same_v) { + REQUIRE(rdiff <= fp_type(2E-11)); + } else { + REQUIRE(rdiff <= fp_type(2E-3)); + } + } + } + // Now let's update the exact accelerations on the tracked particles. + for (size_type i = 0; i < track_idx.size(); ++i) { + exact_accs[i] = t.exact_acc_o(track_idx[i]); + } + // Let's remove the delta. + t.update_particles_u([delta](const auto &r) { + for (auto i = 0u; i < s; ++i) { + r[0].first[i] -= delta; + r[1].first[i] -= delta; + r[2].first[i] -= delta; + } + }); + // Compute the tree accelerations. + t.accs_o(t_accs, theta); + // Verify they are close to the original ones (i.e., close in all components, + // not only absolute value). + for (size_type i = 0; i < track_idx.size(); ++i) { + const auto [ex, ey, ez] = exact_accs[i]; + const auto xrdiff = std::abs((ex - t_accs[0][track_idx[i]]) / ex); + const auto yrdiff = std::abs((ey - t_accs[1][track_idx[i]]) / ey); + const auto zrdiff = std::abs((ez - t_accs[2][track_idx[i]]) / ez); + if (std::numeric_limits::is_iec559) { + if (std::is_same_v) { + REQUIRE(xrdiff <= fp_type(2E-11)); + REQUIRE(yrdiff <= fp_type(2E-11)); + REQUIRE(zrdiff <= fp_type(2E-11)); + } else { + REQUIRE(xrdiff <= fp_type(2E-3)); + REQUIRE(yrdiff <= fp_type(2E-3)); + REQUIRE(zrdiff <= fp_type(2E-3)); + } + } + } + }); +} \ No newline at end of file diff --git a/test/accuracy.cpp b/test/accuracy.cpp index 8a60575..a496fff 100644 --- a/test/accuracy.cpp +++ b/test/accuracy.cpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -96,7 +97,7 @@ TEST_CASE("accuracy") std::cout << "\n\n\ntot_max_x_diff=" << tot_max_x_diff << '\n'; std::cout << "tot_max_y_diff=" << tot_max_y_diff << '\n'; std::cout << "tot_max_z_diff=" << tot_max_z_diff << "\n\n\n"; - if constexpr (std::is_same_v) { + if constexpr (std::is_same_v && std::numeric_limits::is_iec559) { // These numbers are, of course, totally arbitrary, based // on the fact that 'double' is actually double-precision, // and derived experimentally.