dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set ts=4 sw=2 et sts=2: 00003 #ifndef DUNE_NOVLPSCHWARZ_HH 00004 #define DUNE_NOVLPSCHWARZ_HH 00005 00006 #include <iostream> // for input/output to shell 00007 #include <fstream> // for input/output to files 00008 #include <vector> // STL vector class 00009 #include <sstream> 00010 00011 #include <cmath> // Yes, we do some math here 00012 00013 #include <dune/common/timer.hh> 00014 00015 #include"io.hh" 00016 #include"bvector.hh" 00017 #include"vbvector.hh" 00018 #include"bcrsmatrix.hh" 00019 #include"io.hh" 00020 #include"gsetc.hh" 00021 #include"ilu.hh" 00022 #include"operators.hh" 00023 #include"solvers.hh" 00024 #include"preconditioners.hh" 00025 #include"scalarproducts.hh" 00026 #include"owneroverlapcopy.hh" 00027 00028 namespace Dune { 00029 00058 template<class M, class X, class Y, class C> 00059 class NonoverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y> 00060 { 00061 public: 00063 typedef M matrix_type; 00065 typedef X domain_type; 00067 typedef Y range_type; 00069 typedef typename X::field_type field_type; 00071 typedef C communication_type; 00072 00073 typedef typename C::PIS PIS; 00074 typedef typename C::RI RI; 00075 typedef typename RI::RemoteIndexList RIL; 00076 typedef typename RI::const_iterator RIIterator; 00077 typedef typename RIL::const_iterator RILIterator; 00078 typedef typename M::ConstColIterator ColIterator; 00079 typedef typename M::ConstRowIterator RowIterator; 00080 typedef std::multimap<int,int> MM; 00081 typedef std::multimap<int,std::pair<int,RILIterator> > RIMap; 00082 typedef typename RIMap::iterator RIMapit; 00083 00084 enum { 00086 category=SolverCategory::nonoverlapping 00087 }; 00088 00096 NonoverlappingSchwarzOperator (const matrix_type& A, const communication_type& com) 00097 : _A_(A), communication(com), buildcomm(true) 00098 {} 00099 00101 virtual void apply (const X& x, Y& y) const 00102 { 00103 y = 0; 00104 novlp_op_apply(x,y,1); 00105 communication.addOwnerCopyToOwnerCopy(y,y); 00106 } 00107 00109 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const 00110 { 00111 // only apply communication to alpha*A*x to make it consistent, 00112 // y already has to be consistent. 00113 Y y1(y); 00114 y = 0; 00115 novlp_op_apply(x,y,alpha); 00116 communication.addOwnerCopyToOwnerCopy(y,y); 00117 y += y1; 00118 } 00119 00121 virtual const matrix_type& getmat () const 00122 { 00123 return _A_; 00124 } 00125 00126 void novlp_op_apply (const X& x, Y& y, field_type alpha) const 00127 { 00128 //get index sets 00129 const PIS& pis=communication.indexSet(); 00130 const RI& ri = communication.remoteIndices(); 00131 00132 // at the beginning make a multimap "bordercontribution". 00133 // process has i and j as border dofs but is not the owner 00134 // => only contribute to Ax if i,j is in bordercontribution 00135 if (buildcomm == true){ 00136 00137 // set up mask vector 00138 if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())){ 00139 mask.resize(x.size()); 00140 for (typename std::vector<double>::size_type i=0; i<mask.size(); i++) 00141 mask[i] = 1; 00142 for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i) 00143 if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy) 00144 mask[i->local().local()] = 0; 00145 else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap) 00146 mask[i->local().local()] = 2; 00147 } 00148 00149 for (MM::iterator iter = bordercontribution.begin(); 00150 iter != bordercontribution.end(); ++iter) 00151 bordercontribution.erase(iter); 00152 std::map<int,int> owner; //key: local index i, value: process, that owns i 00153 RIMap rimap; 00154 00155 // for each local index make multimap rimap: 00156 // key: local index i, data: pair of process that knows i and pointer to RI entry 00157 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i) 00158 if (mask[i.index()] == 0) 00159 for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote){ 00160 RIL& ril = *(remote->second.first); 00161 for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex) 00162 if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap) 00163 if (rindex->localIndexPair().local().local() == i.index()){ 00164 rimap.insert 00165 (std::make_pair(i.index(), 00166 std::pair<int,RILIterator>(remote->first, rindex))); 00167 if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner) 00168 owner.insert(std::make_pair(i.index(),remote->first)); 00169 } 00170 } 00171 00172 int iowner = 0; 00173 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i){ 00174 if (mask[i.index()] == 0){ 00175 std::map<int,int>::iterator it = owner.find(i.index()); 00176 iowner = it->second; 00177 std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index()); 00178 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j){ 00179 if (mask[j.index()] == 0){ 00180 bool flag = true; 00181 for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi){ 00182 std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index()); 00183 for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj) 00184 if (foundj->second.first == foundi->second.first) 00185 if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner 00186 || foundj->second.first == iowner 00187 || foundj->second.first < communication.communicator().rank()){ 00188 flag = false; 00189 continue; 00190 } 00191 if (flag == false) 00192 continue; 00193 } 00194 // don“t contribute to Ax if 00195 // 1. the owner of j has i as interior/border dof 00196 // 2. iowner has j as interior/border dof 00197 // 3. there is another process with smaller rank that has i and j 00198 // as interor/border dofs 00199 // if the owner of j does not have i as interior/border dof, 00200 // it will not be taken into account 00201 if (flag==true) 00202 bordercontribution.insert(std::pair<int,int>(i.index(),j.index())); 00203 } 00204 } 00205 } 00206 } 00207 buildcomm = false; 00208 } 00209 00210 //compute alpha*A*x nonoverlapping case 00211 for (RowIterator i = _A_.begin(); i != _A_.end(); ++i){ 00212 if (mask[i.index()] == 0){ 00213 //dof doesn't belong to process but is border (not ghost) 00214 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j){ 00215 if (mask[j.index()] == 1) //j is owner => then sum entries 00216 (*j).usmv(alpha,x[j.index()],y[i.index()]); 00217 else if (mask[j.index()] == 0){ 00218 std::pair<MM::iterator, MM::iterator> itp = 00219 bordercontribution.equal_range(i.index()); 00220 for (MM::iterator it = itp.first; it != itp.second; ++it) 00221 if ((*it).second == (int)j.index()) 00222 (*j).usmv(alpha,x[j.index()],y[i.index()]); 00223 } 00224 } 00225 } 00226 else if (mask[i.index()] == 1){ 00227 for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j) 00228 if (mask[j.index()] != 2) 00229 (*j).usmv(alpha,x[j.index()],y[i.index()]); 00230 } 00231 } 00232 } 00233 00234 private: 00235 const matrix_type& _A_; 00236 const communication_type& communication; 00237 mutable bool buildcomm; 00238 mutable std::vector<double> mask; 00239 mutable std::multimap<int,int> bordercontribution; 00240 }; 00241 00253 template<class X, class C> 00254 class NonoverlappingSchwarzScalarProduct : public ScalarProduct<X> 00255 { 00256 public: 00258 typedef X domain_type; 00260 typedef typename X::field_type field_type; 00262 typedef C communication_type; 00263 00265 enum {category=SolverCategory::nonoverlapping}; 00266 00271 NonoverlappingSchwarzScalarProduct (const communication_type& com) 00272 : communication(com) 00273 {} 00274 00279 virtual field_type dot (const X& x, const X& y) 00280 { 00281 field_type result; 00282 communication.dot(x,y,result); 00283 return result; 00284 } 00285 00289 virtual double norm (const X& x) 00290 { 00291 return communication.norm(x); 00292 } 00293 00296 void make_consistent (X& x) const 00297 { 00298 communication.copyOwnerToAll(x,x); 00299 } 00300 00301 private: 00302 const communication_type& communication; 00303 }; 00304 00305 template<class X, class C> 00306 struct ScalarProductChooser<X,C,SolverCategory::nonoverlapping> 00307 { 00309 typedef NonoverlappingSchwarzScalarProduct<X,C> ScalarProduct; 00311 typedef C communication_type; 00312 00313 enum{ 00315 solverCategory=SolverCategory::nonoverlapping 00316 }; 00317 00318 static ScalarProduct* construct(const communication_type& comm) 00319 { 00320 return new ScalarProduct(comm); 00321 } 00322 }; 00323 00324 namespace Amg 00325 { 00326 template<class T> class ConstructionTraits; 00327 } 00328 00338 template<class C, class P> 00339 class NonoverlappingBlockPreconditioner 00340 : public Dune::Preconditioner<typename P::domain_type,typename P::range_type> { 00341 friend class Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >; 00342 public: 00344 typedef typename P::domain_type domain_type; 00346 typedef typename P::range_type range_type; 00348 typedef C communication_type; 00349 00350 // define the category 00351 enum { 00353 category=SolverCategory::nonoverlapping}; 00354 00362 NonoverlappingBlockPreconditioner (P& prec, const communication_type& c) 00363 : preconditioner(prec), communication(c) 00364 {} 00365 00371 virtual void pre (domain_type& x, range_type& b) 00372 { 00373 preconditioner.pre(x,b); 00374 } 00375 00381 virtual void apply (domain_type& v, const range_type& d) 00382 { 00383 // block preconditioner equivalent to WrappedPreconditioner from 00384 // pdelab/backend/ovlpistsolverbackend.hh, 00385 // but not to BlockPreconditioner from schwarz.hh 00386 preconditioner.apply(v,d); 00387 communication.addOwnerCopyToOwnerCopy(v,v); 00388 } 00389 00395 virtual void post (domain_type& x) 00396 { 00397 preconditioner.post(x); 00398 } 00399 00400 private: 00402 P& preconditioner; 00403 00405 const communication_type& communication; 00406 }; 00407 00410 } // end namespace 00411 00412 #endif