3#ifndef DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
4#define DUNE_PDELAB_BACKEND_ISTL_PARALLELHELPER_HH
12#if HAVE_DUNE_UGGRID && PDELAB_SEQUENTIAL_UG
18#include <dune/istl/solvercategory.hh>
24#include <dune/istl/paamg/pinfo.hh>
28#include <dune/pdelab/constraints/common/constraints.hh>
29#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
30#include <dune/pdelab/backend/interface.hh>
31#include <dune/pdelab/backend/istl/vector.hh>
32#include <dune/pdelab/backend/istl/utility.hh>
33#include <dune/pdelab/gridfunctionspace/tags.hh>
49 template<
typename GFS>
54 typedef int RankIndex;
57 using RankVector = Dune::PDELab::Backend::Vector<GFS,RankIndex>;
59 using GhostVector = Dune::PDELab::Backend::Vector<GFS,bool>;
62 typedef typename GFS::Ordering::Traits::ContainerIndex ContainerIndex;
66 ParallelHelper (
const GFS& gfs,
int verbose = 1)
68 , _rank(gfs.gridView().comm().rank())
69 , _rank_partition(gfs,_rank)
91 if (_gfs.gridView().comm().size()>1)
96 GhostDataHandle<GFS,GhostVector>
97 gdh(_gfs,_ghosts,
false);
103 DisjointPartitioningDataHandle<GFS,RankVector> pdh(_gfs,_rank_partition);
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);
119 const std::vector<RankIndex>& getNeighborRanks()
const
121 return _neighbor_ranks;
126 void maskForeignDOFs(X& x)
const
128 using Backend::native;
130 maskForeignDOFs(ISTL::container_tag(native(x)),native(x),native(_rank_partition));
136 template<
typename X,
typename Mask>
137 void maskForeignDOFs(ISTL::tags::block_vector, X& x,
const Mask&
mask)
const
139 typename Mask::const_iterator mask_it =
mask.begin();
140 for (
typename X::iterator it = x.begin(),
144 maskForeignDOFs(ISTL::container_tag(*it),*it,*mask_it);
148 template<
typename X,
typename Mask>
149 void maskForeignDOFs(ISTL::tags::field_vector, X& x,
const Mask&
mask)
const
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;
168 bool isGhost(
const ContainerIndex& i)
const
174 template<
typename X,
typename Y>
175 typename PromotionTraits<
176 typename X::field_type,
177 typename Y::field_type
179 disjointDot(
const X& x,
const Y& y)
const
181 using Backend::native;
182 return disjointDot(ISTL::container_tag(native(x)),
185 native(_rank_partition)
193 template<
typename X,
typename Y,
typename Mask>
194 typename PromotionTraits<
195 typename X::field_type,
196 typename Y::field_type
198 disjointDot(ISTL::tags::block_vector,
const X& x,
const Y& y,
const Mask&
mask)
const
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)
213 r += disjointDot(ISTL::container_tag(*x_it),*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
225 disjointDot(ISTL::tags::field_vector,
const X& x,
const Y& y,
const Mask&
mask)
const
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));
248 RankIndex rank()
const
272 template<
typename MatrixType,
typename Comm>
273 void createIndexSetAndProjectForAMG(MatrixType& m, Comm& c);
280 bool owned_for_amg(std::size_t i)
const
282 return Backend::native(_rank_partition)[i][0] == _rank;
290 const RankIndex _rank;
291 RankVector _rank_partition;
292 std::vector<RankIndex> _neighbor_ranks;
305 template<
typename GFS>
306 template<
typename M,
typename C>
307 void ParallelHelper<GFS>::createIndexSetAndProjectForAMG(M& m, C& c)
310 using Backend::native;
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);
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));
372 using GIVector = Dune::PDELab::Backend::Vector<GFS,GlobalIndex>;
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);
390 c.indexSet().beginResize();
391 for (size_type i=0; i<scalarIndices.N(); ++i)
393 Dune::OwnerOverlapCopyAttributeSet::AttributeSet attr;
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);
420 c.remoteIndices().setNeighbours(neighbors);
421 c.remoteIndices().template rebuild<false>();
426 template<
int s,
bool isFakeMPIHelper>
429 typedef Dune::Amg::SequentialInformation type;
437 struct CommSelector<s,false>
439 typedef OwnerOverlapCopyCommunication<bigunsignedint<s>,
int> type;
445 void assertParallelUG(T comm)
448#if HAVE_DUNE_UGGRID && PDELAB_SEQUENTIAL_UG
Collective communication interface and sequential default implementation.
Definition: communication.hh:100
Front-end for the grid manager of the finite element toolbox UG3.
Definition: uggrid.hh:208
Traits for type conversions and type information.
Some generic functions for pretty printing vectors and matrices.
Implementations of the inverse operator interface.
auto dot(const A &a, const B &b) -> typename std::enable_if< IsNumber< A >::value &&!IsVector< A >::value &&!std::is_same< typename FieldTraits< A >::field_type, typename FieldTraits< A >::real_type > ::value, decltype(conj(a) *b)>::type
computes the dot product for fundamental data types according to Petsc's VectDot function: dot(a,...
Definition: dotproduct.hh:42
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition: defaults.hh:153
Rebind< bool, V > Mask
Mask type type of some SIMD type.
Definition: interface.hh:289
DVerbType dverb(std::cout)
Singleton of verbose debug stream.
Definition: stdstreams.hh:117
Helpers for dealing with MPI.
constexpr InteriorBorder interiorBorder
PartitionSet for the interior and border partitions.
Definition: partitionset.hh:286
Dune namespace.
Definition: alignedallocator.hh:13
Define general, extensible interface for operators. The available implementation wraps a matrix.
Classes providing communication interfaces for overlapping Schwarz methods.
Define general preconditioner interface.
Define base class for scalar product and norm.
Standard Dune debug streams.
template which always yields a false value
Definition: typetraits.hh:124
Classes for using SuperLU with ISTL matrices.