Line data    Source code

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