1 #ifndef DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH 2 #define DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH 3 4 #include <dune/common/parametertree.hh> 5 #include <dune/common/power.hh> 6 7 #include <dune/istl/matrixmatrix.hh> 8 9 #include <dune/grid/common/datahandleif.hh> 10 #include <dune/pdelab/backend/istl/vector.hh> 11 #include <dune/pdelab/backend/istl/bcrsmatrix.hh> 12 #include <dune/pdelab/backend/istl/bcrsmatrixbackend.hh> 13 #include <dune/pdelab/backend/istl/ovlpistlsolverbackend.hh> 14 #include <dune/pdelab/gridoperator/gridoperator.hh> 15 #include <dune/pdelab/localoperator/flags.hh> 16 #include <dune/pdelab/localoperator/idefault.hh> 17 #include <dune/pdelab/localoperator/pattern.hh> 18 #include <dune/pdelab/localoperator/defaultimp.hh> 19 20 namespace Dune { 21 namespace PDELab { 22 //*********************************************************** 23 // A data handle / function to communicate matrix entries 24 // in the overlap. It turned out that it is actually not 25 // required to do that, but anyway it is there now. 26 //*********************************************************** 27 28 /** Data handle to build up local index to global id and global id to local index map for codim 0 in overlap 29 */ 30 template<class GFS> 31 class LocalGlobalMapDataHandle 32 : public Dune::CommDataHandleIF<LocalGlobalMapDataHandle<GFS>,int> 33 { 34 public: 35 // some types from the grid view 36 typedef typename GFS::Traits::GridView GV; 37 typedef typename GV::IndexSet IndexSet; 38 typedef typename IndexSet::IndexType IndexType; 39 typedef typename GV::Grid Grid; 40 typedef typename Grid::Traits::GlobalIdSet GlobalIdSet; 41 typedef typename GlobalIdSet::IdType IdType; 42 43 //! export type of data for message buffer 44 typedef int DataType; 45 46 // map types 47 typedef std::map<IndexType,IdType> LocalToGlobalMap; 48 typedef std::map<IdType,IndexType> GlobalToLocalMap; 49 50 //! returns true if data for this codim should be communicated contains(int dim,int codim) const51 bool contains (int dim, int codim) const 52 { 53 return (codim==0); 54 } 55 56 //! returns true if size per entity of given dim and codim is a constant fixedSize(int dim,int codim) const57 bool fixedSize (int dim, int codim) const 58 { 59 return true; 60 } 61 62 /*! how many objects of type DataType have to be sent for a given entity 63 64 Note: Only the sender side needs to know this size. 65 */ 66 template<class EntityType> size(EntityType & e) const67 size_t size (EntityType& e) const 68 { 69 return 1; 70 } 71 72 //! pack data from user to message buffer 73 template<class MessageBuffer, class EntityType> gather(MessageBuffer & buff,const EntityType & e) const74 void gather (MessageBuffer& buff, const EntityType& e) const 75 { 76 // fill map 77 IndexType myindex = indexSet.index(e); 78 IdType myid = globalIdSet.id(e); 79 l2g[myindex] = myid; 80 g2l[myid] = myindex; 81 //std::cout << gv.comm().rank() << ": gather local=" << myindex << " global=" << myid << std::endl; 82 83 buff.write(0); // this is a dummy, we are not interested in the data 84 } 85 86 /*! unpack data from message buffer to user 87 88 n is the number of objects sent by the sender 89 */ 90 template<class MessageBuffer, class EntityType> scatter(MessageBuffer & buff,const EntityType & e,size_t n)91 void scatter (MessageBuffer& buff, const EntityType& e, size_t n) 92 { 93 DataType x; 94 buff.read(x); // this is a dummy, we are not interested in the data 95 96 IndexType myindex = indexSet.index(e); 97 IdType myid = globalIdSet.id(e); 98 l2g[myindex] = myid; 99 g2l[myid] = myindex; 100 //std::cout << gv.comm().rank() << ": scatter local=" << myindex << " global=" << myid << std::endl; 101 } 102 103 //! constructor LocalGlobalMapDataHandle(const GFS & gfs_,LocalToGlobalMap & l2g_,GlobalToLocalMap & g2l_)104 LocalGlobalMapDataHandle (const GFS& gfs_, LocalToGlobalMap& l2g_, GlobalToLocalMap& g2l_) 105 : gfs(gfs_), gv(gfs.gridView()), indexSet(gv.indexSet()), 106 grid(gv.grid()), globalIdSet(grid.globalIdSet()), 107 l2g(l2g_), g2l(g2l_) 108 { 109 } 110 111 private: 112 const GFS& gfs; 113 GV gv; 114 const IndexSet& indexSet; 115 const Grid& grid; 116 const GlobalIdSet& globalIdSet; 117 LocalToGlobalMap& l2g; 118 GlobalToLocalMap& g2l; 119 }; 120 121 122 // A DataHandle class to exchange rows of a sparse matrix 123 template<class GFS, class M> // mapper type and vector type 124 class MatrixExchangeDataHandle 125 : public Dune::CommDataHandleIF<MatrixExchangeDataHandle<GFS,M>, 126 std::pair<typename GFS::Traits::GridView::Grid::Traits::GlobalIdSet::IdType, 127 typename M::block_type> > 128 { 129 public: 130 // some types from the grid view 131 typedef typename GFS::Traits::GridView GV; 132 typedef typename GV::IndexSet IndexSet; 133 typedef typename IndexSet::IndexType IndexType; 134 typedef typename GV::Grid Grid; 135 typedef typename Grid::Traits::GlobalIdSet GlobalIdSet; 136 typedef typename GlobalIdSet::IdType IdType; 137 138 // some types from the matrix 139 typedef typename M::block_type B; 140 141 //! export type of data for message buffer 142 typedef std::pair<IdType,B> DataType; 143 144 // map types 145 typedef std::map<IndexType,IdType> LocalToGlobalMap; 146 typedef std::map<IdType,IndexType> GlobalToLocalMap; 147 148 //! returns true if data for this codim should be communicated contains(int dim,int codim) const149 bool contains (int dim, int codim) const 150 { 151 return (codim==0); 152 } 153 154 //! returns true if size per entity of given dim and codim is a constant fixedSize(int dim,int codim) const155 bool fixedSize (int dim, int codim) const 156 { 157 return false; 158 } 159 160 /*! how many objects of type DataType have to be sent for a given entity 161 162 Note: Only the sender side needs to know this size. 163 */ 164 template<class EntityType> size(EntityType & e) const165 size_t size (EntityType& e) const 166 { 167 IndexType myindex = indexSet.index(e); 168 typename M::row_type myrow = m[myindex]; 169 typename M::row_type::iterator endit=myrow.end(); 170 size_t count=0; 171 for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it) 172 { 173 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index()); 174 if (find!=l2g.end()) 175 { 176 count++; 177 } 178 } 179 //std::cout << gv.comm().rank() << ": row=" << myindex << " size=" << count << std::endl; 180 return count; 181 } 182 183 //! pack data from user to message buffer 184 template<class MessageBuffer, class EntityType> gather(MessageBuffer & buff,const EntityType & e) const185 void gather (MessageBuffer& buff, const EntityType& e) const 186 { 187 IndexType myindex = indexSet.index(e); 188 //std::cout << gv.comm().rank() << ": begin gather local=" << myindex << std::endl; 189 typename M::row_type myrow = m[myindex]; 190 typename M::row_type::iterator endit=myrow.end(); 191 size_t count = 0; 192 for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it) 193 { 194 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index()); 195 if (find!=l2g.end()) 196 { 197 buff.write(std::make_pair(find->second,*it)); 198 //std::cout << gv.comm().rank() << ": ==> local=" << find->first << " global=" << find->second << std::endl; 199 count++; 200 } 201 } 202 //std::cout << gv.comm().rank() << ": end gather row=" << myindex << " size=" << count << std::endl; 203 } 204 205 /*! unpack data from message buffer to user 206 207 n is the number of objects sent by the sender 208 */ 209 template<class MessageBuffer, class EntityType> scatter(MessageBuffer & buff,const EntityType & e,size_t n)210 void scatter (MessageBuffer& buff, const EntityType& e, size_t n) 211 { 212 IndexType myindex = indexSet.index(e); 213 std::cout << gv.comm().rank() << ": begin scatter local=" << myindex << " size=" << n << std::endl; 214 DataType x; 215 size_t count = 0; 216 for (size_t i=0; i<n; ++i) 217 { 218 buff.read(x); 219 std::cout << gv.comm().rank() << ": --> received global=" << x.first << std::endl; 220 typename GlobalToLocalMap::const_iterator find=g2l.find(x.first); 221 if (find!=g2l.end()) 222 { 223 IndexType nbindex=find->second; 224 if (m.exists(myindex,nbindex)) 225 { 226 m[myindex][nbindex] = x.second; 227 B block(x.second); 228 block -= m2[myindex][nbindex]; 229 std::cout << gv.comm().rank() << ": compare i=" << myindex << " j=" << nbindex 230 << " norm=" << block.infinity_norm() << std::endl; 231 232 count++; 233 //std::cout << gv.comm().rank() << ": --> local=" << find->first << " global=" << find->second << std::endl; 234 } 235 } 236 } 237 //std::cout << gv.comm().rank() << ": end scatter row=" << myindex << " size=" << count << std::endl; 238 } 239 240 //! constructor MatrixExchangeDataHandle(const GFS & gfs_,M & m_,const LocalToGlobalMap & l2g_,const GlobalToLocalMap & g2l_,M & m2_)241 MatrixExchangeDataHandle (const GFS& gfs_, M& m_, const LocalToGlobalMap& l2g_, const GlobalToLocalMap& g2l_,M& m2_) 242 : gfs(gfs_), m(m_), gv(gfs.gridView()), indexSet(gv.indexSet()), 243 grid(gv.grid()), globalIdSet(grid.globalIdSet()), 244 l2g(l2g_), g2l(g2l_), m2(m2_) 245 { 246 } 247 248 private: 249 const GFS& gfs; 250 M& m; 251 GV gv; 252 const IndexSet& indexSet; 253 const Grid& grid; 254 const GlobalIdSet& globalIdSet; 255 const LocalToGlobalMap& l2g; 256 const GlobalToLocalMap& g2l; 257 M& m2; 258 }; 259 260 261 /** A function to communicate matrix entries 262 */ 263 template<class GFS, class T, class A, int n, int m> restore_overlap_entries(const GFS & gfs,Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A> & matrix,Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A> & matrix2)264 void restore_overlap_entries (const GFS& gfs, Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix, 265 Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix2) 266 { 267 typedef Dune::FieldMatrix<T,n,m> B; 268 typedef Dune::BCRSMatrix<B,A> M; // m is of type M 269 typedef typename LocalGlobalMapDataHandle<GFS>::LocalToGlobalMap LocalToGlobalMap; 270 typedef typename LocalGlobalMapDataHandle<GFS>::GlobalToLocalMap GlobalToLocalMap; 271 272 // build up two maps to associate local indices and global ids 273 LocalToGlobalMap l2g; 274 GlobalToLocalMap g2l; 275 LocalGlobalMapDataHandle<GFS> lgdh(gfs,l2g,g2l); 276 if (gfs.gridView().comm().size()>1) 277 gfs.gridView().communicate(lgdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication); 278 279 // exchange matrix entries 280 MatrixExchangeDataHandle<GFS,M> mexdh(gfs,matrix,l2g,g2l,matrix2); 281 if (gfs.gridView().comm().size()>1) 282 gfs.gridView().communicate(mexdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication); 283 } 284 285 //*********************************************************** 286 // The DG/AMG preconditioner in the overlapping case 287 //*********************************************************** 288 289 /** An ISTL preconditioner for DG based on AMG applied to CG subspace 290 291 The template parameters are: 292 DGGFS DG space 293 DGMatrix BCRSMatrix assembled with DG 294 DGPrec preconditioner to be used for DG 295 CGPrec preconditioner to be used on CG subspace 296 P BCRSMatrix for grid transfer 297 */ 298 template<class DGGFS, class DGMatrix, class DGPrec, class DGCC, 299 class CGGFS, class CGPrec, class CGCC, 300 class P, class DGHelper, class Comm> 301 class OvlpDGAMGPrec 302 : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>, 303 Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>> 304 { 305 const DGGFS& dggfs; 306 DGMatrix& dgmatrix; 307 DGPrec& dgprec; 308 const DGCC& dgcc; 309 const CGGFS& cggfs; 310 CGPrec& cgprec; 311 const CGCC& cgcc; 312 P& p; 313 const DGHelper& dghelper; 314 const Comm& comm; 315 int n1,n2; 316 317 public: 318 using V = Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>; 319 using W = Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>; 320 using X = Backend::Native<V>; 321 using Y = Backend::Native<W>; 322 using CGV = Dune::PDELab::Backend::Vector<CGGFS,typename CGPrec::domain_type::field_type>; 323 using CGW = Dune::PDELab::Backend::Vector<CGGFS,typename CGPrec::range_type::field_type>; 324 325 // define the category category() const326 SolverCategory::Category category() const override 327 { 328 return SolverCategory::overlapping; 329 } 330 331 /*! \brief Constructor. 332 333 Constructor gets all parameters to operate the prec. 334 \param A The matrix to operate on. 335 \param n The number of iterations to perform. 336 \param w The relaxation factor. 337 */ OvlpDGAMGPrec(const DGGFS & dggfs_,DGMatrix & dgmatrix_,DGPrec & dgprec_,const DGCC & dgcc_,const CGGFS & cggfs_,CGPrec & cgprec_,const CGCC & cgcc_,P & p_,const DGHelper & dghelper_,const Comm & comm_,int n1_,int n2_)338 OvlpDGAMGPrec (const DGGFS& dggfs_, DGMatrix& dgmatrix_, DGPrec& dgprec_, const DGCC& dgcc_, 339 const CGGFS& cggfs_, CGPrec& cgprec_, const CGCC& cgcc_, P& p_, 340 const DGHelper& dghelper_, const Comm& comm_, int n1_, int n2_) 341 : dggfs(dggfs_), dgmatrix(dgmatrix_), dgprec(dgprec_), dgcc(dgcc_), 342 cggfs(cggfs_), cgprec(cgprec_), cgcc(cgcc_), p(p_), dghelper(dghelper_), 343 comm(comm_), n1(n1_), n2(n2_) 344 { 345 } 346 347 /*! 348 \brief Prepare the preconditioner. 349 350 \copydoc Preconditioner::pre(X&,Y&) 351 */ pre(V & x,W & b)352 virtual void pre (V& x, W& b) override 353 { 354 using Backend::native; 355 dgprec.pre(native(x),native(b)); 356 CGW cgd(cggfs,0.0); 357 CGV cgv(cggfs,0.0); 358 cgprec.pre(native(cgv),native(cgd)); 359 } 360 361 /*! 362 \brief Apply the precondioner. 363 364 \copydoc Preconditioner::apply(X&,const Y&) 365 */ apply(V & x,const W & b)366 virtual void apply (V& x, const W& b) override 367 { 368 using Backend::native; 369 // need local copies to store defect and solution 370 W d(b); 371 Dune::PDELab::set_constrained_dofs(dgcc,0.0,d); 372 V v(x); // only to get correct size 373 374 // pre-smoothing on DG matrix 375 for (int i=0; i<n1; i++) 376 { 377 using Backend::native; 378 v = 0.0; 379 dgprec.apply(native(v),native(d)); 380 Dune::PDELab::AddDataHandle<DGGFS,V> adddh(dggfs,v); 381 if (dggfs.gridView().comm().size()>1) 382 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); 383 dgmatrix.mmv(native(v),native(d)); 384 Dune::PDELab::set_constrained_dofs(dgcc,0.0,d); 385 x += v; 386 } 387 388 // restrict defect to CG subspace 389 dghelper.maskForeignDOFs(d); // DG defect is additive for overlap 1, but in case we use more 390 CGW cgd(cggfs,0.0); 391 p.mtv(native(d),native(cgd)); 392 Dune::PDELab::AddDataHandle<CGGFS,CGW> adddh(cggfs,cgd); 393 if (cggfs.gridView().comm().size()>1) 394 cggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); // now we have consistent defect on coarse grid 395 Dune::PDELab::set_constrained_dofs(cgcc,0.0,cgd); 396 comm.project(native(cgd)); 397 CGV cgv(cggfs,0.0); 398 399 400 // call preconditioner 401 cgprec.apply(native(cgv),native(cgd)); 402 403 // prolongate correction 404 p.mv(native(cgv),native(v)); 405 dgmatrix.mmv(native(v),native(d)); 406 Dune::PDELab::set_constrained_dofs(dgcc,0.0,d); 407 x += v; 408 409 // post-smoothing on DG matrix 410 for (int i=0; i<n2; i++) 411 { 412 v = 0.0; 413 dgprec.apply(native(v),native(d)); 414 Dune::PDELab::AddDataHandle<DGGFS,V> adddh(dggfs,v); 415 if (dggfs.gridView().comm().size()>1) 416 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); 417 dgmatrix.mmv(native(v),native(d)); 418 Dune::PDELab::set_constrained_dofs(dgcc,0.0,d); 419 x += v; 420 } 421 } 422 423 /*! 424 \brief Clean up. 425 426 \copydoc Preconditioner::post(X&) 427 */ post(V & x)428 virtual void post (V& x) override 429 { 430 dgprec.post(Backend::native(x)); 431 CGV cgv(cggfs,0.0); 432 cgprec.post(Backend::native(cgv)); 433 } 434 }; 435 436 //*********************************************************** 437 // The DG/AMG linear solver backend in the overlapping case 438 //*********************************************************** 439 440 /** Overlapping solver backend for using AMG for DG in PDELab 441 442 The template parameters are: 443 DGGO GridOperator for DG discretization, allows access to matrix, vector and grid function space 444 DGCC constraints container for DG problem 445 CGGFS grid function space for CG subspace 446 CGCC constraints container for CG problem 447 TransferLOP local operator to assemble prolongation from CGGFS to DGGFS 448 DGPrec preconditioner for DG problem 449 Solver solver to be used on the complete problem 450 int s size of global index to be used in AMG 451 452 Note: The subspace matrix is calculated by a triple matrix product with the 453 fine space DG matrix passed into the apply method. This only works if this 454 DG matrix does not contain P0ParallelConstraints. In order to achieve this 455 behaviour you should use this solver backend with a StationaryLinearProblem 456 or Newton solver that was created using a grid operator with empty 457 constranints. See the overlapping test case in 458 459 dune-pdelab/dune/pdelab/test-dg-amg.cc 460 461 for an example. 462 */ 463 template<class DGGO, class DGCC, class CGGFS, class CGCC, class TransferLOP, 464 template<class,class,class,int> class DGPrec, template<class> class Solver, int s=96> 465 class ISTLBackend_OVLP_AMG_4_DG : 466 public Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>, 467 public Dune::PDELab::LinearResultStorage 468 { 469 public: 470 // DG grid function space 471 typedef typename DGGO::Traits::TrialGridFunctionSpace GFS; 472 473 // vectors and matrices on DG level 474 typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix 475 typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector 476 typedef Backend::Native<M> Matrix; // istl DG matrix 477 typedef Backend::Native<V> Vector; // istl DG vector 478 typedef typename Vector::field_type field_type; 479 480 // vectors and matrices on CG level 481 using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector 482 typedef Backend::Native<CGV> CGVector; // istl CG vector 483 484 // prolongation matrix 485 typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE; 486 typedef Dune::PDELab::EmptyTransformation CC; 487 typedef TransferLOP CGTODGLOP; // local operator 488 typedef Dune::PDELab::GridOperator<CGGFS,GFS,CGTODGLOP,MBE,field_type,field_type,field_type,CC,CC> PGO; 489 typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix 490 typedef Backend::Native<PMatrix> P; // ISTL prolongation matrix 491 492 // CG subspace matrix 493 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG; 494 typedef typename Dune::MatMultMatResult<PTADG,P>::type CGMatrix; // istl coarse space matrix 495 496 // AMG in CG-subspace 497 typedef typename Dune::PDELab::ISTL::CommSelector<s,Dune::MPIHelper::isFake>::type Comm; 498 typedef Dune::OverlappingSchwarzOperator<CGMatrix,CGVector,CGVector,Comm> ParCGOperator; 499 typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother; 500 typedef Dune::BlockPreconditioner<CGVector,CGVector,Comm,Smoother> ParSmoother; 501 typedef Dune::Amg::AMG<ParCGOperator,CGVector,ParSmoother,Comm> AMG; 502 typedef Dune::Amg::Parameters Parameters; 503 504 private: 505 506 const GFS& gfs; 507 DGGO& dggo; 508 const DGCC& dgcc; 509 CGGFS& cggfs; 510 const CGCC& cgcc; 511 std::shared_ptr<AMG> amg; 512 Parameters amg_parameters; 513 unsigned maxiter; 514 int verbose; 515 bool reuse; 516 bool firstapply; 517 bool usesuperlu; 518 std::size_t low_order_space_entries_per_row; 519 520 // Parameters for DG smoother 521 double smoother_relaxation = 0.92; // Relaxation parameter of DG smoother 522 int n1=3; // Number of DG pre-smoothing steps 523 int n2=3; // Number of DG post-smoothing steps 524 525 CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix 526 PGO pgo; // grid operator to assemble prolongation matrix 527 PMatrix pmatrix; // wrapped prolongation matrix 528 529 /** an empty local operator to assemble processor boundary constraints 530 */ 531 class EmptyLop : public Dune::PDELab::NumericalJacobianApplyVolume<EmptyLop>, 532 public Dune::PDELab::FullVolumePattern, 533 public Dune::PDELab::LocalOperatorDefaultFlags, 534 public Dune::PDELab::InstationaryLocalOperatorDefaultMethods<double> 535 { 536 }; 537 538 // grid operator with empty local operator for constraints assembly in parallel 539 typedef Dune::PDELab::GridOperator<CGGFS,CGGFS,EmptyLop,MBE,field_type,field_type,field_type,CGCC,CGCC> CGGO; 540 // CG-subspace matrix 541 typedef typename CGGO::Jacobian CGM; 542 CGM acg; 543 544 public: 545 546 // access to prolongation matrix prolongation_matrix()547 PMatrix& prolongation_matrix () 548 { 549 return pmatrix; 550 } 551 552 /*! \brief set AMG parameters 553 554 \param[in] amg_parameters_ a parameter object of Type Dune::Amg::Parameters 555 */ setParameters(const Parameters & amg_parameters_)556 void setParameters(const Parameters& amg_parameters_) 557 { 558 amg_parameters = amg_parameters_; 559 } 560 561 /** 562 * @brief Get the parameters describing the behaviuour of AMG. 563 * 564 * The returned object can be adjusted to ones needs and then can be 565 * reset using setParameters. 566 * @return The object holding the parameters of AMG. 567 */ parameters() const568 const Parameters& parameters() const 569 { 570 return amg_parameters; 571 } 572 573 //! Set whether the AMG should be reused again during call to apply(). setReuse(bool reuse_)574 void setReuse(bool reuse_) 575 { 576 reuse = reuse_; 577 } 578 579 //! Return whether the AMG is reused during call to apply() getReuse() const580 bool getReuse() const 581 { 582 return reuse; 583 } 584 585 /** make backend object 586 */ ISTLBackend_OVLP_AMG_4_DG(DGGO & dggo_,const DGCC & dgcc_,CGGFS & cggfs_,const CGCC & cgcc_,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)587 ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_, 588 unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, 589 bool usesuperlu_=true) 590 : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace()) 591 , gfs(dggo_.trialGridFunctionSpace()) 592 , dggo(dggo_) 593 , dgcc(dgcc_) 594 , cggfs(cggfs_) 595 , cgcc(cgcc_) 596 , amg_parameters(15,2000) 597 , maxiter(maxiter_) 598 , verbose(verbose_) 599 , reuse(reuse_) 600 , firstapply(true) 601 , usesuperlu(usesuperlu_) 602 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power) 603 , cgtodglop() 604 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row)) 605 , pmatrix(pgo) 606 , acg(Backend::attached_container()) 607 { 608 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension); 609 amg_parameters.setDebugLevel(verbose_); 610 #if !HAVE_SUPERLU 611 if (usesuperlu == true) 612 { 613 if (gfs.gridView().comm().rank()==0) 614 std::cout << "WARNING: You are using AMG without SuperLU!" 615 << " Please consider installing SuperLU," 616 << " or set the usesuperlu flag to false" 617 << " to suppress this warning." << std::endl; 618 } 619 #endif 620 621 // assemble prolongation matrix; this will not change from one apply to the next 622 pmatrix = 0.0; 623 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl; 624 CGV cgx(cggfs,0.0); // need vector to call jacobian 625 pgo.jacobian(cgx,pmatrix); 626 } 627 628 629 /** make backend object 630 */ ISTLBackend_OVLP_AMG_4_DG(DGGO & dggo_,const DGCC & dgcc_,CGGFS & cggfs_,const CGCC & cgcc_,const ParameterTree & params)631 ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_, 632 const ParameterTree& params) 633 : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace()) 634 , gfs(dggo_.trialGridFunctionSpace()) 635 , dggo(dggo_) 636 , dgcc(dgcc_) 637 , cggfs(cggfs_) 638 , cgcc(cgcc_) 639 , maxiter(params.get<int>("max_iterations",5000)) 640 , amg_parameters(15,2000) 641 , verbose(params.get<int>("verbose",1)) 642 , reuse(params.get<bool>("reuse",false)) 643 , firstapply(true) 644 , usesuperlu(params.get<bool>("use_superlu",true)) 645 , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power)) 646 , cgtodglop() 647 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row)) 648 , pmatrix(pgo) 649 , acg(Backend::attached_container()) 650 { 651 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension); 652 amg_parameters.setDebugLevel(params.get<int>("verbose",1)); 653 #if !HAVE_SUPERLU 654 if (usesuperlu == true) 655 { 656 if (gfs.gridView().comm().rank()==0) 657 std::cout << "WARNING: You are using AMG without SuperLU!" 658 << " Please consider installing SuperLU," 659 << " or set the usesuperlu flag to false" 660 << " to suppress this warning." << std::endl; 661 } 662 #endif 663 664 // assemble prolongation matrix; this will not change from one apply to the next 665 pmatrix = 0.0; 666 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl; 667 CGV cgx(cggfs,0.0); // need vector to call jacobian 668 pgo.jacobian(cgx,pmatrix); 669 } 670 671 672 //! set number of presmoothing steps on the DG level setDGSmootherRelaxation(double relaxation_)673 void setDGSmootherRelaxation(double relaxation_) 674 { 675 smoother_relaxation = relaxation_; 676 } 677 678 //! set number of presmoothing steps on the DG level setNoDGPreSmoothSteps(int n1_)679 void setNoDGPreSmoothSteps(int n1_) 680 { 681 n1 = n1_; 682 } 683 684 //! set number of postsmoothing steps on the DG level setNoDGPostSmoothSteps(int n2_)685 void setNoDGPostSmoothSteps(int n2_) 686 { 687 n2 = n2_; 688 } 689 690 691 /*! \brief solve the given linear system 692 693 \param[in] A the given matrix 694 \param[out] z the solution vector to be computed 695 \param[in] r right hand side 696 \param[in] reduction to be achieved 697 */ apply(M & A,V & z,V & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)698 void apply (M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 699 { 700 using Backend::native; 701 // make operator and scalar product for overlapping solver 702 typedef Dune::PDELab::OverlappingOperator<DGCC,M,V,V> POP; 703 POP pop(dgcc,A); 704 typedef Dune::PDELab::OVLPScalarProduct<GFS,V> PSP; 705 PSP psp(*this); 706 707 // compute CG matrix 708 // make grid operator with empty local operator => matrix data type and constraints assembly 709 EmptyLop emptylop; 710 CGGO cggo(cggfs,cgcc,cggfs,cgcc,emptylop,MBE(low_order_space_entries_per_row)); 711 712 // do triple matrix product ACG = P^T ADG P; this is purely local 713 Dune::Timer watch; 714 watch.reset(); 715 // only do triple matrix product if the matrix changes 716 double triple_product_time = 0.0; 717 // no need to set acg here back to zero, this is done in matMultmat 718 if(reuse == false || firstapply == true) { 719 PTADG ptadg; 720 Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A)); // 1a 721 //Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A2)); // 1b 722 Dune::matMultMat(native(acg),ptadg,native(pmatrix)); 723 triple_product_time = watch.elapsed(); 724 if (verbose>0 && gfs.gridView().comm().rank()==0) 725 std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl; 726 //Dune::printmatrix(std::cout,native(acg),"triple product matrix","row",10,2); 727 CGV cgx(cggfs,0.0); // need vector to call jacobian 728 cggo.jacobian(cgx,acg); // insert trivial rows at processor boundaries 729 //std::cout << "CG constraints: " << cgcc.size() << " out of " << cggfs.globalSize() << std::endl; 730 } 731 else if(verbose>0 && gfs.gridView().comm().rank()==0) 732 std::cout << "=== reuse CG matrix, SKIPPING triple matrix product " << std::endl; 733 734 // NOW we need to insert the processor boundary conditions in DG matrix 735 typedef Dune::PDELab::GridOperator<GFS,GFS,EmptyLop,MBE,field_type,field_type,field_type,DGCC,DGCC> DGGOEmpty; 736 DGGOEmpty dggoempty(gfs,dgcc,gfs,dgcc,emptylop,MBE(1 << GFS::Traits::GridView::dimension)); 737 dggoempty.jacobian(z,A); 738 739 // and in the residual 740 Dune::PDELab::set_constrained_dofs(dgcc,0.0,r); 741 742 // now set up parallel AMG solver for the CG subspace 743 Comm oocc(gfs.gridView().comm()); 744 typedef Dune::PDELab::ISTL::ParallelHelper<CGGFS> CGHELPER; 745 CGHELPER cghelper(cggfs,2); 746 cghelper.createIndexSetAndProjectForAMG(acg,oocc); 747 ParCGOperator paroop(native(acg),oocc); 748 Dune::OverlappingSchwarzScalarProduct<CGVector,Comm> sp(oocc); 749 750 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs; 751 SmootherArgs smootherArgs; 752 smootherArgs.iterations = 1; 753 smootherArgs.relaxationFactor = 1.0; 754 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion; 755 Criterion criterion(amg_parameters); 756 watch.reset(); 757 758 // only construct a new AMG for the CG-subspace if the matrix changes 759 double amg_setup_time = 0.0; 760 if(reuse == false || firstapply == true) { 761 amg.reset(new AMG(paroop,criterion,smootherArgs,oocc)); 762 firstapply = false; 763 amg_setup_time = watch.elapsed(); 764 if (verbose>0 && gfs.gridView().comm().rank()==0) 765 std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl; 766 } 767 else if (verbose>0 && gfs.gridView().comm().rank()==0) 768 std::cout << "=== reuse CG matrix, SKIPPING AMG setup " << std::endl; 769 770 // set up hybrid DG/CG preconditioner 771 typedef DGPrec<Matrix,Vector,Vector,1> DGPrecType; 772 DGPrecType dgprec(native(A),1,smoother_relaxation); 773 typedef Dune::PDELab::ISTL::ParallelHelper<GFS> DGHELPER; 774 typedef OvlpDGAMGPrec<GFS,Matrix,DGPrecType,DGCC,CGGFS,AMG,CGCC,P,DGHELPER,Comm> HybridPrec; 775 HybridPrec hybridprec(gfs,native(A),dgprec,dgcc,cggfs,*amg,cgcc,native(pmatrix), 776 this->parallelHelper(),oocc,n1,n2); 777 778 // /********/ 779 // /* Test */ 780 // /********/ 781 // Solver<CGVector> testsolver(paroop,sp,amg,1e-8,100,2); 782 // CGV cgxx(cggfs,0.0); 783 // CGV cgdd(cggfs,1.0); 784 // Dune::InverseOperatorResult statstat; 785 // testsolver.apply(native(cgxx),native(cgdd),statstat); 786 // /********/ 787 788 // set up solver 789 int verb=verbose; 790 if (gfs.gridView().comm().rank()>0) verb=0; 791 Solver<V> solver(pop,psp,hybridprec,reduction,maxiter,verb); 792 793 // solve 794 Dune::InverseOperatorResult stat; 795 watch.reset(); 796 solver.apply(z,r,stat); 797 double amg_solve_time = watch.elapsed(); 798 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl; 799 res.converged = stat.converged; 800 res.iterations = stat.iterations; 801 res.elapsed = amg_solve_time+amg_setup_time+triple_product_time; 802 res.reduction = stat.reduction; 803 res.conv_rate = stat.conv_rate; 804 } 805 806 }; 807 } 808 } 809 #endif // DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH 810