3 #ifndef DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
4 #define DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
8 #include <dune/common/parallel/mpihelper.hh>
9 #include <dune/common/stdstreams.hh>
10 #include <dune/common/typetraits.hh>
12 #if HAVE_UG && PDELAB_SEQUENTIAL_UG
14 #include <dune/grid/uggrid.hh>
17 #include <dune/istl/owneroverlapcopy.hh>
18 #include <dune/istl/solvercategory.hh>
19 #include <dune/istl/operators.hh>
20 #include <dune/istl/solvers.hh>
21 #include <dune/istl/preconditioners.hh>
22 #include <dune/istl/scalarproducts.hh>
23 #include <dune/istl/paamg/amg.hh>
24 #include <dune/istl/paamg/pinfo.hh>
25 #include <dune/istl/io.hh>
26 #include <dune/istl/superlu.hh>
49 template<
typename GFS>
54 typedef int RankIndex;
62 typedef typename GFS::Ordering::Traits::ContainerIndex ContainerIndex;
68 , _rank(gfs.gridView().comm().
rank())
69 , _rank_partition(gfs,_rank)
77 if (gfs.entitySet().partitions().value == Partitions::interiorBorder.value)
81 _interiorBorder_all_interface = InteriorBorder_InteriorBorder_Interface;
82 _all_all_interface = InteriorBorder_InteriorBorder_Interface;
87 _interiorBorder_all_interface = InteriorBorder_All_Interface;
88 _all_all_interface = All_All_Interface;
91 if (_gfs.gridView().comm().size()>1)
97 gdh(_gfs,_ghosts,
false);
98 _gfs.gridView().communicate(gdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
104 _gfs.gridView().communicate(pdh,_interiorBorder_all_interface,Dune::ForwardCommunication);
109 std::set<RankIndex> rank_set;
110 for (RankIndex
rank : _rank_partition)
112 rank_set.insert(
rank);
114 for (RankIndex
rank : rank_set)
115 _neighbor_ranks.push_back(
rank);
121 return _neighbor_ranks;
136 template<
typename X,
typename Mask>
139 typename Mask::const_iterator mask_it = mask.begin();
140 for (
typename X::iterator it = x.begin(),
148 template<
typename X,
typename Mask>
151 typename Mask::const_iterator mask_it = mask.begin();
152 for (
typename X::iterator it = x.begin(),
156 *it = (*mask_it == _rank ? *it :
typename X::field_type(0));
162 bool owned(
const ContainerIndex& i)
const
164 return _rank_partition[i] == _rank;
174 template<
typename X,
typename Y>
175 typename PromotionTraits<
176 typename X::field_type,
177 typename Y::field_type
193 template<
typename X,
typename Y,
typename Mask>
194 typename PromotionTraits<
195 typename X::field_type,
196 typename Y::field_type
200 typedef typename PromotionTraits<
201 typename X::field_type,
202 typename Y::field_type
203 >::PromotedType result_type;
207 typename Y::const_iterator y_it = y.begin();
208 typename Mask::const_iterator mask_it = mask.begin();
209 for (
typename X::const_iterator x_it = x.begin(),
212 ++x_it, ++y_it, ++mask_it)
220 template<
typename X,
typename Y,
typename Mask>
221 typename PromotionTraits<
222 typename X::field_type,
223 typename Y::field_type
227 typedef typename PromotionTraits<
228 typename X::field_type,
229 typename Y::field_type
230 >::PromotedType result_type;
234 typename Y::const_iterator y_it = y.begin();
235 typename Mask::const_iterator mask_it = mask.begin();
236 for (
typename X::const_iterator x_it = x.begin(),
239 ++x_it, ++y_it, ++mask_it)
240 r += (*mask_it == _rank ? Dune::dot(*x_it,*y_it) : result_type(0));
272 template<
typename MatrixType,
typename Comm>
273 void createIndexSetAndProjectForAMG(MatrixType& m, Comm& c);
280 bool owned_for_amg(std::size_t i)
const
290 const RankIndex _rank;
291 RankVector _rank_partition;
292 std::vector<RankIndex> _neighbor_ranks;
297 InterfaceType _interiorBorder_all_interface;
300 InterfaceType _all_all_interface;
305 template<
typename GFS>
306 template<
typename M,
typename C>
307 void ParallelHelper<GFS>::createIndexSetAndProjectForAMG(M& m, C& c)
312 const bool is_bcrs_matrix =
314 typename ISTL::tags::container<
317 ISTL::tags::bcrs_matrix
320 const bool block_type_is_field_matrix =
322 typename ISTL::tags::container<
323 typename Backend::Native<M>::block_type
325 ISTL::tags::field_matrix
329 static_assert(is_bcrs_matrix && block_type_is_field_matrix,
"matrix structure not compatible with AMG");
338 typedef typename GFS::Traits::GridViewType GV;
339 typedef typename RankVector::size_type size_type;
340 const GV& gv = _gfs.gridView();
343 const bool need_communication = _gfs.gridView().comm().size() > 1;
346 using BoolVector = Backend::Vector<GFS,bool>;
347 BoolVector sharedDOF(_gfs,
false);
349 if (need_communication)
351 SharedDOFDataHandle<GFS,BoolVector> data_handle(_gfs,sharedDOF,
false);
352 _gfs.gridView().communicate(data_handle,_all_all_interface,Dune::ForwardCommunication);
356 typedef typename C::ParallelIndexSet::GlobalIndex GlobalIndex;
357 GlobalIndex count = 0;
359 for (size_type i = 0; i < sharedDOF.N(); ++i)
360 if (owned_for_amg(i) &&
native(sharedDOF)[i][0])
363 dverb << gv.comm().rank() <<
": shared block count is " << count.touint() << std::endl;
366 std::vector<GlobalIndex> counts(_gfs.gridView().comm().size());
367 _gfs.gridView().comm().allgather(&count, 1, &(counts[0]));
370 GlobalIndex start = std::accumulate(counts.begin(),counts.begin() + _rank,GlobalIndex(0));
373 GIVector scalarIndices(_gfs, std::numeric_limits<GlobalIndex>::max());
375 for (size_type i = 0; i < sharedDOF.N(); ++i)
376 if (owned_for_amg(i) &&
native(sharedDOF)[i][0])
378 native(scalarIndices)[i][0] = start;
383 if (need_communication)
385 MinDataHandle<GFS,GIVector> data_handle(_gfs,scalarIndices);
386 _gfs.gridView().communicate(data_handle,_interiorBorder_all_interface,Dune::ForwardCommunication);
390 c.indexSet().beginResize();
391 for (size_type i=0; i<scalarIndices.N(); ++i)
393 Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
394 if(
native(scalarIndices)[i][0] != std::numeric_limits<GlobalIndex>::max())
397 if (owned_for_amg(i))
400 attr = Dune::OwnerOverlapCopyAttributeSet::owner;
404 attr = Dune::OwnerOverlapCopyAttributeSet::copy;
406 c.indexSet().add(
native(scalarIndices)[i][0],
typename C::ParallelIndexSet::LocalIndex(i,attr));
409 c.indexSet().endResize();
412 std::set<int> neighbors;
414 if (need_communication)
416 GFSNeighborDataHandle<GFS,int> data_handle(_gfs,_rank,neighbors);
417 _gfs.gridView().communicate(data_handle,_all_all_interface,Dune::ForwardCommunication);
420 c.remoteIndices().setNeighbours(neighbors);
421 c.remoteIndices().template rebuild<false>();
426 template<
int s,
bool isFakeMPIHelper>
429 typedef Dune::Amg::SequentialInformation
type;
439 typedef OwnerOverlapCopyCommunication<bigunsignedint<s>,
int>
type;
448 #if HAVE_UG && PDELAB_SEQUENTIAL_UG
450 void assertParallelUG(Dune::CollectiveCommunication<Dune::UGGrid<dim> > comm)
452 static_assert(Dune::AlwaysFalse<Dune::UGGrid<dim> >::
value,
"Using sequential UG in parallel environment");
461 #endif // DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH