3 #ifndef DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
8 #include <dune/common/deprecated.hh>
9 #include <dune/common/parallel/mpihelper.hh>
10 #include <dune/common/static_assert.hh>
11 #include <dune/common/stdstreams.hh>
13 #include <dune/istl/owneroverlapcopy.hh>
14 #include <dune/istl/solvercategory.hh>
15 #include <dune/istl/operators.hh>
16 #include <dune/istl/solvers.hh>
17 #include <dune/istl/preconditioners.hh>
18 #include <dune/istl/scalarproducts.hh>
19 #include <dune/istl/paamg/amg.hh>
20 #include <dune/istl/paamg/pinfo.hh>
21 #include <dune/istl/io.hh>
22 #include <dune/istl/superlu.hh>
44 template<
typename GFS>
49 typedef int RankIndex;
57 typedef typename GFS::Ordering::Traits::ContainerIndex ContainerIndex;
63 , _rank(gfs.gridView().comm().
rank())
71 if (gfs.ordering().containedPartitions() == NonOverlappingPartitionSelector::partition_mask)
75 _interiorBorder_all_interface = InteriorBorder_InteriorBorder_Interface;
76 _all_all_interface = InteriorBorder_InteriorBorder_Interface;
81 _interiorBorder_all_interface = InteriorBorder_All_Interface;
82 _all_all_interface = All_All_Interface;
85 if (_gfs.gridView().comm().size()>1)
91 gdh(_gfs,_ghosts,
false);
92 _gfs.gridView().communicate(gdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
98 _gfs.gridView().communicate(pdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
115 template<
typename X,
typename Mask>
118 typename Mask::const_iterator mask_it = mask.begin();
119 for (
typename X::iterator it = x.begin(),
127 template<
typename X,
typename Mask>
130 typename Mask::const_iterator mask_it = mask.begin();
131 for (
typename X::iterator it = x.begin(),
135 *it = (*mask_it == _rank ? *it :
typename X::field_type(0));
141 bool owned(
const ContainerIndex& i)
const
143 return _ranks[i] == _rank;
153 template<
typename X,
typename Y>
154 typename PromotionTraits<
155 typename X::field_type,
156 typename Y::field_type
171 template<
typename X,
typename Y,
typename Mask>
172 typename PromotionTraits<
173 typename X::field_type,
174 typename Y::field_type
178 typedef typename PromotionTraits<
179 typename X::field_type,
180 typename Y::field_type
181 >::PromotedType result_type;
185 typename Y::const_iterator y_it = y.begin();
186 typename Mask::const_iterator mask_it = mask.begin();
187 for (
typename X::const_iterator x_it = x.begin(),
190 ++x_it, ++y_it, ++mask_it)
198 template<
typename X,
typename Y,
typename Mask>
199 typename PromotionTraits<
200 typename X::field_type,
201 typename Y::field_type
205 typedef typename PromotionTraits<
206 typename X::field_type,
207 typename Y::field_type
208 >::PromotedType result_type;
212 typename Y::const_iterator y_it = y.begin();
213 typename Mask::const_iterator mask_it = mask.begin();
214 for (
typename X::const_iterator x_it = x.begin(),
217 ++x_it, ++y_it, ++mask_it)
218 r += (*mask_it == _rank ? Dune::dot(*x_it,*y_it) : result_type(0));
250 template<
typename MatrixType,
typename Comm>
258 bool owned_for_amg(std::size_t i)
const
260 return _ranks.base()[i][0] == _rank;
266 bool is_ghost_for_amg(std::size_t i)
const
268 return _ghosts.base()[i][0];
276 const RankIndex _rank;
282 InterfaceType _interiorBorder_all_interface;
285 InterfaceType _all_all_interface;
290 template<
typename GFS>
291 template<
typename M,
typename C>
295 const bool is_bcrs_matrix =
303 const bool block_type_is_field_matrix =
312 dune_static_assert(is_bcrs_matrix && block_type_is_field_matrix,
"matrix structure not compatible with AMG");
321 typedef typename GFS::Traits::GridViewType GV;
322 typedef typename RankVector::size_type size_type;
323 const GV& gv = _gfs.gridView();
326 const bool need_communication = _gfs.gridView().comm().size() > 1;
330 BoolVector sharedDOF(_gfs,
false);
332 if (need_communication)
335 _gfs.gridView().communicate(data_handle,_all_all_interface,Dune::ForwardCommunication);
339 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
340 GlobalIndex count = 0;
342 for (size_type i = 0; i < sharedDOF.N(); ++i)
343 if (owned_for_amg(i) && sharedDOF.base()[i][0])
346 dverb << gv.comm().rank() <<
": shared block count is " << count.touint() << std::endl;
349 std::vector<GlobalIndex> counts(_gfs.gridView().comm().size());
350 _gfs.gridView().comm().allgather(&count, 1, &(counts[0]));
353 GlobalIndex start = std::accumulate(counts.begin(),counts.begin() + _rank,GlobalIndex(0));
356 GIVector scalarIndices(_gfs, std::numeric_limits<GlobalIndex>::max());
358 for (size_type i = 0; i < sharedDOF.N(); ++i)
359 if (owned_for_amg(i) && sharedDOF.base()[i][0])
361 scalarIndices.base()[i][0] = start;
366 if (need_communication)
369 _gfs.gridView().communicate(data_handle,_interiorBorder_all_interface,Dune::ForwardCommunication);
373 c.indexSet().beginResize();
374 for (size_type i=0; i<scalarIndices.N(); ++i)
376 Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
377 if(scalarIndices.base()[i][0] != std::numeric_limits<GlobalIndex>::max())
380 if (owned_for_amg(i))
383 attr = Dune::OwnerOverlapCopyAttributeSet::owner;
385 else if (is_ghost_for_amg(i) && c.getSolverCategory() ==
static_cast<int>(SolverCategory::nonoverlapping))
388 attr = Dune::OwnerOverlapCopyAttributeSet::overlap;
392 attr = Dune::OwnerOverlapCopyAttributeSet::copy;
394 c.indexSet().add(scalarIndices.base()[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;
438 #endif // DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
Extracts the container tag from T.
Definition: backend/istl/tags.hh:143
void maskForeignDOFs(X &x) const
Mask out all DOFs not owned by the current process with 0.
Definition: parallelhelper.hh:106
Returns the raw ISTL type associated with C, or C itself it is already an ISTL type.
Definition: backend/istl/utility.hh:40
Tag describing a BCRSMatrix.
Definition: backend/istl/tags.hh:61
Data handle for marking shared DOFs.
Definition: genericdatahandle.hh:1015
ParallelHelper(const GFS &gfs, int verbose=1)
Definition: parallelhelper.hh:61
Tag describing a BlockVector.
Definition: backend/istl/tags.hh:24
bool isGhost(const ContainerIndex &i) const
Tests whether the given index belongs to a ghost DOF.
Definition: parallelhelper.hh:147
Definition: parallelhelper.hh:415
GatherScatter data handle for creating a disjoint DOF partitioning.
Definition: genericdatahandle.hh:931
bool owned(const ContainerIndex &i) const
Tests whether the given index is owned by this process.
Definition: parallelhelper.hh:141
C type
Definition: backend/istl/utility.hh:42
V & raw(V &v)
Returns the raw ISTL object associated with v, or v itself it is already an ISTL object.
Definition: backend/istl/utility.hh:26
Definition: parallelhelper.hh:45
Definition: backendselector.hh:11
void createIndexSetAndProjectForAMG(MatrixType &m, Comm &c)
Makes the matrix consistent and creates the parallel information for AMG.
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
BackendVectorSelectorHelper< Backend, GridFunctionSpace, FieldType >::Type Type
Definition: backendselector.hh:14
Definition: genericdatahandle.hh:717
Data handle for marking ghost DOFs.
Definition: genericdatahandle.hh:813
Tag describing an arbitrary FieldVector.
Definition: backend/istl/tags.hh:44
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:158
tags::container< T >::type container_tag(const T &)
Gets instance of container tag associated with T.
Definition: backend/istl/tags.hh:247
OwnerOverlapCopyCommunication< bigunsignedint< s >, int > type
Definition: parallelhelper.hh:427
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:417
RankIndex rank() const
Returns the MPI rank of this process.
Definition: parallelhelper.hh:226
Data handle for collecting set of neighboring MPI ranks.
Definition: genericdatahandle.hh:1061
const std::string s
Definition: function.hh:1103
Tag describing an arbitrary FieldMatrix.
Definition: backend/istl/tags.hh:81