21 #ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
22 #define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
25 #include <opm/grid/UnstructuredGrid.h>
26 #include <opm/common/ErrorMacros.hpp>
34 #include <type_traits>
37 #if HAVE_MPI && HAVE_DUNE_ISTL
39 #include <opm/common/utility/platform_dependent/disable_warnings.h>
41 #include <dune/istl/owneroverlapcopy.hh>
42 #include <dune/common/parallel/interface.hh>
43 #include <dune/common/parallel/communicator.hh>
44 #include <dune/common/enumset.hh>
45 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
47 #include <opm/simulators/utils/ParallelCommunication.hpp>
56 : std::integral_constant<bool, false>
58 template<
typename... T>
59 struct is_tuple<std::tuple<T...> >
60 : std::integral_constant<bool, true>
66 class ParallelISTLInformation
70 typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
72 typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
75 ParallelISTLInformation()
76 : indexSet_(new ParallelIndexSet),
77 remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
78 communicator_(MPI_COMM_WORLD)
82 ParallelISTLInformation(MPI_Comm communicator)
83 : indexSet_(new ParallelIndexSet),
84 remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
85 communicator_(communicator)
91 ParallelISTLInformation(
const std::shared_ptr<ParallelIndexSet>& indexSet,
92 const std::shared_ptr<RemoteIndices>& remoteIndices,
93 MPI_Comm communicator)
94 : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
99 ParallelISTLInformation(
const ParallelISTLInformation& other)
100 : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
101 communicator_(other.communicator_)
104 std::shared_ptr<ParallelIndexSet> indexSet()
const
109 std::shared_ptr<RemoteIndices> remoteIndices()
const
111 return remoteIndices_;
114 Parallel::Communication communicator()
const
116 return communicator_;
121 void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices,
122 std::size_t local_component_size = 0, std::size_t num_components = 1)
const
124 ParallelIndexSet::GlobalIndex global_component_size = local_component_size;
125 if ( num_components > 1 )
127 ParallelIndexSet::GlobalIndex max_gi = 0;
129 for(
auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
131 max_gi = std::max(max_gi, i->global());
133 global_component_size = max_gi+1;
134 global_component_size = communicator_.max(global_component_size);
136 indexSet.beginResize();
137 IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
138 local_component_size, num_components);
139 std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
140 indexSet.endResize();
141 remoteIndices.rebuild<
false>();
147 void copyOwnerToAll (
const T& source, T& dest)
const
149 typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
150 typedef Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner> OwnerSet;
151 typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
152 OwnerSet sourceFlags;
154 Dune::Interface interface(communicator_);
155 if( !remoteIndices_->isSynced() )
157 remoteIndices_->rebuild<
false>();
159 interface.build(*remoteIndices_,sourceFlags,destFlags);
160 Dune::BufferedCommunicator communicator;
161 communicator.template build<T>(interface);
162 communicator.template forward<CopyGatherScatter<T> >(source,dest);
166 const std::vector<double>& updateOwnerMask(
const T& container)
const
170 OPM_THROW(std::runtime_error,
"Trying to update owner mask without parallel information!");
172 if(
static_cast<std::size_t
>(container.size())!= ownerMask_.size() )
174 ownerMask_.resize(container.size(), 1.);
175 for(
auto i=indexSet_->begin(), end=indexSet_->end(); i!=end; ++i )
177 if (i->local().attribute()!=Dune::OwnerOverlapCopyAttributeSet::owner)
179 ownerMask_[i->local().local()] = 0.;
191 const std::vector<double>& getOwnerMask()
const
215 template<
typename Container,
typename BinaryOperator,
typename T>
216 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
219 computeReduction(container, binaryOperator, value, is_tuple<Container>());
225 template<
typename Container,
typename BinaryOperator,
typename T>
226 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
227 T& value, std::integral_constant<bool,true>)
const
229 computeTupleReduction(container, binaryOperator, value);
234 template<
typename Container,
typename BinaryOperator,
typename T>
235 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
236 T& value, std::integral_constant<bool,false>)
const
238 std::tuple<const Container&> containers=std::tuple<const Container&>(container);
239 auto values=std::make_tuple(value);
240 auto operators=std::make_tuple(binaryOperator);
241 computeTupleReduction(containers, operators, values);
242 value=std::get<0>(values);
245 template<
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
246 void computeTupleReduction(
const std::tuple<Containers...>& containers,
247 std::tuple<BinaryOperators...>& operators,
248 std::tuple<ReturnValues...>& values)
const
250 static_assert(std::tuple_size<std::tuple<Containers...> >::value==
251 std::tuple_size<std::tuple<BinaryOperators...> >::value,
252 "We need the same number of containers and binary operators");
253 static_assert(std::tuple_size<std::tuple<Containers...> >::value==
254 std::tuple_size<std::tuple<ReturnValues...> >::value,
255 "We need the same number of containers and return values");
256 if( std::tuple_size<std::tuple<Containers...> >::value==0 )
261 std::tuple<ReturnValues...> init=values;
262 updateOwnerMask(std::get<0>(containers));
263 computeLocalReduction(containers, operators, values);
264 std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
265 communicator_.allgather(&values, 1, &(receivedValues[0]));
267 for(
auto rvals=receivedValues.begin(), endvals=receivedValues.end(); rvals!=endvals;
270 computeGlobalReduction(*rvals, operators, values);
276 template<
int I=0,
typename... BinaryOperators,
typename... ReturnValues>
277 typename std::enable_if<I ==
sizeof...(BinaryOperators),
void>::type
278 computeGlobalReduction(
const std::tuple<ReturnValues...>&,
279 std::tuple<BinaryOperators...>&,
280 std::tuple<ReturnValues...>&)
const
283 template<
int I=0,
typename... BinaryOperators,
typename... ReturnValues>
284 typename std::enable_if<I !=
sizeof...(BinaryOperators),
void>::type
285 computeGlobalReduction(
const std::tuple<ReturnValues...>& receivedValues,
286 std::tuple<BinaryOperators...>& operators,
287 std::tuple<ReturnValues...>& values)
const
289 auto& val=std::get<I>(values);
290 val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
291 computeGlobalReduction<I+1>(receivedValues, operators, values);
296 template<
int I=0,
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
297 typename std::enable_if<I==
sizeof...(Containers),
void>::type
298 computeLocalReduction(
const std::tuple<Containers...>&,
299 std::tuple<BinaryOperators...>&,
300 std::tuple<ReturnValues...>&)
const
303 template<
int I=0,
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
304 typename std::enable_if<I!=
sizeof...(Containers),
void>::type
305 computeLocalReduction(
const std::tuple<Containers...>& containers,
306 std::tuple<BinaryOperators...>& operators,
307 std::tuple<ReturnValues...>& values)
const
309 const auto& container = std::get<I>(containers);
310 if( container.size() )
312 auto& reduceOperator = std::get<I>(operators);
319 auto mask = ownerMask_.begin();
320 auto& value = std::get<I>(values);
321 value = reduceOperator.getInitialValue();
323 for(
auto endVal=ownerMask_.end(); mask!=endVal;
326 value = reduceOperator(value, container[mask-ownerMask_.begin()], *mask);
329 computeLocalReduction<I+1>(containers, operators, values);
333 struct CopyGatherScatter
335 typedef typename Dune::CommPolicy<T>::IndexedType V;
337 static V gather(
const T& a, std::size_t i)
342 static void scatter(T& a, V v, std::size_t i)
348 class IndexSetInserter
351 typedef T ParallelIndexSet;
352 typedef typename ParallelIndexSet::LocalIndex LocalIndex;
353 typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
355 IndexSetInserter(ParallelIndexSet& indexSet,
const GlobalIndex& component_size,
356 std::size_t local_component_size, std::size_t num_components)
357 : indexSet_(&indexSet), component_size_(component_size),
358 local_component_size_(local_component_size),
359 num_components_(num_components)
361 void operator()(
const typename ParallelIndexSet::IndexPair& pair)
363 for(std::size_t i = 0; i < num_components_; i++)
364 indexSet_->add(i * component_size_ + pair.global(),
365 LocalIndex(i * local_component_size_ + pair.local(),
366 pair.local().attribute()));
369 ParallelIndexSet* indexSet_;
371 GlobalIndex component_size_;
373 std::size_t local_component_size_;
375 std::size_t num_components_;
377 std::shared_ptr<ParallelIndexSet> indexSet_;
378 std::shared_ptr<RemoteIndices> remoteIndices_;
379 Dune::CollectiveCommunication<MPI_Comm> communicator_;
380 mutable std::vector<double> ownerMask_;
390 template<
typename BinaryOperator>
391 struct MaskIDOperator
395 typedef typename std::remove_cv<
396 typename std::remove_reference<typename BinaryOperator::result_type>::type
404 template<
class T,
class T1>
405 T operator()(
const T& t1,
const T& t2,
const T1& mask)
407 return b_(t1, maskValue(t2, mask));
409 template<
class T,
class T1>
410 T maskValue(
const T& t,
const T1& mask)
414 BinaryOperator& localOperator()
418 Result getInitialValue()
428 struct InnerProductFunctor
437 T operator()(
const T& t1,
const T& t2,
const T1& mask)
439 T masked = maskValue(t2, mask);
440 return t1 + masked * masked;
443 T maskValue(
const T& t,
const T1& mask)
447 std::plus<T> localOperator()
449 return std::plus<T>();
462 template<
typename BinaryOperator>
463 struct MaskToMinOperator
467 typedef typename std::remove_reference<
468 typename std::remove_const<typename BinaryOperator::result_type>::type
471 MaskToMinOperator(BinaryOperator b)
480 template<
class T,
class T1>
481 T operator()(
const T& t1,
const T& t2,
const T1& mask)
483 return b_(t1, maskValue(t2, mask));
485 template<
class T,
class T1>
486 T maskValue(
const T& t,
const T1& mask)
494 return getInitialValue();
497 Result getInitialValue()
502 if( std::is_integral<Result>::value )
504 return std::numeric_limits<Result>::min();
508 return -std::numeric_limits<Result>::max();
515 BinaryOperator& localOperator()
526 template<
typename BinaryOperator>
527 struct MaskToMaxOperator
531 typedef typename std::remove_cv<
532 typename std::remove_reference<typename BinaryOperator::result_type>::type
535 MaskToMaxOperator(BinaryOperator b)
544 template<
class T,
class T1>
545 T operator()(
const T& t1,
const T& t2,
const T1& mask)
547 return b_(t1, maskValue(t2, mask));
549 template<
class T,
class T1>
550 T maskValue(
const T& t,
const T1& mask)
558 return std::numeric_limits<T>::max();
561 BinaryOperator& localOperator()
565 Result getInitialValue()
567 return std::numeric_limits<Result>::max();
576 MaskIDOperator<std::plus<T> >
577 makeGlobalSumFunctor()
579 return MaskIDOperator<std::plus<T> >();
585 auto makeGlobalMaxFunctor()
589 using result_type = T;
590 const result_type& operator()(
const T& t1,
const T& t2)
592 return std::max(t1, t2);
595 return MaskToMinOperator(MaxOp());
601 template<
typename T,
typename Enable =
void>
604 using result_type = T;
605 result_type operator()(
const T& t1,
608 return std::max(std::abs(t1), std::abs(t2));
616 struct MaxAbsFunctor<T, typename std::enable_if<std::is_unsigned<T>::value>::type>
618 using result_type = T;
619 result_type operator()(
const T& t1,
622 return std::max(t1, t2);
631 MaskIDOperator<detail::MaxAbsFunctor<T> >
632 makeLInfinityNormFunctor()
634 return MaskIDOperator<detail::MaxAbsFunctor<T> >();
641 makeGlobalMinFunctor()
645 using result_type = T;
646 const result_type& operator()(
const T& t1,
const T& t2)
648 return std::min(t1, t2);
651 return MaskToMaxOperator(MinOp());
654 InnerProductFunctor<T>
655 makeInnerProductFunctor()
657 return InnerProductFunctor<T>();
678 (void)anyComm; (void)grid;
690 -> decltype(container[0]*(*maskContainer)[0])
692 decltype(container[0]*(*maskContainer)[0]) initial = 0;
696 return std::inner_product(container.begin(), container.end(), maskContainer->begin(),
700 return std::accumulate(container.begin(), container.end(), initial);
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: BlackoilPhases.hpp:27
void extractParallelGridInformationToISTL(std::any &anyComm, const UnstructuredGrid &grid)
Extracts the information about the data decomposition from the grid for dune-istl.
Definition: ParallelIstlInformation.hpp:676
auto accumulateMaskedValues(const T1 &container, const std::vector< double > *maskContainer) -> decltype(container[0] *(*maskContainer)[0])
Accumulates entries masked with 1.
Definition: ParallelIstlInformation.hpp:689