Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <algorithm>
       7             : #include <charm++.h>
       8             : #include <cstddef>
       9             : #include <optional>
      10             : #include <unordered_map>
      11             : #include <unordered_set>
      12             : #include <vector>
      13             : 
      14             : #include "Domain/Creators/Tags/Domain.hpp"
      15             : #include "Domain/Creators/Tags/InitialExtents.hpp"
      16             : #include "Domain/Creators/Tags/InitialRefinementLevels.hpp"
      17             : #include "Domain/ElementDistribution.hpp"
      18             : #include "Domain/Structure/ElementId.hpp"
      19             : #include "Domain/Tags/ElementDistribution.hpp"
      20             : #include "Elliptic/DiscontinuousGalerkin/Tags.hpp"
      21             : #include "NumericalAlgorithms/Convergence/Tags.hpp"
      22             : #include "NumericalAlgorithms/Spectral/Basis.hpp"
      23             : #include "Parallel/GlobalCache.hpp"
      24             : #include "Parallel/Local.hpp"
      25             : #include "Parallel/Printf/Printf.hpp"
      26             : #include "Parallel/Protocols/ArrayElementsAllocator.hpp"
      27             : #include "Parallel/Section.hpp"
      28             : #include "Parallel/Tags/Section.hpp"
      29             : #include "ParallelAlgorithms/LinearSolver/Multigrid/Hierarchy.hpp"
      30             : #include "ParallelAlgorithms/LinearSolver/Multigrid/Tags.hpp"
      31             : #include "Utilities/ErrorHandling/Error.hpp"
      32             : #include "Utilities/PrettyType.hpp"
      33             : #include "Utilities/ProtocolHelpers.hpp"
      34             : #include "Utilities/Serialization/Serialize.hpp"
      35             : #include "Utilities/System/ParallelInfo.hpp"
      36             : #include "Utilities/TaggedTuple.hpp"
      37             : 
      38             : namespace LinearSolver::multigrid {
      39             : 
      40             : /*!
      41             :  * \brief A `Parallel::protocols::ArrayElementsAllocator` that creates array
      42             :  * elements to cover the initial computational domain multiple times at
      43             :  * different refinement levels, suitable for the
      44             :  * `LinearSolver::multigrid::Multigrid` algorithm
      45             :  *
      46             :  * The `initial_element_ids` function and the option-created
      47             :  * `domain::Tags::Domain` and `domain::Tags::InitialRefinementLevels` determine
      48             :  * the initial set of element IDs in the domain. This is taken as the finest
      49             :  * grid in the multigrid hierarchy. Coarser grids are determined by successively
      50             :  * applying `LinearSolver::multigrid::coarsen`, up to
      51             :  * `LinearSolver::multigrid::Tags::InitialCoarseLevels` coarser grids.
      52             :  *
      53             :  * Array elements are created for all element IDs on all grids, meaning they all
      54             :  * share the same parallel component, action list etc. Elements are connected to
      55             :  * their neighbors _on the same grid_ by the
      56             :  * `domain::create_initial_element` function (this is independent of the
      57             :  * multigrid code, the function is typically called in an initialization
      58             :  * action). Elements are connected to their parent and children _across grids_
      59             :  * by the multigrid-tags set here and in the multigrid initialization actions
      60             :  * (see `LinearSolver::multigrid::parent_id` and
      61             :  * `LinearSolver::multigrid::child_ids`).
      62             :  *
      63             :  * Once the multigrid hierarchy is created, the AMR infrastructure registers
      64             :  * the grids and creates the following sections (see `Parallel::Section`).
      65             :  * - `Parallel::Tags::Section<ElementArray, amr::Tags::GridIndex>`: One section
      66             :  *   per grid. All elements are part of a section with this tag.
      67             :  * - `Parallel::Tags::Section<ElementArray, amr::Tags::IsFinestGrid>`: A single
      68             :  *   section that holds only the elements on the finest grid. Holds
      69             :  *   `std::nullopt` on all other elements.
      70             :  *
      71             :  * The elements are distributed on processors using the
      72             :  * `domain::BlockZCurveProcDistribution` for every grid independently. An
      73             :  * unordered set of `size_t`s can be passed to the `apply` function which
      74             :  * represents physical processors to avoid placing elements on.
      75             :  */
      76             : template <size_t Dim, typename OptionsGroup>
      77           1 : struct ElementsAllocator
      78             :     : tt::ConformsTo<Parallel::protocols::ArrayElementsAllocator> {
      79             :   template <typename ElementArray>
      80           0 :   using array_allocation_tags = tmpl::list<>;
      81             : 
      82             :   template <typename ElementArray, typename Metavariables,
      83             :             typename... InitializationTags>
      84           0 :   static void apply(Parallel::CProxy_GlobalCache<Metavariables>& global_cache,
      85             :                     const tuples::TaggedTuple<InitializationTags...>&
      86             :                         original_initialization_items,
      87             :                     const tuples::tagged_tuple_from_typelist<
      88             :                         typename ElementArray::
      89             :                             array_allocation_tags>& /*array_allocation_items*/
      90             :                     = {},
      91             :                     const std::unordered_set<size_t>& procs_to_ignore = {}) {
      92             :     // Copy the initialization items so we can adjust them on each refinement
      93             :     // level
      94             :     auto initialization_items =
      95             :         deserialize<tuples::TaggedTuple<InitializationTags...>>(
      96             :             serialize<tuples::TaggedTuple<InitializationTags...>>(
      97             :                 original_initialization_items)
      98             :                 .data());
      99             :     auto& local_cache = *Parallel::local_branch(global_cache);
     100             :     auto& element_array =
     101             :         Parallel::get_parallel_component<ElementArray>(local_cache);
     102             :     const auto& domain = get<domain::Tags::Domain<Dim>>(local_cache);
     103             :     const auto& initial_extents =
     104             :         get<domain::Tags::InitialExtents<Dim>>(initialization_items);
     105             :     auto& initial_refinement_levels =
     106             :         get<domain::Tags::InitialRefinementLevels<Dim>>(initialization_items);
     107             :     auto& children_refinement_levels =
     108             :         get<Tags::ChildrenRefinementLevels<Dim>>(initialization_items);
     109             :     auto& parent_refinement_levels =
     110             :         get<Tags::ParentRefinementLevels<Dim>>(initialization_items);
     111             :     const auto basis = Spectral::Basis::Legendre;
     112             :     const auto& quadrature =
     113             :         Parallel::get<elliptic::dg::Tags::Quadrature>(local_cache);
     114             :     const std::optional<domain::ElementWeight>& element_weight =
     115             :         get<domain::Tags::ElementDistribution>(local_cache);
     116             :     std::optional<size_t> max_coarse_levels =
     117             :         get<Tags::InitialCoarseLevels<OptionsGroup>>(local_cache);
     118             :     const size_t number_of_procs =
     119             :         Parallel::number_of_procs<size_t>(local_cache);
     120             :     const size_t num_iterations =
     121             :         get<Convergence::Tags::Iterations<OptionsGroup>>(local_cache);
     122             :     if (UNLIKELY(num_iterations == 0)) {
     123             :       Parallel::printf(
     124             :           "%s is disabled (zero iterations). Only creating a single grid.\n",
     125             :           pretty_type::name<OptionsGroup>());
     126             :       max_coarse_levels = 0;
     127             :     }
     128             :     const auto& blocks = domain.blocks();
     129             : 
     130             :     // Determine number of initial multigrid levels
     131             :     size_t num_levels = 1;
     132             :     for (const auto& ref_levs : initial_refinement_levels) {
     133             :       num_levels = std::max(
     134             :           num_levels, *std::max_element(ref_levs.begin(), ref_levs.end()) + 1);
     135             :     }
     136             :     if (max_coarse_levels.has_value()) {
     137             :       num_levels = std::min(num_levels, max_coarse_levels.value() + 1);
     138             :     }
     139             : 
     140             :     // Create the initial grids
     141             :     for (int multigrid_level = static_cast<int>(num_levels) - 1;
     142             :          multigrid_level >= 0; --multigrid_level) {
     143             :       // Store the current grid as child grid before coarsening it
     144             :       children_refinement_levels = initial_refinement_levels;
     145             :       initial_refinement_levels = parent_refinement_levels;
     146             :       // Construct coarsened (parent) grid
     147             :       if (multigrid_level > 0) {
     148             :         parent_refinement_levels =
     149             :             LinearSolver::multigrid::coarsen(initial_refinement_levels);
     150             :       }
     151             :       // Create element IDs for all elements on this level
     152             :       std::vector<ElementId<Dim>> element_ids{};
     153             :       for (const auto& block : blocks) {
     154             :         const std::vector<ElementId<Dim>> block_element_ids =
     155             :             initial_element_ids(block.id(),
     156             :                                 initial_refinement_levels[block.id()],
     157             :                                 static_cast<size_t>(multigrid_level));
     158             :         element_ids.insert(element_ids.begin(), block_element_ids.begin(),
     159             :                            block_element_ids.end());
     160             :       }
     161             :       // Create the elements for this refinement level and distribute them among
     162             :       // processors
     163             :       const size_t num_of_procs_to_use =
     164             :           static_cast<size_t>(sys::number_of_procs()) - procs_to_ignore.size();
     165             :       // Distributed with weighted space filling curve
     166             :       if (element_weight.has_value()) {
     167             :         const std::unordered_map<ElementId<Dim>, double> element_costs =
     168             :             domain::get_element_costs(blocks, initial_refinement_levels,
     169             :                                       initial_extents, element_weight.value(),
     170             :                                       basis, quadrature);
     171             :         const domain::BlockZCurveProcDistribution<Dim> element_distribution{
     172             :             element_costs,   num_of_procs_to_use,
     173             :             blocks,          initial_refinement_levels,
     174             :             initial_extents, procs_to_ignore};
     175             : 
     176             :         for (const auto& element_id : element_ids) {
     177             :           const size_t target_proc =
     178             :               element_distribution.get_proc_for_element(element_id);
     179             :           element_array(element_id)
     180             :               .insert(global_cache, initialization_items, target_proc);
     181             :         }
     182             :       } else {
     183             :         // Distributed with round-robin
     184             :         size_t which_proc = 0;
     185             :         for (const auto& element_id : element_ids) {
     186             :           while (procs_to_ignore.find(which_proc) != procs_to_ignore.end()) {
     187             :             which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
     188             :           }
     189             : 
     190             :           element_array(element_id)
     191             :               .insert(global_cache, initialization_items, which_proc);
     192             : 
     193             :           which_proc = which_proc + 1 == number_of_procs ? 0 : which_proc + 1;
     194             :         }
     195             :       }
     196             :     }
     197             :     element_array.doneInserting();
     198             :   }
     199             : };
     200             : 
     201             : }  // namespace LinearSolver::multigrid