test/device/propagator_test.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s).
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/units.hpp"
0014 #include "detray/navigation/caching_navigator.hpp"
0015 #include "detray/propagator/actors.hpp"
0016 #include "detray/propagator/actors/parameter_updater.hpp"
0017 #include "detray/propagator/base_actor.hpp"
0018 #include "detray/propagator/propagator.hpp"
0019 #include "detray/propagator/rk_stepper.hpp"
0020 #include "detray/tracks/tracks.hpp"
0021 #include "detray/utils/logging.hpp"
0022
0023 // Detray test include(s)
0024 #include "detray/test/common/track_generators.hpp"
0025 #include "detray/test/framework/types.hpp"
0026 #include "detray/test/utils/inspectors.hpp"
0027 #include "detray/test/utils/step_tracer.hpp"
0028
0029 // Vecmem include(s)
0030 #include <vecmem/memory/memory_resource.hpp>
0031
0032 // Covfie include(s)
0033 #include <covfie/core/field.hpp>
0034
0035 // GTest include(s).
0036 #include <gtest/gtest.h>
0037
0038 // System include(s)
0039 #include <stdexcept>
0040
0041 namespace detray {
0042
0043 // These types are identical in host and device code for all bfield types
0044 using test_algebra = test::algebra;
0045 using scalar = dscalar<test_algebra>;
0046 using vector3 = dvector3D<test_algebra>;
0047 using point3 = dpoint3D<test_algebra>;
0048 using test_track = free_track_parameters<test_algebra>;
0049 using free_matrix_t = free_matrix<test_algebra>;
0050
0051 constexpr std::size_t cache_size{navigation::default_cache_size};
0052
0053 // Navigator
0054 template <typename detector_t, typename inspector_t>
0055 using navigator_w_insp_t =
0056     caching_navigator<detector_t, cache_size, inspector_t>;
0057 template <typename detector_t>
0058 using navigator_t = navigator_w_insp_t<detector_t, navigation::void_inspector>;
0059 template <typename detector_t>
0060 using intersection_t = typename navigator_t<detector_t>::intersection_type;
0061
0062 // Stepper
0063 using constraints_t = constrained_step<scalar>;
0064 template <typename bfield_view_t>
0065 using rk_stepper_t = rk_stepper<bfield_view_t, test_algebra, constraints_t>;
0066
0067 // Track generator
0068 using generator_t = uniform_track_generator<test_track>;
0069
0070 /// Test tolerance
0071 constexpr scalar is_close{1e-4f};
0072
0073 /// Test configuration
0074 struct propagator_test_config {
0075   generator_t::configuration track_generator;
0076   propagation::config propagation;
0077 };
0078
0079 // Assemble actor chain type
0080 using step_tracer_host_t = actor::step_tracer<test_algebra, vecmem::vector>;
0081 using step_tracer_device_t =
0082     actor::step_tracer<test_algebra, vecmem::device_vector>;
0083 using pathlimit_aborter_t = actor::pathlimit_aborter<scalar>;
0084 using parameter_setter_t = actor::parameter_setter<test_algebra>;
0085 using actor_chain_host_t = actor_chain<
0086     pathlimit_aborter_t,
0087     actor::parameter_updater<test_algebra,
0088                              actor::pointwise_material_interactor<test_algebra>,
0089                              step_tracer_host_t>>;
0090 using actor_chain_device_t = actor_chain<
0091     pathlimit_aborter_t,
0092     actor::parameter_updater<test_algebra,
0093                              actor::pointwise_material_interactor<test_algebra>,
0094                              step_tracer_device_t>>;
0095
0096 /// Precompute the tracks
0097 template <typename track_generator_t = uniform_track_generator<test_track>>
0098 inline auto generate_tracks(
0099     vecmem::memory_resource *mr,
0100     const typename track_generator_t::configuration &cfg = {}) {
0101   // Track collection
0102   vecmem::vector<typename track_generator_t::track_type> tracks(mr);
0103
0104   // Iterate through uniformly distributed momentum directions
0105   for (auto track : track_generator_t{cfg}) {
0106     // Put it into vector of trajectories
0107     tracks.push_back(track);
0108   }
0109
0110   return tracks;
0111 }
0112
0113 /// Test function for propagator on the host
0114 template <typename bfield_bknd_t, typename host_detector_t>
0115 inline auto run_propagation_host(vecmem::memory_resource *mr,
0116                                  const host_detector_t &det,
0117                                  const propagation::config &cfg,
0118                                  const covfie::field<bfield_bknd_t> &field,
0119                                  const vecmem::vector<test_track> &tracks)
0120     -> vecmem::jagged_vector<detail::step_data<test_algebra>> {
0121   // Construct propagator from stepper and navigator
0122   using host_stepper_t =
0123       rk_stepper_t<typename covfie::field<bfield_bknd_t>::view_t>;
0124   using host_navigator_t =
0125       navigator_w_insp_t<host_detector_t, navigation::print_inspector>;
0126   using propagator_host_t =
0127       propagator<host_stepper_t, host_navigator_t, actor_chain_host_t>;
0128
0129   propagator_host_t p{cfg};
0130
0131   // Create vector for track recording
0132   vecmem::jagged_vector<detail::step_data<test_algebra>> host_steps(mr);
0133
0134   for (const auto &trk : tracks) {
0135     // Create the propagator state
0136     step_tracer_host_t::state tracer_state{*mr};
0137     typename pathlimit_aborter_t::state pathlimit_state{
0138         cfg.stepping.path_limit};
0139     actor::parameter_updater_state<test_algebra> updater_state{cfg};
0140     actor::pointwise_material_interactor<test_algebra>::state
0141         interactor_state{};
0142
0143     auto actor_states = detray::tie(tracer_state, pathlimit_state,
0144                                     updater_state, interactor_state);
0145
0146     typename covfie::field<bfield_bknd_t>::view_t bfield_view(field);
0147     typename propagator_host_t::state state(trk, bfield_view, det);
0148
0149     state.stepping().template set_constraint<step::constraint::e_accuracy>(
0150         cfg.stepping.step_constraint);
0151
0152     // Run propagation
0153     if (!p.propagate(state, actor_states)) {
0154       DETRAY_FATAL_HOST(state.navigation().inspector().to_string());
0155       throw std::runtime_error("Host propagation failed");
0156     } else if (tracer_state.get_step_data().empty()) {
0157       DETRAY_FATAL_HOST(state.navigation().inspector().to_string());
0158       throw std::runtime_error(
0159           "Host propagation did not record reference data correctly");
0160     }
0161
0162     host_steps.push_back(std::move(tracer_state.get_step_data()));
0163   }
0164
0165   return host_steps;
0166 }
0167
0168 /// Compare the results between host and device propagation
0169 inline void compare_propagation_results(
0170     const vecmem::jagged_vector<detail::step_data<test_algebra>> &host_steps,
0171     const vecmem::jagged_vector<detail::step_data<test_algebra>>
0172         &device_steps) {
0173   // Make sure the same number of tracks were tested on both backends
0174   ASSERT_EQ(host_steps.size(), device_steps.size());
0175
0176   // Compare the positions and pathlengths
0177   for (unsigned int i = 0u; i < host_steps.size(); i++) {
0178     EXPECT_TRUE(host_steps[i].size() > 0u)
0179         << "No surfaces forund on track " << i;
0180
0181     // Recorded as many positions as path lengths
0182     EXPECT_EQ(host_steps[i].size(), device_steps[i].size());
0183
0184     for (unsigned int j = 0u;
0185          j < math::min(host_steps[i].size(), device_steps[i].size()); j++) {
0186       // Compare recorded path lengths along track
0187       const auto &host_step = host_steps[i][j];
0188       const auto &device_step = device_steps[i][j];
0189
0190       EXPECT_NEAR(host_step.path_length, device_step.path_length,
0191                   host_step.path_length * is_close)
0192           << "ERROR: Path length at track " << i << " step " << j << std::endl;
0193
0194       // Compare recorded positions along track
0195       const point3 &host_pos = host_step.track_params.pos();
0196       const point3 &device_pos = device_step.track_params.pos();
0197
0198       auto relative_error = static_cast<point3>(1.f / host_step.path_length *
0199                                                 (host_pos - device_pos));
0200
0201       EXPECT_NEAR(vector::norm(relative_error), 0.f, is_close)
0202           << "ERROR: Position at track " << i << " step " << j << ": ["
0203           << host_pos[0] << ", " << host_pos[1] << ", " << host_pos[2]
0204           << "] (host), [" << device_pos[0] << ", " << device_pos[1] << ", "
0205           << device_pos[2] << "] (device)" << std::endl;
0206
0207       // Compare the Jacobians
0208       const free_matrix_t &host_J = host_step.jacobian;
0209       const free_matrix_t &device_J = device_step.jacobian;
0210
0211       for (std::size_t row = 0u; row < e_free_size; row++) {
0212         for (std::size_t col = 0u; col < e_free_size; col++) {
0213           scalar host_val = getter::element(host_J, row, col);
0214
0215           scalar device_val = getter::element(device_J, row, col);
0216
0217           ASSERT_NEAR((host_val - device_val) / host_step.path_length, 0.f,
0218                       is_close)
0219               << "ERROR: matrix element mismatch at row " << row << ", col "
0220               << col << std::endl;
0221         }
0222       }
0223     }
0224   }
0225 }
0226
0227 }  // namespace detray