Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <charm++.h>
       7             : #include <cstddef>
       8             : #include <functional>
       9             : #include <optional>
      10             : #include <unordered_map>
      11             : #include <unordered_set>
      12             : #include <vector>
      13             : 
      14             : #include "DataStructures/TaggedTuple.hpp"
      15             : #include "Domain/Block.hpp"
      16             : #include "Domain/Creators/DomainCreator.hpp"
      17             : #include "Domain/Creators/OptionTags.hpp"
      18             : #include "Domain/Creators/Tags/Domain.hpp"
      19             : #include "Domain/Creators/Tags/InitialExtents.hpp"
      20             : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
      21             : #include "Domain/Domain.hpp"
      22             : #include "Domain/ElementDistribution.hpp"
      23             : #include "Domain/Structure/ElementId.hpp"
      24             : #include "Domain/Structure/InitialElementIds.hpp"
      25             : #include "Domain/Tags.hpp"
      26             : #include "Domain/Tags/ElementDistribution.hpp"
      27             : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
      28             : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      29             : #include "NumericalAlgorithms/Spectral/Quadrature.hpp"
      30             : #include "Parallel/Algorithms/AlgorithmArray.hpp"
      31             : #include "Parallel/DomainDiagnosticInfo.hpp"
      32             : #include "Parallel/GlobalCache.hpp"
      33             : #include "Parallel/Info.hpp"
      34             : #include "Parallel/Local.hpp"
      35             : #include "Parallel/ParallelComponentHelpers.hpp"
      36             : #include "Parallel/Phase.hpp"
      37             : #include "Parallel/Printf/Printf.hpp"
      38             : #include "Parallel/Protocols/ArrayElementsAllocator.hpp"
      39             : #include "Utilities/Literals.hpp"
      40             : #include "Utilities/Numeric.hpp"
      41             : #include "Utilities/ProtocolHelpers.hpp"
      42             : #include "Utilities/System/ParallelInfo.hpp"
      43             : #include "Utilities/TMPL.hpp"
      44             : 
      45             : namespace elliptic {
      46             : 
      47             : /*!
      48             :  * \brief A `Parallel::protocols::ArrayElementsAllocator` that creates array
      49             :  * elements to cover the initial computational domain
      50             :  *
      51             :  * An element is created for every element ID in every block, determined by the
      52             :  * `initial_element_ids` function and the option-created `domain::Tags::Domain`
      53             :  * and `domain::Tags::InitialRefinementLevels`. The elements are distributed
      54             :  * on processors using the `domain::BlockZCurveProcDistribution`. In both cases,
      55             :  * an unordered set of `size_t`s can be passed to the `allocate_array` function
      56             :  * which represents physical processors to avoid placing elements on. `Element`s
      57             :  * are distributed to processors according to their computational costs
      58             :  * determined by the number of grid points.
      59             :  */
      60             : template <size_t Dim>
      61           1 : struct DefaultElementsAllocator
      62             :     : tt::ConformsTo<Parallel::protocols::ArrayElementsAllocator> {
      63             :   template <typename ParallelComponent>
      64           0 :   using array_allocation_tags = tmpl::list<>;
      65             : 
      66             :   template <typename ParallelComponent, typename Metavariables,
      67             :             typename... InitializationTags>
      68           0 :   static void apply(
      69             :       Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
      70             :       const tuples::TaggedTuple<InitializationTags...>& initialization_items,
      71             :       const tuples::tagged_tuple_from_typelist<
      72             :           typename ParallelComponent::array_allocation_tags>&
      73             :       /*array_allocation_items*/
      74             :       = {},
      75             :       const std::unordered_set<size_t>& procs_to_ignore = {}) {
      76             :     auto& local_cache = *Parallel::local_branch(global_cache);
      77             :     auto& element_array =
      78             :         Parallel::get_parallel_component<ParallelComponent>(local_cache);
      79             :     const auto& initial_extents =
      80             :         get<domain::Tags::InitialExtents<Dim>>(initialization_items);
      81             :     const auto basis = Spectral::Basis::Legendre;
      82             :     const auto& quadrature =
      83             :         Parallel::get<elliptic::dg::Tags::Quadrature>(local_cache);
      84             : 
      85             :     const auto& domain = Parallel::get<domain::Tags::Domain<Dim>>(local_cache);
      86             :     const auto& initial_refinement_levels =
      87             :         get<domain::Tags::InitialRefinementLevels<Dim>>(initialization_items);
      88             : 
      89             :     const size_t number_of_procs =
      90             :         Parallel::number_of_procs<size_t>(local_cache);
      91             :     const size_t number_of_nodes =
      92             :         Parallel::number_of_nodes<size_t>(local_cache);
      93             :     const size_t num_of_procs_to_use = number_of_procs - procs_to_ignore.size();
      94             : 
      95             :     const auto& blocks = domain.blocks();
      96             : 
      97             :     const std::optional<domain::ElementWeight>& element_weight =
      98             :         get<domain::Tags::ElementDistribution>(local_cache);
      99             : 
     100             :     domain::BlockZCurveProcDistribution<Dim> element_distribution{};
     101             :     if (element_weight.has_value()) {
     102             :       const std::unordered_map<ElementId<Dim>, double> element_costs =
     103             :           domain::get_element_costs(blocks, initial_refinement_levels,
     104             :                                     initial_extents, element_weight.value(),
     105             :                                     basis, quadrature);
     106             :       element_distribution = domain::BlockZCurveProcDistribution<Dim>{
     107             :           element_costs,   num_of_procs_to_use,
     108             :           blocks,          initial_refinement_levels,
     109             :           initial_extents, procs_to_ignore};
     110             :     }
     111             : 
     112             :     // Will be used to print domain diagnostic info
     113             :     std::vector<size_t> elements_per_core(number_of_procs, 0_st);
     114             :     std::vector<size_t> elements_per_node(number_of_nodes, 0_st);
     115             :     std::vector<size_t> grid_points_per_core(number_of_procs, 0_st);
     116             :     std::vector<size_t> grid_points_per_node(number_of_nodes, 0_st);
     117             : 
     118             :     size_t which_proc = 0;
     119             :     for (const auto& block : blocks) {
     120             :       const size_t grid_points_per_element = alg::accumulate(
     121             :           initial_extents[block.id()], 1_st, std::multiplies<size_t>());
     122             : 
     123             :       const std::vector<ElementId<Dim>> element_ids = initial_element_ids(
     124             :           block.id(), initial_refinement_levels[block.id()]);
     125             : 
     126             :       // Distributed with weighted space filling curve
     127             :       if (element_weight.has_value()) {
     128             :         for (const auto& element_id : element_ids) {
     129             :           const size_t target_proc =
     130             :               element_distribution.get_proc_for_element(element_id);
     131             :           element_array(element_id)
     132             :               .insert(global_cache, initialization_items, target_proc);
     133             : 
     134             :           const size_t target_node =
     135             :               Parallel::node_of<size_t>(target_proc, local_cache);
     136             :           ++elements_per_core[target_proc];
     137             :           ++elements_per_node[target_node];
     138             :           grid_points_per_core[target_proc] += grid_points_per_element;
     139             :           grid_points_per_node[target_node] += grid_points_per_element;
     140             :         }
     141             :       } else {
     142             :         // Distributed with round-robin
     143             :         for (size_t i = 0; i < element_ids.size(); ++i) {
     144             :           while (procs_to_ignore.find(which_proc) != procs_to_ignore.end()) {
     145             :             which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
     146             :           }
     147             : 
     148             :           element_array(ElementId<Dim>(element_ids[i]))
     149             :               .insert(global_cache, initialization_items, which_proc);
     150             : 
     151             :           const size_t target_node =
     152             :               Parallel::node_of<size_t>(which_proc, local_cache);
     153             :           ++elements_per_core[which_proc];
     154             :           ++elements_per_node[target_node];
     155             :           grid_points_per_core[which_proc] += grid_points_per_element;
     156             :           grid_points_per_node[target_node] += grid_points_per_element;
     157             : 
     158             :           which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
     159             :         }
     160             :       }
     161             :     }
     162             :     element_array.doneInserting();
     163             : 
     164             :     Parallel::printf("\n%s\n", domain::diagnostic_info(
     165             :                                    domain.blocks().size(), local_cache,
     166             :                                    elements_per_core, elements_per_node,
     167             :                                    grid_points_per_core, grid_points_per_node));
     168             :   }
     169             : };
     170             : 
     171             : /*!
     172             :  * \brief The parallel component responsible for managing the DG elements that
     173             :  * compose the computational domain
     174             :  *
     175             :  * This parallel component will perform the actions specified by the
     176             :  * `PhaseDepActionList`.
     177             :  *
     178             :  * \note This parallel component is nearly identical to
     179             :  * `Evolution/DiscontinuousGalerkin/DgElementArray.hpp` right now, but will
     180             :  * likely diverge in the future, for instance to support a multigrid domain.
     181             :  *
     182             :  */
     183             : template <typename Metavariables, typename PhaseDepActionList,
     184             :           typename ElementsAllocator =
     185             :               DefaultElementsAllocator<Metavariables::volume_dim>>
     186           1 : struct DgElementArray {
     187           0 :   static constexpr size_t volume_dim = Metavariables::volume_dim;
     188             :   static_assert(
     189             :       tt::assert_conforms_to_v<ElementsAllocator,
     190             :                                Parallel::protocols::ArrayElementsAllocator>);
     191             : 
     192           0 :   using chare_type = Parallel::Algorithms::Array;
     193           0 :   static constexpr bool checkpoint_data = true;
     194           0 :   using metavariables = Metavariables;
     195           0 :   using phase_dependent_action_list = PhaseDepActionList;
     196           0 :   using array_index = ElementId<volume_dim>;
     197             : 
     198           0 :   using const_global_cache_tags = tmpl::list<domain::Tags::Domain<volume_dim>,
     199             :                                              domain::Tags::ElementDistribution>;
     200             : 
     201           0 :   using array_allocation_tags =
     202             :       typename ElementsAllocator::template array_allocation_tags<
     203             :           DgElementArray>;
     204             : 
     205           0 :   using simple_tags_from_options =
     206             :       tmpl::append<Parallel::get_simple_tags_from_options<
     207             :                        Parallel::get_initialization_actions_list<
     208             :                            phase_dependent_action_list>>,
     209             :                    array_allocation_tags>;
     210             : 
     211           0 :   static void allocate_array(
     212             :       Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
     213             :       const tuples::tagged_tuple_from_typelist<simple_tags_from_options>&
     214             :           initialization_items,
     215             :       const tuples::tagged_tuple_from_typelist<array_allocation_tags>&
     216             :           array_allocation_items = {},
     217             :       const std::unordered_set<size_t>& procs_to_ignore = {}) {
     218             :     ElementsAllocator::template apply<DgElementArray>(
     219             :         global_cache, initialization_items, array_allocation_items,
     220             :         procs_to_ignore);
     221             :   }
     222             : 
     223           0 :   static void execute_next_phase(
     224             :       const Parallel::Phase next_phase,
     225             :       Parallel::CProxy_GlobalCache<Metavariables>& global_cache) {
     226             :     auto& local_cache = *Parallel::local_branch(global_cache);
     227             :     Parallel::get_parallel_component<DgElementArray>(local_cache)
     228             :         .start_phase(next_phase);
     229             :   }
     230             : };
     231             : 
     232             : }  // namespace elliptic