21#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
22#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
25#include <opm/core/grid.h>
26#include <opm/common/ErrorMacros.hpp>
27#include <boost/any.hpp>
36#if HAVE_MPI && HAVE_DUNE_ISTL
38#include <opm/common/utility/platform_dependent/disable_warnings.h>
40#include <dune/istl/owneroverlapcopy.hh>
41#include <dune/common/parallel/interface.hh>
42#include <dune/common/parallel/communicator.hh>
43#include <dune/common/enumset.hh>
44#include <opm/common/utility/platform_dependent/reenable_warnings.h>
53 : std::integral_constant<bool, false>
55 template<
typename... T>
56 struct is_tuple<std::tuple<T...> >
57 : std::integral_constant<bool, true>
63class ParallelISTLInformation
67 typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
69 typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
72 ParallelISTLInformation()
73 : indexSet_(new ParallelIndexSet),
74 remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
75 communicator_(MPI_COMM_WORLD)
79 ParallelISTLInformation(MPI_Comm communicator)
80 : indexSet_(new ParallelIndexSet),
81 remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
82 communicator_(communicator)
88 ParallelISTLInformation(
const std::shared_ptr<ParallelIndexSet>& indexSet,
89 const std::shared_ptr<RemoteIndices>& remoteIndices,
90 MPI_Comm communicator)
91 : indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
96 ParallelISTLInformation(
const ParallelISTLInformation& other)
97 : indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
98 communicator_(other.communicator_)
101 std::shared_ptr<ParallelIndexSet> indexSet()
const
106 std::shared_ptr<RemoteIndices> remoteIndices()
const
108 return remoteIndices_;
111 Dune::CollectiveCommunication<MPI_Comm> communicator()
const
113 return communicator_;
118 void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices,
119 std::size_t local_component_size = 0, std::size_t num_components = 1)
const
121 ParallelIndexSet::GlobalIndex global_component_size = local_component_size;
122 if ( num_components > 1 )
124 ParallelIndexSet::GlobalIndex max_gi = 0;
126 for(
auto i = indexSet_->begin(), end = indexSet_->end(); i != end; ++i )
128 max_gi = std::max(max_gi, i->global());
130 global_component_size = max_gi+1;
131 global_component_size = communicator_.max(global_component_size);
133 indexSet.beginResize();
134 IndexSetInserter<ParallelIndexSet> inserter(indexSet, global_component_size,
135 local_component_size, num_components);
136 std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
137 indexSet.endResize();
138 remoteIndices.rebuild<
false>();
144 void copyOwnerToAll (
const T& source, T& dest)
const
146 typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
147 typedef Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner> OwnerSet;
148 typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
149 OwnerSet sourceFlags;
151 Dune::Interface interface(communicator_);
152 if( !remoteIndices_->isSynced() )
154 remoteIndices_->rebuild<
false>();
156 interface.build(*remoteIndices_,sourceFlags,destFlags);
157 Dune::BufferedCommunicator communicator;
158 communicator.template build<T>(interface);
159 communicator.template forward<CopyGatherScatter<T> >(source,dest);
163 const std::vector<double>& updateOwnerMask(
const T& container)
const
167 OPM_THROW(std::runtime_error,
"Trying to update owner mask without parallel information!");
169 if(
static_cast<std::size_t
>(container.size())!= ownerMask_.size() )
171 ownerMask_.resize(container.size(), 1.);
172 for(
auto i=indexSet_->begin(), end=indexSet_->end(); i!=end; ++i )
174 if (i->local().attribute()!=Dune::OwnerOverlapCopyAttributeSet::owner)
176 ownerMask_[i->local().local()] = 0.;
188 const std::vector<double>& getOwnerMask()
const
212 template<
typename Container,
typename BinaryOperator,
typename T>
213 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
216 computeReduction(container, binaryOperator, value, is_tuple<Container>());
222 template<
typename Container,
typename BinaryOperator,
typename T>
223 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
224 T& value, std::integral_constant<bool,true>)
const
226 computeTupleReduction(container, binaryOperator, value);
231 template<
typename Container,
typename BinaryOperator,
typename T>
232 void computeReduction(
const Container& container, BinaryOperator binaryOperator,
233 T& value, std::integral_constant<bool,false>)
const
235 std::tuple<const Container&> containers=std::tuple<const Container&>(container);
236 auto values=std::make_tuple(value);
237 auto operators=std::make_tuple(binaryOperator);
238 computeTupleReduction(containers, operators, values);
239 value=std::get<0>(values);
242 template<
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
243 void computeTupleReduction(
const std::tuple<Containers...>& containers,
244 std::tuple<BinaryOperators...>& operators,
245 std::tuple<ReturnValues...>& values)
const
247 static_assert(std::tuple_size<std::tuple<Containers...> >::value==
248 std::tuple_size<std::tuple<BinaryOperators...> >::value,
249 "We need the same number of containers and binary operators");
250 static_assert(std::tuple_size<std::tuple<Containers...> >::value==
251 std::tuple_size<std::tuple<ReturnValues...> >::value,
252 "We need the same number of containers and return values");
253 if( std::tuple_size<std::tuple<Containers...> >::value==0 )
258 std::tuple<ReturnValues...> init=values;
259 updateOwnerMask(std::get<0>(containers));
260 computeLocalReduction(containers, operators, values);
261 std::vector<std::tuple<ReturnValues...> > receivedValues(communicator_.size());
262 communicator_.allgather(&values, 1, &(receivedValues[0]));
264 for(
auto rvals=receivedValues.begin(), endvals=receivedValues.end(); rvals!=endvals;
267 computeGlobalReduction(*rvals, operators, values);
273 template<
int I=0,
typename... BinaryOperators,
typename... ReturnValues>
274 typename std::enable_if<I ==
sizeof...(BinaryOperators),
void>::type
275 computeGlobalReduction(
const std::tuple<ReturnValues...>&,
276 std::tuple<BinaryOperators...>&,
277 std::tuple<ReturnValues...>&)
const
280 template<
int I=0,
typename... BinaryOperators,
typename... ReturnValues>
281 typename std::enable_if<I !=
sizeof...(BinaryOperators),
void>::type
282 computeGlobalReduction(
const std::tuple<ReturnValues...>& receivedValues,
283 std::tuple<BinaryOperators...>& operators,
284 std::tuple<ReturnValues...>& values)
const
286 auto& val=std::get<I>(values);
287 val = std::get<I>(operators).localOperator()(val, std::get<I>(receivedValues));
288 computeGlobalReduction<I+1>(receivedValues, operators, values);
293 template<
int I=0,
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
294 typename std::enable_if<I==
sizeof...(Containers),
void>::type
295 computeLocalReduction(
const std::tuple<Containers...>&,
296 std::tuple<BinaryOperators...>&,
297 std::tuple<ReturnValues...>&)
const
300 template<
int I=0,
typename... Containers,
typename... BinaryOperators,
typename... ReturnValues>
301 typename std::enable_if<I!=
sizeof...(Containers),
void>::type
302 computeLocalReduction(
const std::tuple<Containers...>& containers,
303 std::tuple<BinaryOperators...>& operators,
304 std::tuple<ReturnValues...>& values)
const
306 const auto& container = std::get<I>(containers);
307 if( container.size() )
309 auto& reduceOperator = std::get<I>(operators);
316 auto mask = ownerMask_.begin();
317 auto& value = std::get<I>(values);
318 value = reduceOperator.getInitialValue();
320 for(
auto endVal=ownerMask_.end(); mask!=endVal;
323 value = reduceOperator(value, container[mask-ownerMask_.begin()], *mask);
326 computeLocalReduction<I+1>(containers, operators, values);
330 struct CopyGatherScatter
332 typedef typename Dune::CommPolicy<T>::IndexedType V;
334 static V gather(
const T& a, std::size_t i)
339 static void scatter(T& a, V v, std::size_t i)
345 class IndexSetInserter
348 typedef T ParallelIndexSet;
349 typedef typename ParallelIndexSet::LocalIndex LocalIndex;
350 typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
352 IndexSetInserter(ParallelIndexSet& indexSet,
const GlobalIndex& component_size,
353 std::size_t local_component_size, std::size_t num_components)
354 : indexSet_(&indexSet), component_size_(component_size),
355 local_component_size_(local_component_size),
356 num_components_(num_components)
358 void operator()(
const typename ParallelIndexSet::IndexPair& pair)
360 for(std::size_t i = 0; i < num_components_; i++)
361 indexSet_->add(i * component_size_ + pair.global(),
362 LocalIndex(i * local_component_size_ + pair.local(),
363 pair.local().attribute()));
366 ParallelIndexSet* indexSet_;
368 GlobalIndex component_size_;
370 std::size_t local_component_size_;
372 std::size_t num_components_;
374 std::shared_ptr<ParallelIndexSet> indexSet_;
375 std::shared_ptr<RemoteIndices> remoteIndices_;
376 Dune::CollectiveCommunication<MPI_Comm> communicator_;
377 mutable std::vector<double> ownerMask_;
387 template<
typename BinaryOperator>
388 struct MaskIDOperator
392 typedef typename std::remove_cv<
393 typename std::remove_reference<typename BinaryOperator::result_type>::type
401 template<
class T,
class T1>
402 T operator()(
const T& t1,
const T& t2,
const T1& mask)
404 return b_(t1, maskValue(t2, mask));
406 template<
class T,
class T1>
407 T maskValue(
const T& t,
const T1& mask)
411 BinaryOperator& localOperator()
415 Result getInitialValue()
425 struct InnerProductFunctor
434 T operator()(
const T& t1,
const T& t2,
const T1& mask)
436 T masked = maskValue(t2, mask);
437 return t1 + masked * masked;
440 T maskValue(
const T& t,
const T1& mask)
444 std::plus<T> localOperator()
446 return std::plus<T>();
459 template<
typename BinaryOperator>
460 struct MaskToMinOperator
464 typedef typename std::remove_reference<
465 typename std::remove_const<typename BinaryOperator::result_type>::type
468 MaskToMinOperator(BinaryOperator b)
477 template<
class T,
class T1>
478 T operator()(
const T& t1,
const T& t2,
const T1& mask)
480 return b_(t1, maskValue(t2, mask));
482 template<
class T,
class T1>
483 T maskValue(
const T& t,
const T1& mask)
491 return getInitialValue();
494 Result getInitialValue()
499 if( std::is_integral<Result>::value )
501 return std::numeric_limits<Result>::min();
505 return -std::numeric_limits<Result>::max();
512 BinaryOperator& localOperator()
523 template<
typename BinaryOperator>
524 struct MaskToMaxOperator
528 typedef typename std::remove_cv<
529 typename std::remove_reference<typename BinaryOperator::result_type>::type
532 MaskToMaxOperator(BinaryOperator b)
541 template<
class T,
class T1>
542 T operator()(
const T& t1,
const T& t2,
const T1& mask)
544 return b_(t1, maskValue(t2, mask));
546 template<
class T,
class T1>
547 T maskValue(
const T& t,
const T1& mask)
555 return std::numeric_limits<T>::max();
558 BinaryOperator& localOperator()
562 Result getInitialValue()
564 return std::numeric_limits<Result>::max();
573 MaskIDOperator<std::plus<T> >
574 makeGlobalSumFunctor()
576 return MaskIDOperator<std::plus<T> >();
582 MaskToMinOperator<std::pointer_to_binary_function<const T&,const T&,const T&> >
583 makeGlobalMaxFunctor()
585 return MaskToMinOperator<std::pointer_to_binary_function<const T&,const T&,const T&> >
586 (std::pointer_to_binary_function<const T&,const T&,const T&>
587 ((
const T&(*)(
const T&,
const T&))std::max<T>));
593 template<
typename T,
typename Enable =
void>
596 using result_type = T;
597 result_type operator()(
const T& t1,
600 return std::max(std::abs(t1), std::abs(t2));
608 struct MaxAbsFunctor<T, typename std::enable_if<std::is_unsigned<T>::value>::type>
610 using result_type = T;
611 result_type operator()(
const T& t1,
614 return std::max(t1, t2);
623 MaskIDOperator<detail::MaxAbsFunctor<T> >
624 makeLInfinityNormFunctor()
626 return MaskIDOperator<detail::MaxAbsFunctor<T> >();
632 MaskToMaxOperator<std::pointer_to_binary_function<const T&,const T&,const T&> >
633 makeGlobalMinFunctor()
635 return MaskToMaxOperator<std::pointer_to_binary_function<const T&,const T&,const T&> >
636 (std::pointer_to_binary_function<const T&,const T&,const T&>
637 ((
const T&(*)(
const T&,
const T&))std::min<T>));
640 InnerProductFunctor<T>
641 makeInnerProductFunctor()
643 return InnerProductFunctor<T>();
664 (void)anyComm; (void)grid;
676 ->
decltype(container[0]*(*maskContainer)[0])
678 decltype(container[0]*(*maskContainer)[0]) initial = 0;
682 return std::inner_product(container.begin(), container.end(), maskContainer->begin(),
686 return std::accumulate(container.begin(), container.end(), initial);
Definition: AnisotropicEikonal.hpp:44
void extractParallelGridInformationToISTL(boost::any &anyComm, const UnstructuredGrid &grid)
Extracts the information about the data decomposition from the grid for dune-istl.
Definition: ParallelIstlInformation.hpp:662
auto accumulateMaskedValues(const T1 &container, const std::vector< double > *maskContainer) -> decltype(container[0] *(*maskContainer)[0])
Accumulates entries masked with 1.
Definition: ParallelIstlInformation.hpp:675