cpu/builders/material_deduplication.cpp

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 // Project include(s)
0010 #include "detray/builders/detail/material_deduplication.hpp"
0011
0012 #include "detray/builders/cuboid_portal_generator.hpp"
0013 #include "detray/builders/detector_builder.hpp"
0014 #include "detray/builders/homogeneous_material_builder.hpp"
0015 #include "detray/builders/homogeneous_material_factory.hpp"
0016 #include "detray/builders/material_map_builder.hpp"
0017 #include "detray/builders/material_map_factory.hpp"
0018 #include "detray/builders/surface_factory.hpp"
0019 #include "detray/core/detector.hpp"
0020 #include "detray/definitions/indexing.hpp"
0021 #include "detray/material/detail/material_accessor.hpp"
0022 #include "detray/material/predefined_materials.hpp"
0023 #include "detray/utils/consistency_checker.hpp"
0024 #include "detray/utils/ranges.hpp"
0025
0026 // Detray test include(s)
0027 #include "detray/test/framework/types.hpp"
0028
0029 // Vecmem include(s)
0030 #include <vecmem/memory/host_memory_resource.hpp>
0031
0032 // GTest include(s)
0033 #include <gtest/gtest.h>
0034
0035 // System include(s)
0036 #include <array>
0037 #include <memory>
0038 #include <set>
0039 #include <vector>
0040
0041 using namespace detray;
0042
0043 namespace {
0044
0045 using scalar = test::scalar;
0046 using point2 = test::point2;
0047 using point3 = test::point3;
0048
0049 using metadata_t = test::default_metadata;
0050 using detector_builder_t = detector_builder<metadata_t>;
0051 using detector_t = typename detector_builder_t::detector_type;
0052 using transform3 = dtransform3D<typename detector_t::algebra_type>;
0053 using mat_id = typename detector_t::material::id;
0054 using bin_index_t = axis::multi_bin<2u>;
0055
0056 using rectangle_factory = surface_factory<detector_t, rectangle2D>;
0057 using hom_mat_factory_t = homogeneous_material_factory<detector_t>;
0058 using map_factory_t = material_map_factory<detector_t, bin_index_t>;
0059
0060 constexpr dindex n_volumes{3u};
0061
0062 /// Volume local surface indices
0063 /// @{
0064 // Homogeneous material: identical in all volumes
0065 constexpr dindex shared_slab_sf{0u};
0066 // Homogeneous material: different in every volume
0067 constexpr dindex unique_slab_sf{1u};
0068 // Material map: identical in all volumes
0069 constexpr dindex shared_map_sf{2u};
0070 // Material map: identical in all volumes, but on a smaller surface. The grid
0071 // extent is given by the axis spans, so this surface shares the map as well
0072 constexpr dindex shared_map_small_sf{3u};
0073 // Material map: different in every volume
0074 constexpr dindex unique_map_sf{4u};
0075 /// @}
0076
0077 /// Half length of the material maps in x and y
0078 constexpr scalar map_half_length{10.f * unit<scalar>::mm};
0079
0080 /// Add material map data to the factory @param mat_factory for the surface
0081 /// with volume local index @param sf_index
0082 void add_map_data(map_factory_t &mat_factory, std::size_t sf_index, scalar t,
0083                   const material<scalar> &mat) {
0084   typename map_factory_t::data_type mat_data{sf_index};
0085   std::vector<bin_index_t> bins{};
0086   std::vector<std::size_t> n_bins{5u, 10u};
0087   std::vector<std::vector<scalar>> axis_spans = {
0088       {-map_half_length, map_half_length}, {-map_half_length, map_half_length}};
0089
0090   // Add material for every bin
0091   for (auto [i, j] : detray::views::cartesian_product{
0092            detray::views::iota{0u, 5u}, detray::views::iota{0u, 10u}}) {
0093     bins.push_back({i, j});
0094     mat_data.append(t, mat);
0095     t += 0.25f * unit<scalar>::mm;
0096   }
0097
0098   mat_factory.add_material(mat_id::e_rectangle2D_map, std::move(mat_data),
0099                            std::move(n_bins), std::move(axis_spans),
0100                            std::move(bins));
0101 }
0102
0103 /// Build a detector with @c n_volumes cuboid volumes, in which some surfaces
0104 /// carry identical material
0105 detector_t build_detector(vecmem::memory_resource &mr,
0106                           const volume_builder_options &builder_opts) {
0107   detector_builder_t det_builder{};
0108   det_builder.set_name("material_deduplication");
0109
0110   for (dindex v = 0u; v < n_volumes; ++v) {
0111     auto vbuilder = det_builder.new_volume(volume_id::e_cuboid);
0112     vbuilder->add_volume_placement(
0113         point3{0.f, 0.f, static_cast<scalar>(v) * 100.f * unit<scalar>::mm});
0114     const auto vol_idx{vbuilder->vol_index()};
0115
0116     det_builder.template decorate<homogeneous_material_builder<detector_t>>(
0117         vol_idx);
0118     auto mat_builder =
0119         det_builder.template decorate<material_map_builder<detector_t>>(
0120             vol_idx);
0121
0122     const auto vol_link{
0123         static_cast<typename detector_t::surface_type::navigation_link>(
0124             vol_idx)};
0125     const scalar vs{static_cast<scalar>(v)};
0126
0127     // Surfaces with homogeneous material
0128     auto hom_factory = std::make_shared<hom_mat_factory_t>(
0129         std::make_unique<rectangle_factory>());
0130
0131     hom_factory->push_back({surface_id::e_sensitive,
0132                             transform3(point3{0.f, 0.f, -20.f}), vol_link,
0133                             std::vector<scalar>{10.f, 10.f}});
0134     hom_factory->add_material(mat_id::e_material_slab,
0135                               {1.f * unit<scalar>::mm, silicon<scalar>()});
0136
0137     hom_factory->push_back({surface_id::e_sensitive,
0138                             transform3(point3{0.f, 0.f, -10.f}), vol_link,
0139                             std::vector<scalar>{10.f, 10.f}});
0140     hom_factory->add_material(
0141         mat_id::e_material_slab,
0142         {(1.f + vs) * unit<scalar>::mm, tungsten<scalar>()});
0143
0144     // Surfaces with material maps
0145     auto map_factory =
0146         std::make_shared<map_factory_t>(std::make_unique<rectangle_factory>());
0147
0148     map_factory->push_back({surface_id::e_sensitive,
0149                             transform3(point3{0.f, 0.f, 0.f}), vol_link,
0150                             std::vector<scalar>{10.f, 10.f}});
0151     add_map_data(*map_factory, shared_map_sf, 1.f * unit<scalar>::mm,
0152                  silicon<scalar>());
0153
0154     map_factory->push_back({surface_id::e_sensitive,
0155                             transform3(point3{0.f, 0.f, 10.f}), vol_link,
0156                             std::vector<scalar>{6.f, 4.f}});
0157     add_map_data(*map_factory, shared_map_small_sf, 1.f * unit<scalar>::mm,
0158                  silicon<scalar>());
0159
0160     map_factory->push_back({surface_id::e_sensitive,
0161                             transform3(point3{0.f, 0.f, 20.f}), vol_link,
0162                             std::vector<scalar>{10.f, 10.f}});
0163     add_map_data(*map_factory, unique_map_sf, (1.f + vs) * unit<scalar>::mm,
0164                  gold<scalar>());
0165
0166     // Add a portal box around the volume
0167     auto portal_generator =
0168         std::make_shared<cuboid_portal_generator<detector_t>>(0.1f *
0169                                                               unit<scalar>::mm);
0170
0171     mat_builder->add_surfaces(hom_factory);
0172     mat_builder->add_surfaces(map_factory);
0173     mat_builder->add_surfaces(portal_generator);
0174   }
0175
0176   return det_builder.build(mr, builder_opts);
0177 }
0178
0179 /// @returns the material slab of the surface material at a local point
0180 struct get_material_slab {
0181   template <typename mat_coll_t, concepts::index index_t>
0182   auto operator()(const mat_coll_t &mat_coll, const index_t idx,
0183                   const point2 &loc_p) const -> material_slab<scalar> {
0184     using material_t = typename mat_coll_t::value_type;
0185
0186     if constexpr (concepts::material_map<material_t>) {
0187       // The test surfaces only carry 2D material maps
0188       if constexpr (material_t::dim == 2) {
0189         return detail::material_accessor::get(
0190             mat_coll, idx, typename material_t::point_type{loc_p[0], loc_p[1]});
0191       } else {
0192         return {};
0193       }
0194     } else if constexpr (std::same_as<material_t, material_slab<scalar>>) {
0195       return detail::material_accessor::get(mat_coll, idx, loc_p);
0196     } else {
0197       return {};
0198     }
0199   }
0200 };
0201
0202 /// @returns the material link of the surface with volume local index
0203 /// @param sf_idx in the volume @param vol_idx
0204 auto material_link(const detector_t &det, dindex vol_idx, dindex sf_idx) {
0205   // The surfaces with material are the sensitive surfaces, which are added
0206   // first to every volume
0207   const auto &vol_desc = det.volumes()[vol_idx];
0208   const dindex first_sf{
0209       vol_desc.template sf_link<surface_id::e_sensitive>()[0]};
0210
0211   const auto &sf_desc = det.surfaces()[first_sf + sf_idx];
0212   EXPECT_EQ(sf_desc.volume(), vol_idx);
0213
0214   return sf_desc.material();
0215 }
0216
0217 }  // anonymous namespace
0218
0219 /// Test the content comparison of material grids
0220 GTEST_TEST(detray_builders, material_deduplication_grid_comparison) {
0221   vecmem::host_memory_resource host_mr;
0222   volume_builder_options builder_opts{};
0223   builder_opts.deduplicate(false);
0224   const detector_t det = build_detector(host_mr, builder_opts);
0225
0226   const auto &maps =
0227       det.material_store().template get<mat_id::e_rectangle2D_map>();
0228   ASSERT_EQ(maps.size(), 3u * n_volumes);
0229
0230   // The map content is identical for the first two map surfaces
0231   EXPECT_TRUE(detail::is_identical_grid(maps[0], maps[0]));
0232   EXPECT_TRUE(detail::is_identical_grid(maps[0], maps[1]));
0233   EXPECT_TRUE(detail::is_identical_grid(maps[0], maps[3]));
0234   EXPECT_FALSE(detail::is_identical_grid(maps[0], maps[2]));
0235   // Unique map in the first and second volume
0236   EXPECT_FALSE(detail::is_identical_grid(maps[2], maps[5]));
0237
0238   // In contrast, the grid operator== compares the storage position of the
0239   // non-owning grids
0240   EXPECT_FALSE(maps[0] == maps[1]);
0241
0242   // Material slabs
0243   const material_slab<scalar> slab{silicon<scalar>(), 1.f * unit<scalar>::mm};
0244   EXPECT_TRUE(detail::is_identical(slab, slab));
0245   EXPECT_FALSE(detail::is_identical(
0246       slab, material_slab<scalar>{silicon<scalar>(), 2.f * unit<scalar>::mm}));
0247   EXPECT_FALSE(detail::is_identical(
0248       slab, material_slab<scalar>{tungsten<scalar>(), 1.f * unit<scalar>::mm}));
0249   EXPECT_EQ(detail::material_hash(slab),
0250             detail::material_hash(material_slab<scalar>{
0251                 silicon<scalar>(), 1.f * unit<scalar>::mm}));
0252 }
0253
0254 /// Build a detector with identical material on surfaces in different volumes
0255 /// with and without material deduplication
0256 GTEST_TEST(detray_builders, material_deduplication) {
0257   vecmem::host_memory_resource host_mr;
0258   volume_builder_options builder_opts{};
0259
0260   builder_opts.deduplicate(false);
0261   const detector_t ref_det = build_detector(host_mr, builder_opts);
0262   builder_opts.deduplicate(true);
0263   const detector_t det = build_detector(host_mr, builder_opts);
0264
0265   EXPECT_TRUE(detail::check_consistency(ref_det));
0266   EXPECT_TRUE(detail::check_consistency(det));
0267
0268   ASSERT_EQ(ref_det.volumes().size(), n_volumes);
0269   ASSERT_EQ(det.volumes().size(), n_volumes);
0270   ASSERT_EQ(ref_det.surfaces().size(), det.surfaces().size());
0271
0272   // Without deduplication: One material entry per surface
0273   const auto &ref_store = ref_det.material_store();
0274   EXPECT_EQ(ref_store.template size<mat_id::e_material_slab>(), 2u * n_volumes);
0275   EXPECT_EQ(ref_store.template size<mat_id::e_rectangle2D_map>(),
0276             3u * n_volumes);
0277
0278   // With deduplication: The shared material is only present once
0279   const auto &store = det.material_store();
0280   EXPECT_EQ(store.template size<mat_id::e_material_slab>(), 1u + n_volumes);
0281   EXPECT_EQ(store.template size<mat_id::e_rectangle2D_map>(), 1u + n_volumes);
0282
0283   // Check the material links
0284   std::set<dindex> ref_slab_links;
0285   std::set<dindex> ref_map_links;
0286   std::set<dindex> unique_slab_links;
0287   std::set<dindex> unique_map_links;
0288   for (dindex v = 0u; v < n_volumes; ++v) {
0289     // Every surface has its own material entry
0290     for (dindex sf_idx : {shared_slab_sf, unique_slab_sf}) {
0291       const auto link = material_link(ref_det, v, sf_idx);
0292       EXPECT_EQ(link.id(), mat_id::e_material_slab);
0293       EXPECT_TRUE(ref_slab_links.insert(link.index()).second);
0294     }
0295     for (dindex sf_idx : {shared_map_sf, shared_map_small_sf, unique_map_sf}) {
0296       const auto link = material_link(ref_det, v, sf_idx);
0297       EXPECT_EQ(link.id(), mat_id::e_rectangle2D_map);
0298       EXPECT_TRUE(ref_map_links.insert(link.index()).second);
0299     }
0300
0301     // The shared material is the first entry that was added
0302     auto link = material_link(det, v, shared_slab_sf);
0303     EXPECT_EQ(link.id(), mat_id::e_material_slab);
0304     EXPECT_EQ(link.index(), 0u);
0305
0306     link = material_link(det, v, unique_slab_sf);
0307     EXPECT_EQ(link.id(), mat_id::e_material_slab);
0308     EXPECT_TRUE(unique_slab_links.insert(link.index()).second);
0309
0310     for (dindex sf_idx : {shared_map_sf, shared_map_small_sf}) {
0311       link = material_link(det, v, sf_idx);
0312       EXPECT_EQ(link.id(), mat_id::e_rectangle2D_map);
0313       EXPECT_EQ(link.index(), 0u);
0314     }
0315
0316     link = material_link(det, v, unique_map_sf);
0317     EXPECT_EQ(link.id(), mat_id::e_rectangle2D_map);
0318     EXPECT_TRUE(unique_map_links.insert(link.index()).second);
0319   }
0320   EXPECT_FALSE(unique_slab_links.contains(0u));
0321   EXPECT_FALSE(unique_map_links.contains(0u));
0322
0323   // The material lookups are unchanged on every surface
0324   const std::array<point2, 5> loc_points{point2{0.f, 0.f}, point2{-9.f, -9.f},
0325                                          point2{3.5f, -2.f}, point2{9.f, 9.f},
0326                                          point2{-5.5f, 7.f}};
0327
0328   std::size_t n_checked{0u};
0329   for (const auto &[i, ref_sf_desc] :
0330        detray::views::enumerate(ref_det.surfaces())) {
0331     const auto &sf_desc = det.surfaces()[static_cast<dindex>(i)];
0332
0333     ASSERT_EQ(ref_sf_desc.has_material(), sf_desc.has_material());
0334     if (!sf_desc.has_material()) {
0335       continue;
0336     }
0337     ASSERT_EQ(ref_sf_desc.material().id(), sf_desc.material().id());
0338
0339     for (const auto &p : loc_points) {
0340       const auto ref_slab = ref_store.template visit<get_material_slab>(
0341           ref_sf_desc.material(), p);
0342       const auto slab =
0343           store.template visit<get_material_slab>(sf_desc.material(), p);
0344
0345       EXPECT_TRUE(detail::is_identical(ref_slab, slab))
0346           << "surface " << i << ", point (" << p[0] << ", " << p[1] << ")";
0347       ++n_checked;
0348     }
0349   }
0350   EXPECT_EQ(n_checked, 5u * loc_points.size() * n_volumes);
0351 }