3 #ifndef DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH 4 #define DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH 11 #include <dune/common/parallel/mpihelper.hh> 12 #include <dune/common/stdstreams.hh> 13 #include <dune/common/typetraits.hh> 15 #if HAVE_UG && PDELAB_SEQUENTIAL_UG 17 #include <dune/grid/uggrid.hh> 20 #include <dune/istl/owneroverlapcopy.hh> 21 #include <dune/istl/solvercategory.hh> 22 #include <dune/istl/operators.hh> 23 #include <dune/istl/solvers.hh> 24 #include <dune/istl/preconditioners.hh> 25 #include <dune/istl/scalarproducts.hh> 26 #include <dune/istl/paamg/amg.hh> 27 #include <dune/istl/paamg/pinfo.hh> 28 #include <dune/istl/io.hh> 29 #include <dune/istl/superlu.hh> 52 template<
typename GFS>
57 typedef int RankIndex;
65 typedef typename GFS::Ordering::Traits::ContainerIndex ContainerIndex;
71 , _rank(gfs.gridView().comm().
rank())
80 if (gfs.entitySet().partitions().value == Partitions::interiorBorder.value)
84 _interiorBorder_all_interface = InteriorBorder_InteriorBorder_Interface;
85 _all_all_interface = InteriorBorder_InteriorBorder_Interface;
90 _interiorBorder_all_interface = InteriorBorder_All_Interface;
91 _all_all_interface = All_All_Interface;
94 if (_gfs.gridView().comm().size()>1)
100 gdh(_gfs,_ghosts,
false);
101 _gfs.gridView().communicate(gdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
107 _gfs.gridView().communicate(pdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
125 template<
typename X,
typename Mask>
128 typename Mask::const_iterator mask_it = mask.begin();
129 for (
typename X::iterator it = x.begin(),
137 template<
typename X,
typename Mask>
140 typename Mask::const_iterator mask_it = mask.begin();
141 for (
typename X::iterator it = x.begin(),
145 *it = (*mask_it == _rank ? *it :
typename X::field_type(0));
151 bool owned(
const ContainerIndex& i)
const 153 return _ranks[i] == _rank;
163 template<
typename X,
typename Y>
164 typename PromotionTraits<
165 typename X::field_type,
166 typename Y::field_type
182 template<
typename X,
typename Y,
typename Mask>
183 typename PromotionTraits<
184 typename X::field_type,
185 typename Y::field_type
189 typedef typename PromotionTraits<
190 typename X::field_type,
191 typename Y::field_type
192 >::PromotedType result_type;
196 typename Y::const_iterator y_it = y.begin();
197 typename Mask::const_iterator mask_it = mask.begin();
198 for (
typename X::const_iterator x_it = x.begin(),
201 ++x_it, ++y_it, ++mask_it)
209 template<
typename X,
typename Y,
typename Mask>
210 typename PromotionTraits<
211 typename X::field_type,
212 typename Y::field_type
216 typedef typename PromotionTraits<
217 typename X::field_type,
218 typename Y::field_type
219 >::PromotedType result_type;
223 typename Y::const_iterator y_it = y.begin();
224 typename Mask::const_iterator mask_it = mask.begin();
225 for (
typename X::const_iterator x_it = x.begin(),
228 ++x_it, ++y_it, ++mask_it)
229 r += (*mask_it == _rank ? Dune::dot(*x_it,*y_it) : result_type(0));
261 template<
typename MatrixType,
typename Comm>
269 bool owned_for_amg(std::size_t i)
const 279 const RankIndex _rank;
285 InterfaceType _interiorBorder_all_interface;
288 InterfaceType _all_all_interface;
293 template<
typename GFS>
294 template<
typename M,
typename C>
300 const bool is_bcrs_matrix =
308 const bool block_type_is_field_matrix =
317 static_assert(is_bcrs_matrix && block_type_is_field_matrix,
"matrix structure not compatible with AMG");
326 typedef typename GFS::Traits::GridViewType GV;
327 typedef typename RankVector::size_type size_type;
328 const GV& gv = _gfs.gridView();
331 const bool need_communication = _gfs.gridView().comm().size() > 1;
335 BoolVector sharedDOF(_gfs,
false);
337 if (need_communication)
340 _gfs.gridView().communicate(data_handle,_all_all_interface,Dune::ForwardCommunication);
344 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
345 GlobalIndex count = 0;
347 for (size_type i = 0; i < sharedDOF.N(); ++i)
348 if (owned_for_amg(i) &&
native(sharedDOF)[i][0])
351 dverb << gv.comm().rank() <<
": shared block count is " << count.touint() << std::endl;
354 std::vector<GlobalIndex> counts(_gfs.gridView().comm().size());
355 _gfs.gridView().comm().allgather(&count, 1, &(counts[0]));
358 GlobalIndex start = std::accumulate(counts.begin(),counts.begin() + _rank,GlobalIndex(0));
361 GIVector scalarIndices(_gfs, std::numeric_limits<GlobalIndex>::max());
363 for (size_type i = 0; i < sharedDOF.N(); ++i)
364 if (owned_for_amg(i) &&
native(sharedDOF)[i][0])
366 native(scalarIndices)[i][0] = start;
371 if (need_communication)
374 _gfs.gridView().communicate(data_handle,_interiorBorder_all_interface,Dune::ForwardCommunication);
378 c.indexSet().beginResize();
379 for (size_type i=0; i<scalarIndices.N(); ++i)
381 Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
382 if(
native(scalarIndices)[i][0] != std::numeric_limits<GlobalIndex>::max())
385 if (owned_for_amg(i))
388 attr = Dune::OwnerOverlapCopyAttributeSet::owner;
392 attr = Dune::OwnerOverlapCopyAttributeSet::copy;
394 c.indexSet().add(
native(scalarIndices)[i][0],
typename C::ParallelIndexSet::LocalIndex(i,attr));
397 c.indexSet().endResize();
400 std::set<int> neighbors;
402 if (need_communication)
405 _gfs.gridView().communicate(data_handle,_all_all_interface,Dune::ForwardCommunication);
408 c.remoteIndices().setNeighbours(neighbors);
409 c.remoteIndices().template rebuild<false>();
414 template<
int s,
bool isFakeMPIHelper>
417 typedef Dune::Amg::SequentialInformation
type;
427 typedef OwnerOverlapCopyCommunication<bigunsignedint<s>,
int>
type;
436 #if HAVE_UG && PDELAB_SEQUENTIAL_UG 438 void assertParallelUG(Dune::CollectiveCommunication<Dune::UGGrid<dim> > comm)
440 static_assert(Dune::AlwaysFalse<Dune::UGGrid<dim> >::
value,
"Using sequential UG in parallel environment");
449 #endif // DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
Extracts the container tag from T.
Definition: backend/istl/tags.hh:145
Tag describing an arbitrary FieldVector.
Definition: backend/istl/tags.hh:46
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
const std::string s
Definition: function.hh:830
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper...
Definition: backend/interface.hh:176
RankIndex rank() const
Returns the MPI rank of this process.
Definition: parallelhelper.hh:237
void assertParallelUG(T comm)
Definition: parallelhelper.hh:433
Definition: genericdatahandle.hh:716
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:249
bool isGhost(const ContainerIndex &i) const
Tests whether the given index belongs to a ghost DOF.
Definition: parallelhelper.hh:157
void createIndexSetAndProjectForAMG(MatrixType &m, Comm &c)
Makes the matrix consistent and creates the parallel information for AMG.
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
OwnerOverlapCopyCommunication< bigunsignedint< s >, int > type
Definition: parallelhelper.hh:427
Tag describing a BCRSMatrix.
Definition: backend/istl/tags.hh:63
Data handle for marking ghost DOFs.
Definition: genericdatahandle.hh:812
Tag describing an arbitrary FieldMatrix.
Definition: backend/istl/tags.hh:83
bool owned(const ContainerIndex &i) const
Tests whether the given index is owned by this process.
Definition: parallelhelper.hh:151
Data handle for collecting set of neighboring MPI ranks.
Definition: genericdatahandle.hh:1060
ParallelHelper(const GFS &gfs, int verbose=1)
Definition: parallelhelper.hh:69
GatherScatter data handle for creating a disjoint DOF partitioning.
Definition: genericdatahandle.hh:930
void maskForeignDOFs(X &x) const
Mask out all DOFs not owned by the current process with 0.
Definition: parallelhelper.hh:115
PromotionTraits< typename X::field_type, typename Y::field_type >::PromotedType disjointDot(const X &x, const Y &y) const
Calculates the (rank-local) dot product of x and y on the disjoint partition defined by the helper...
Definition: parallelhelper.hh:168
Definition: parallelhelper.hh:415
Data handle for marking shared DOFs.
Definition: genericdatahandle.hh:1014
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
Tag describing a BlockVector.
Definition: backend/istl/tags.hh:26
Definition: parallelhelper.hh:53