Line data    Source code

       1           0 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : #pragma once
       5             : 
       6             : #include <cstddef>
       7             : #include <sstream>
       8             : #include <string>
       9             : 
      10             : #include "DataStructures/DataBox/DataBox.hpp"
      11             : #include "DataStructures/Variables.hpp"
      12             : #include "DataStructures/VariablesTag.hpp"
      13             : #include "Domain/Creators/Tags/Domain.hpp"
      14             : #include "Domain/ElementLogicalCoordinates.hpp"
      15             : #include "Domain/FunctionsOfTime/FunctionOfTime.hpp"
      16             : #include "Domain/TagsTimeDependent.hpp"
      17             : #include "IO/Logging/Verbosity.hpp"
      18             : #include "NumericalAlgorithms/Interpolation/IrregularInterpolant.hpp"
      19             : #include "Parallel/GlobalCache.hpp"
      20             : #include "Parallel/Info.hpp"
      21             : #include "Parallel/Invoke.hpp"
      22             : #include "Parallel/Printf/Printf.hpp"
      23             : #include "ParallelAlgorithms/Interpolation/InterpolationTargetDetail.hpp"
      24             : #include "ParallelAlgorithms/Interpolation/Tags.hpp"
      25             : #include "Utilities/Algorithm.hpp"
      26             : #include "Utilities/Gsl.hpp"
      27             : #include "Utilities/PrettyType.hpp"
      28             : #include "Utilities/StdHelpers.hpp"
      29             : #include "Utilities/System/ParallelInfo.hpp"
      30             : #include "Utilities/TMPL.hpp"
      31             : #include "Utilities/TaggedTuple.hpp"
      32             : 
      33             : /// \cond
      34             : namespace intrp {
      35             : template <typename Metavariables, typename InterpolationTargetTag>
      36             : struct InterpolationTarget;
      37             : namespace Actions {
      38             : template <typename InterpolationTargetTag>
      39             : struct InterpolationTargetReceiveVars;
      40             : }  // namespace Actions
      41             : }  // namespace intrp
      42             : /// \endcond
      43             : 
      44             : namespace intrp {
      45             : 
      46             : namespace interpolator_detail {
      47             : 
      48             : // Interpolates data onto a set of points desired by an InterpolationTarget.
      49             : template <typename InterpolationTargetTag, typename Metavariables,
      50             :           typename DbTags>
      51             : void interpolate_data(
      52             :     const gsl::not_null<db::DataBox<DbTags>*> box,
      53             :     Parallel::GlobalCache<Metavariables>& cache,
      54             :     const typename InterpolationTargetTag::temporal_id::type& temporal_id) {
      55             :   db::mutate_apply<
      56             :       tmpl::list<
      57             :           ::intrp::Tags::InterpolatedVarsHolders<Metavariables>,
      58             :           ::intrp::Tags::VolumeVarsInfo<
      59             :               Metavariables, typename InterpolationTargetTag::temporal_id>>,
      60             :       tmpl::list<domain::Tags::Domain<Metavariables::volume_dim>>>(
      61             :       [&cache, &temporal_id](
      62             :           const gsl::not_null<typename ::intrp::Tags::InterpolatedVarsHolders<
      63             :               Metavariables>::type*>
      64             :               holders,
      65             :           const gsl::not_null<typename ::intrp::Tags::VolumeVarsInfo<
      66             :               Metavariables,
      67             :               typename InterpolationTargetTag::temporal_id>::type*>
      68             :               volume_vars_info,
      69             :           const Domain<Metavariables::volume_dim>& domain) {
      70             :         auto& interp_info =
      71             :             get<Vars::HolderTag<InterpolationTargetTag, Metavariables>>(
      72             :                 *holders)
      73             :                 .infos.at(temporal_id);
      74             : 
      75             :         // Avoid compiler warning for unused variable in some 'if
      76             :         // constexpr' branches.
      77             :         (void)cache;
      78             : 
      79             :         for (auto& volume_info_outer : *volume_vars_info) {
      80             :           // Are we at the right time?
      81             :           if (volume_info_outer.first != temporal_id) {
      82             :             continue;
      83             :           }
      84             : 
      85             :           // Get list of ElementIds that have the correct temporal_id and that
      86             :           // have not yet been interpolated.
      87             :           std::vector<ElementId<Metavariables::volume_dim>> element_ids;
      88             : 
      89             :           for (const auto& volume_info_inner : volume_info_outer.second) {
      90             :             // Have we interpolated this element before?
      91             :             if (interp_info.interpolation_is_done_for_these_elements.find(
      92             :                     volume_info_inner.first) ==
      93             :                 interp_info.interpolation_is_done_for_these_elements.end()) {
      94             :               interp_info.interpolation_is_done_for_these_elements.emplace(
      95             :                   volume_info_inner.first);
      96             :               element_ids.push_back(volume_info_inner.first);
      97             :             }
      98             :           }
      99             : 
     100             :           // Get element logical coordinates.
     101             :           const auto element_coord_holders = element_logical_coordinates(
     102             :               element_ids, interp_info.block_coord_holders);
     103             : 
     104             :           // Construct local vars and interpolate.
     105             :           for (const auto& element_coord_pair : element_coord_holders) {
     106             :             const auto& element_id = element_coord_pair.first;
     107             :             auto& volume_info = volume_info_outer.second.at(element_id);
     108             :             auto& vars_to_interpolate =
     109             :                 get<::intrp::Tags::VarsToInterpolateToTarget<
     110             :                     InterpolationTargetTag>>(volume_info.vars_to_interpolate);
     111             : 
     112             :             if constexpr (InterpolationTarget_detail::
     113             :                               has_compute_vars_to_interpolate_v<
     114             :                                   InterpolationTargetTag>) {
     115             :               if (vars_to_interpolate.size() == 0) {
     116             :                 // vars_to_interpolate has not been filled for
     117             :                 // this element at this temporal_id.  So fill it.
     118             :                 vars_to_interpolate.initialize(
     119             :                     volume_info.source_vars_from_element
     120             :                         .number_of_grid_points());
     121             : 
     122             :                 InterpolationTarget_detail::compute_dest_vars_from_source_vars<
     123             :                     InterpolationTargetTag>(
     124             :                     make_not_null(&vars_to_interpolate),
     125             :                     volume_info.source_vars_from_element, domain,
     126             :                     volume_info.mesh, element_id, cache, temporal_id);
     127             :               }
     128             :             }
     129             : 
     130             :             // Now interpolate.
     131             :             const auto& element_coord_holder = element_coord_pair.second;
     132             :             intrp::Irregular<Metavariables::volume_dim> interpolator(
     133             :                 volume_info.mesh, element_coord_holder.element_logical_coords);
     134             :             // This first branch is used if compute_vars_to_interpolate exists
     135             :             // or if the vars_to_interpolate_to_target is a subset of the
     136             :             // interpolator_source_vars.
     137             :             if constexpr (
     138             :                 InterpolationTarget_detail::has_compute_vars_to_interpolate_v<
     139             :                     InterpolationTargetTag> or
     140             :                 not std::is_same_v<
     141             :                     tmpl::list_difference<
     142             :                         typename Metavariables::interpolator_source_vars,
     143             :                         typename InterpolationTargetTag::
     144             :                             vars_to_interpolate_to_target>,
     145             :                     tmpl::list<>>) {
     146             :               interp_info.vars.emplace_back(
     147             :                   interpolator.interpolate(vars_to_interpolate));
     148             :             } else {
     149             :               // If compute_vars_to_interpolate does not exist and
     150             :               // vars_to_interpolate_to_target isn't a subset of
     151             :               // interpolator_source_vars, then
     152             :               // volume_info.source_vars_from_element is the same as
     153             :               // volume_info.vars_to_interpolate.
     154             :               interp_info.vars.emplace_back(interpolator.interpolate(
     155             :                   volume_info.source_vars_from_element));
     156             :             }
     157             :             interp_info.global_offsets.emplace_back(
     158             :                 element_coord_holder.offsets);
     159             :           }
     160             :         }
     161             :       },
     162             :       box);
     163             : }
     164             : }  // namespace interpolator_detail
     165             : 
     166             : /// Check if we have enough information to interpolate.  If so, do the
     167             : /// interpolation and send data to the InterpolationTarget.
     168             : template <typename InterpolationTargetTag, typename Metavariables,
     169             :           typename DbTags>
     170           1 : void try_to_interpolate(
     171             :     const gsl::not_null<db::DataBox<DbTags>*> box,
     172             :     const gsl::not_null<Parallel::GlobalCache<Metavariables>*> cache,
     173             :     const typename InterpolationTargetTag::temporal_id::type& temporal_id) {
     174             :   const auto& holders =
     175             :       db::get<Tags::InterpolatedVarsHolders<Metavariables>>(*box);
     176             :   const auto& vars_infos =
     177             :       get<Vars::HolderTag<InterpolationTargetTag, Metavariables>>(holders)
     178             :           .infos;
     179             :   std::stringstream ss{};
     180             :   const ::Verbosity& verbosity = Parallel::get<intrp::Tags::Verbosity>(*cache);
     181             :   const bool debug_print = verbosity >= ::Verbosity::Debug;
     182             : 
     183             :   if (debug_print) {
     184             :     ss << InterpolationTarget_detail::interpolator_output_prefix<
     185             :               InterpolationTargetTag>(temporal_id)
     186             :        << ", ";
     187             :   }
     188             : 
     189             :   // If we don't yet have any points for this InterpolationTarget at
     190             :   // this temporal_id, we should exit (we can't interpolate anyway).
     191             :   if (vars_infos.count(temporal_id) == 0) {
     192             :     if (debug_print) {
     193             :       ss << "no points received yet. Can't do anything.\n";
     194             :       Parallel::printf("%s", ss.str());
     195             :     }
     196             :     return;
     197             :   }
     198             : 
     199             :   interpolator_detail::interpolate_data<InterpolationTargetTag, Metavariables>(
     200             :       box, *cache, temporal_id);
     201             : 
     202             :   // Send interpolated data only if interpolation has been done on all
     203             :   // of the local elements.
     204             :   const auto& all_num_elements =
     205             :       db::get<Tags::NumberOfElements<Metavariables::volume_dim>>(*box);
     206             :   const std::string& target_name = pretty_type::name<InterpolationTargetTag>();
     207             :   if (not all_num_elements.contains(target_name)) {
     208             :     ERROR("Target name " << target_name << " not in NumberOfElements: "
     209             :                          << keys_of(all_num_elements));
     210             :   }
     211             :   const auto& num_elements = all_num_elements.at(target_name);
     212             :   // Since some targets may not need all elements, we must specifically check if
     213             :   // we have all the necessary elements
     214             :   if (alg::all_of(num_elements, [&](const auto& element) {
     215             :         return vars_infos.at(temporal_id)
     216             :             .interpolation_is_done_for_these_elements.contains(element);
     217             :       })) {
     218             :     // Send data to InterpolationTarget, but only if the list of points is
     219             :     // non-empty.
     220             :     if (debug_print) {
     221             :       ss << "finished interpolation on all " << num_elements.size()
     222             :          << " elements on core " << Parallel::my_proc<int>(*cache);
     223             :     }
     224             : 
     225             :     if (not vars_infos.at(temporal_id).global_offsets.empty()) {
     226             :       const auto& info = vars_infos.at(temporal_id);
     227             :       auto& receiver_proxy = Parallel::get_parallel_component<
     228             :           InterpolationTarget<Metavariables, InterpolationTargetTag>>(*cache);
     229             :       Parallel::simple_action<
     230             :           Actions::InterpolationTargetReceiveVars<InterpolationTargetTag>>(
     231             :           receiver_proxy, info.vars, info.global_offsets, temporal_id);
     232             : 
     233             :       if (debug_print) {
     234             :         ss << " Sending interpolated data.\n";
     235             :       }
     236             :     } else if (debug_print) {
     237             :       ss << " Global offsets empty. Nothing to interpolate.";
     238             :     }
     239             : 
     240             :     // Clear interpolated data, since we don't need it anymore.
     241             :     db::mutate<Tags::InterpolatedVarsHolders<Metavariables>>(
     242             :         [&temporal_id](
     243             :             const gsl::not_null<
     244             :                 typename Tags::InterpolatedVarsHolders<Metavariables>::type*>
     245             :                 holders_l) {
     246             :           get<Vars::HolderTag<InterpolationTargetTag, Metavariables>>(
     247             :               *holders_l)
     248             :               .infos.erase(temporal_id);
     249             :         },
     250             :         box);
     251             :   } else if (debug_print) {
     252             :     ss << "interpolation not finished on all local elements of core "
     253             :        << Parallel::my_proc<int>(*cache) << ". Expected " << num_elements.size()
     254             :        << ", received "
     255             :        << vars_infos.at(temporal_id)
     256             :               .interpolation_is_done_for_these_elements.size();
     257             :   }
     258             : 
     259             :   if (debug_print) {
     260             :     Parallel::printf("%s\n", ss.str());
     261             :   }
     262             : }
     263             : 
     264             : }  // namespace intrp