1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=8 sw=2 sts=2: 3 #ifndef DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH 4 #define DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH 5 6 #include <cstddef> 7 8 #include <dune/common/deprecated.hh> 9 #include <dune/common/parallel/mpihelper.hh> 10 11 #include <dune/grid/common/gridenums.hh> 12 13 #include <dune/istl/io.hh> 14 #include <dune/istl/operators.hh> 15 #include <dune/istl/owneroverlapcopy.hh> 16 #include <dune/istl/paamg/amg.hh> 17 #include <dune/istl/paamg/pinfo.hh> 18 #include <dune/istl/preconditioners.hh> 19 #include <dune/istl/scalarproducts.hh> 20 #include <dune/istl/solvercategory.hh> 21 #include <dune/istl/solvers.hh> 22 #include <dune/istl/superlu.hh> 23 24 #include <dune/pdelab/constraints/common/constraints.hh> 25 #include <dune/pdelab/gridfunctionspace/genericdatahandle.hh> 26 #include <dune/pdelab/backend/istl/vector.hh> 27 #include <dune/pdelab/backend/istl/bcrsmatrix.hh> 28 #include <dune/pdelab/backend/istl/blockmatrixdiagonal.hh> 29 #include <dune/pdelab/backend/istl/parallelhelper.hh> 30 #include <dune/pdelab/backend/istl/seqistlsolverbackend.hh> 31 32 namespace Dune { 33 namespace PDELab { 34 35 //! \addtogroup Backend 36 //! \ingroup PDELab 37 //! \{ 38 39 //======================================================== 40 // Generic support for nonoverlapping grids 41 //======================================================== 42 43 //! Operator for the non-overlapping parallel case 44 /** 45 * Calculate \f$y:=Ax\f$. 46 * 47 * \tparam GFS The GridFunctionSpace the vectors apply to. 48 * \tparam M Type of the matrix. Should be one of the ISTL matrix types. 49 * \tparam X Type of the vectors the matrix is applied to. 50 * \tparam Y Type of the result vectors. 51 */ 52 template<typename GFS, typename M, typename X, typename Y> 53 class NonoverlappingOperator 54 : public Dune::AssembledLinearOperator<M,X,Y> 55 { 56 public: 57 //! export type of matrix 58 using matrix_type = Backend::Native<M>; 59 //! export type of vectors the matrix is applied to 60 using domain_type = Backend::Native<X>; 61 //! export type of result vectors 62 using range_type = Backend::Native<Y>; 63 //! export type of the entries for x 64 typedef typename X::field_type field_type; 65 66 //! Construct a non-overlapping operator 67 /** 68 * \param gfs_ GridFunctionsSpace for the vectors. 69 * \param A Matrix for this operator. This should be the locally 70 * assembled matrix. 71 * 72 * \note The constructed object stores references to all the objects 73 * given as parameters here. They should be valid for as long as 74 * the constructed object is used. They are not needed to 75 * destruct the constructed object. 76 */ NonoverlappingOperator(const GFS & gfs_,const M & A)77 NonoverlappingOperator (const GFS& gfs_, const M& A) 78 : gfs(gfs_), _A_(A) 79 { } 80 81 //! apply operator 82 /** 83 * Compute \f$y:=A(x)\f$ on this process, then make y consistent (sum up 84 * corresponding entries of y on the different processes and store the 85 * result back in y on each process). 86 */ apply(const X & x,Y & y) const87 virtual void apply (const X& x, Y& y) const override 88 { 89 using Backend::native; 90 // apply local operator; now we have sum y_p = sequential y 91 native(_A_).mv(native(x),native(y)); 92 93 // accumulate y on border 94 Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y); 95 if (gfs.gridView().comm().size()>1) 96 gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication); 97 } 98 99 //! apply operator to x, scale and add: \f$ y = y + \alpha A(x) \f$ 100 /** 101 * Compute \f$y:=\alpha A(x)\f$ on this process, then make y consistent 102 * (sum up corresponding entries of y on the different processes and 103 * store the result back in y on each process). 104 */ applyscaleadd(field_type alpha,const X & x,Y & y) const105 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const override 106 { 107 using Backend::native; 108 // apply local operator; now we have sum y_p = sequential y 109 native(_A_).usmv(alpha,native(x),native(y)); 110 111 // accumulate y on border 112 Dune::PDELab::AddDataHandle<GFS,Y> adddh(gfs,y); 113 if (gfs.gridView().comm().size()>1) 114 gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication); 115 } 116 category() const117 SolverCategory::Category category() const override 118 { 119 return SolverCategory::nonoverlapping; 120 } 121 122 //! extract the matrix getmat() const123 virtual const M& getmat () const override 124 { 125 return _A_; 126 } 127 128 private: 129 const GFS& gfs; 130 const M& _A_; 131 }; 132 133 // parallel scalar product assuming no overlap 134 template<class GFS, class X> 135 class NonoverlappingScalarProduct : public Dune::ScalarProduct<X> 136 { 137 public: 138 //! export types 139 typedef X domain_type; 140 typedef typename X::ElementType field_type; 141 category() const142 SolverCategory::Category category() const override 143 { 144 return SolverCategory::nonoverlapping; 145 } 146 147 /*! \brief Constructor needs to know the grid function space 148 */ NonoverlappingScalarProduct(const GFS & gfs_,const ISTL::ParallelHelper<GFS> & helper_)149 NonoverlappingScalarProduct (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_) 150 : gfs(gfs_), helper(helper_) 151 {} 152 153 /*! \brief Dot product of two vectors. 154 It is assumed that the vectors are consistent on the interior+border 155 partition. 156 */ dot(const X & x,const X & y) const157 virtual field_type dot (const X& x, const X& y) const override 158 { 159 // do local scalar product on unique partition 160 field_type sum = helper.disjointDot(x,y); 161 162 // do global communication 163 return gfs.gridView().comm().sum(sum); 164 } 165 166 /*! \brief Norm of a right-hand side vector. 167 The vector must be consistent on the interior+border partition 168 */ norm(const X & x) const169 virtual double norm (const X& x) const override 170 { 171 return sqrt(static_cast<double>(this->dot(x,x))); 172 } 173 174 /*! \brief make additive vector consistent 175 */ make_consistent(X & x) const176 void make_consistent (X& x) const 177 { 178 Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,x); 179 if (gfs.gridView().comm().size()>1) 180 gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication); 181 } 182 183 private: 184 const GFS& gfs; 185 const ISTL::ParallelHelper<GFS>& helper; 186 }; 187 188 // parallel Richardson preconditioner 189 template<class GFS, class X, class Y> 190 class NonoverlappingRichardson : public Dune::Preconditioner<X,Y> 191 { 192 public: 193 //! \brief The domain type of the preconditioner. 194 typedef X domain_type; 195 //! \brief The range type of the preconditioner. 196 typedef Y range_type; 197 //! \brief The field type of the preconditioner. 198 typedef typename X::ElementType field_type; 199 200 // define the category category() const201 SolverCategory::Category category() const override 202 { 203 return SolverCategory::nonoverlapping; 204 } 205 206 //! \brief Constructor. NonoverlappingRichardson(const GFS & gfs_,const ISTL::ParallelHelper<GFS> & helper_)207 NonoverlappingRichardson (const GFS& gfs_, const ISTL::ParallelHelper<GFS>& helper_) 208 : gfs(gfs_), helper(helper_) 209 { 210 } 211 212 /*! 213 \brief Prepare the preconditioner. 214 */ pre(X & x,Y & b) const215 virtual void pre (X& x, Y& b) const override {} 216 217 /*! 218 \brief Apply the precondioner. 219 */ apply(X & v,const Y & d) const220 virtual void apply (X& v, const Y& d) const override 221 { 222 v = d; 223 } 224 225 /*! 226 \brief Clean up. 227 */ post(X & x)228 virtual void post (X& x) override {} 229 230 private: 231 const GFS& gfs; 232 const ISTL::ParallelHelper<GFS>& helper; 233 }; 234 235 //! parallel non-overlapping Jacobi preconditioner 236 /** 237 * \tparam Diagonal Vector type used to store the diagonal of the matrix 238 * \tparam X Vector type used to store the result of applying the 239 * preconditioner. 240 * \tparam Y Vector type used to store the defect. 241 * \tparam A The matrix type to be used. 242 * 243 * The Jacobi preconditioner approximates the inverse of a matrix M by 244 * taking the diagonal diag(M) and inverting that. In the parallel case 245 * the matrix M is assumed to be inconsistent, so diagonal entries for 246 * dofs on the border are summed up over all relevant processes by this 247 * precoditioner before the inverse is computed. 248 */ 249 template<typename A, typename X, typename Y> 250 class NonoverlappingJacobi 251 : public Dune::Preconditioner<X,Y> 252 { 253 254 typedef typename ISTL::BlockMatrixDiagonal<A>::MatrixElementVector Diagonal; 255 256 Diagonal _inverse_diagonal; 257 258 public: 259 //! The domain type of the operator. 260 /** 261 * The preconditioner is an inverse operator, so this is the output type 262 * of the preconditioner. 263 */ 264 typedef X domain_type; 265 //! \brief The range type of the operator. 266 /** 267 * The preconditioner is an inverse operator, so this is the input type 268 * of the preconditioner. 269 */ 270 typedef Y range_type; 271 //! \brief The field type of the preconditioner. 272 typedef typename X::ElementType field_type; 273 category() const274 SolverCategory::Category category() const override 275 { 276 return SolverCategory::nonoverlapping; 277 } 278 279 //! \brief Constructor. 280 /** 281 * \param gfs The GridFunctionSpace the matrix and the vectors live on. 282 * \param m The matrix whose inverse the preconditioner should 283 * estimate. m is assumed to be inconsistent (i.e. rows for 284 * dofs on the border only contain the contribution of the 285 * local process). 286 * 287 * The preconditioner does not store any reference to the gfs or the 288 * matrix m. The diagonal of m is copied, since it has to be made 289 * consistent. 290 */ 291 template<typename GFS> NonoverlappingJacobi(const GFS & gfs,const A & m)292 NonoverlappingJacobi(const GFS& gfs, const A &m) 293 : _inverse_diagonal(m) 294 { 295 // make the diagonal consistent... 296 typename ISTL::BlockMatrixDiagonal<A>::template AddMatrixElementVectorDataHandle<GFS> addDH(gfs, _inverse_diagonal); 297 gfs.gridView().communicate(addDH, 298 InteriorBorder_InteriorBorder_Interface, 299 ForwardCommunication); 300 301 // ... and then invert it 302 _inverse_diagonal.invert(); 303 304 } 305 306 //! Prepare the preconditioner. pre(X & x,Y & b)307 virtual void pre (X& x, Y& b) override {} 308 309 //! Apply the precondioner. 310 /* 311 * For this preconditioner, this method works with both consistent and 312 * inconsistent vectors: if d is consistent, v will be consistent, if d 313 * is inconsistent, v will be inconsistent. 314 */ apply(X & v,const Y & d)315 virtual void apply (X& v, const Y& d) override 316 { 317 _inverse_diagonal.mv(d,v); 318 } 319 320 //! Clean up. post(X & x)321 virtual void post (X& x) override {} 322 }; 323 324 //! \addtogroup PDELab_novlpsolvers Nonoverlapping Solvers 325 //! \{ 326 327 //! \brief Nonoverlapping parallel CG solver without preconditioner 328 template<class GFS> 329 class ISTLBackend_NOVLP_CG_NOPREC 330 { 331 typedef ISTL::ParallelHelper<GFS> PHELPER; 332 333 public: 334 /*! \brief make a linear solver object 335 336 \param[in] gfs_ a grid function space 337 \param[in] maxiter_ maximum number of iterations to do 338 \param[in] verbose_ print messages if true 339 */ ISTLBackend_NOVLP_CG_NOPREC(const GFS & gfs_,unsigned maxiter_=5000,int verbose_=1)340 explicit ISTLBackend_NOVLP_CG_NOPREC (const GFS& gfs_, 341 unsigned maxiter_=5000, 342 int verbose_=1) 343 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_) 344 {} 345 346 /*! \brief compute global norm of a vector 347 348 \param[in] v the given vector 349 */ 350 template<class V> norm(const V & v) const351 typename V::ElementType norm (const V& v) const 352 { 353 V x(v); // make a copy because it has to be made consistent 354 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 355 PSP psp(gfs,phelper); 356 psp.make_consistent(x); 357 return psp.norm(x); 358 } 359 360 /*! \brief solve the given linear system 361 362 \param[in] A the given matrix 363 \param[out] z the solution vector to be computed 364 \param[in] r right hand side 365 \param[in] reduction to be achieved 366 */ 367 template<class M, class V, class W> apply(M & A,V & z,W & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)368 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 369 { 370 typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP; 371 POP pop(gfs,A); 372 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 373 PSP psp(gfs,phelper); 374 typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH; 375 PRICH prich(gfs,phelper); 376 int verb=0; 377 if (gfs.gridView().comm().rank()==0) verb=verbose; 378 Dune::CGSolver<V> solver(pop,psp,prich,reduction,maxiter,verb); 379 Dune::InverseOperatorResult stat; 380 solver.apply(z,r,stat); 381 res.converged = stat.converged; 382 res.iterations = stat.iterations; 383 res.elapsed = stat.elapsed; 384 res.reduction = stat.reduction; 385 res.conv_rate = stat.conv_rate; 386 } 387 388 /*! \brief Return access to result data */ result() const389 const Dune::PDELab::LinearSolverResult<double>& result() const 390 { 391 return res; 392 } 393 394 private: 395 const GFS& gfs; 396 PHELPER phelper; 397 Dune::PDELab::LinearSolverResult<double> res; 398 unsigned maxiter; 399 int verbose; 400 }; 401 402 //! \brief Nonoverlapping parallel CG solver with Jacobi preconditioner 403 template<class GFS> 404 class ISTLBackend_NOVLP_CG_Jacobi 405 { 406 typedef ISTL::ParallelHelper<GFS> PHELPER; 407 408 const GFS& gfs; 409 PHELPER phelper; 410 LinearSolverResult<double> res; 411 unsigned maxiter; 412 int verbose; 413 414 public: 415 //! make a linear solver object 416 /** 417 * \param gfs_ A grid function space 418 * \param maxiter_ Maximum number of iterations to do. 419 * \param verbose_ Verbosity level, directly handed to the CGSolver. 420 */ ISTLBackend_NOVLP_CG_Jacobi(const GFS & gfs_,unsigned maxiter_=5000,int verbose_=1)421 explicit ISTLBackend_NOVLP_CG_Jacobi(const GFS& gfs_, 422 unsigned maxiter_ = 5000, 423 int verbose_ = 1) : 424 gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_) 425 {} 426 427 //! compute global norm of a vector 428 /** 429 * \param v The vector to compute the norm of. Should be an 430 * inconsistent vector (i.e. the entries corresponding a DoF on 431 * the border should only contain the summand of this process). 432 */ 433 template<class V> norm(const V & v) const434 typename V::ElementType norm (const V& v) const 435 { 436 V x(v); // make a copy because it has to be made consistent 437 typedef NonoverlappingScalarProduct<GFS,V> PSP; 438 PSP psp(gfs,phelper); 439 psp.make_consistent(x); 440 return psp.norm(x); 441 } 442 443 //! solve the given linear system 444 /** 445 * \param A The matrix to solve. Should be a matrix from one of 446 * PDELabs ISTL backends (only ISTLBCRSMatrixBackend at 447 * the moment). 448 * \param z The solution vector to be computed 449 * \param r Right hand side 450 * \param reduction to be achieved 451 * 452 * Solve the linear system A*z=r such that 453 * norm(A*z0-r)/norm(A*z-r) < reduction where z0 is the initial value of 454 * z. 455 */ 456 template<class M, class V, class W> apply(M & A,V & z,W & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)457 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 458 { 459 typedef NonoverlappingOperator<GFS,M,V,W> POP; 460 POP pop(gfs,A); 461 typedef NonoverlappingScalarProduct<GFS,V> PSP; 462 PSP psp(gfs,phelper); 463 464 typedef NonoverlappingJacobi<M,V,W> PPre; 465 PPre ppre(gfs,Backend::native(A)); 466 467 int verb=0; 468 if (gfs.gridView().comm().rank()==0) verb=verbose; 469 CGSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb); 470 InverseOperatorResult stat; 471 solver.apply(z,r,stat); 472 res.converged = stat.converged; 473 res.iterations = stat.iterations; 474 res.elapsed = stat.elapsed; 475 res.reduction = stat.reduction; 476 res.conv_rate = stat.conv_rate; 477 } 478 479 //! Return access to result data result() const480 const LinearSolverResult<double>& result() const 481 { return res; } 482 }; 483 484 //! \brief Nonoverlapping parallel BiCGStab solver without preconditioner 485 template<class GFS> 486 class ISTLBackend_NOVLP_BCGS_NOPREC 487 { 488 typedef ISTL::ParallelHelper<GFS> PHELPER; 489 490 public: 491 /*! \brief make a linear solver object 492 493 \param[in] gfs_ a grid function space 494 \param[in] maxiter_ maximum number of iterations to do 495 \param[in] verbose_ print messages if true 496 */ ISTLBackend_NOVLP_BCGS_NOPREC(const GFS & gfs_,unsigned maxiter_=5000,int verbose_=1)497 explicit ISTLBackend_NOVLP_BCGS_NOPREC (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1) 498 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_) 499 {} 500 501 /*! \brief compute global norm of a vector 502 503 \param[in] v the given vector 504 */ 505 template<class V> norm(const V & v) const506 typename V::ElementType norm (const V& v) const 507 { 508 V x(v); // make a copy because it has to be made consistent 509 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 510 PSP psp(gfs,phelper); 511 psp.make_consistent(x); 512 return psp.norm(x); 513 } 514 515 /*! \brief solve the given linear system 516 517 \param[in] A the given matrix 518 \param[out] z the solution vector to be computed 519 \param[in] r right hand side 520 \param[in] reduction to be achieved 521 */ 522 template<class M, class V, class W> apply(M & A,V & z,W & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)523 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 524 { 525 typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP; 526 POP pop(gfs,A); 527 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 528 PSP psp(gfs,phelper); 529 typedef Dune::PDELab::NonoverlappingRichardson<GFS,V,W> PRICH; 530 PRICH prich(gfs,phelper); 531 int verb=0; 532 if (gfs.gridView().comm().rank()==0) verb=verbose; 533 Dune::BiCGSTABSolver<V> solver(pop,psp,prich,reduction,maxiter,verb); 534 Dune::InverseOperatorResult stat; 535 solver.apply(z,r,stat); 536 res.converged = stat.converged; 537 res.iterations = stat.iterations; 538 res.elapsed = stat.elapsed; 539 res.reduction = stat.reduction; 540 res.conv_rate = stat.conv_rate; 541 } 542 543 /*! \brief Return access to result data */ result() const544 const Dune::PDELab::LinearSolverResult<double>& result() const 545 { 546 return res; 547 } 548 549 private: 550 const GFS& gfs; 551 PHELPER phelper; 552 Dune::PDELab::LinearSolverResult<double> res; 553 unsigned maxiter; 554 int verbose; 555 }; 556 557 558 //! \brief Nonoverlapping parallel BiCGStab solver with Jacobi preconditioner 559 template<class GFS> 560 class ISTLBackend_NOVLP_BCGS_Jacobi 561 { 562 typedef ISTL::ParallelHelper<GFS> PHELPER; 563 564 public: 565 /*! \brief make a linear solver object 566 567 \param[in] gfs_ a grid function space 568 \param[in] maxiter_ maximum number of iterations to do 569 \param[in] verbose_ print messages if true 570 */ ISTLBackend_NOVLP_BCGS_Jacobi(const GFS & gfs_,unsigned maxiter_=5000,int verbose_=1)571 explicit ISTLBackend_NOVLP_BCGS_Jacobi (const GFS& gfs_, unsigned maxiter_=5000, int verbose_=1) 572 : gfs(gfs_), phelper(gfs,verbose_), maxiter(maxiter_), verbose(verbose_) 573 {} 574 575 /*! \brief compute global norm of a vector 576 577 \param[in] v the given vector 578 */ 579 template<class V> norm(const V & v) const580 typename V::ElementType norm (const V& v) const 581 { 582 V x(v); // make a copy because it has to be made consistent 583 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 584 PSP psp(gfs,phelper); 585 psp.make_consistent(x); 586 return psp.norm(x); 587 } 588 589 /*! \brief solve the given linear system 590 591 \param[in] A the given matrix 592 \param[out] z the solution vector to be computed 593 \param[in] r right hand side 594 \param[in] reduction to be achieved 595 */ 596 template<class M, class V, class W> apply(M & A,V & z,W & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)597 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 598 { 599 typedef Dune::PDELab::NonoverlappingOperator<GFS,M,V,W> POP; 600 POP pop(gfs,A); 601 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 602 PSP psp(gfs,phelper); 603 604 typedef NonoverlappingJacobi<M,V,W> PPre; 605 PPre ppre(gfs,A); 606 607 int verb=0; 608 if (gfs.gridView().comm().rank()==0) verb=verbose; 609 Dune::BiCGSTABSolver<V> solver(pop,psp,ppre,reduction,maxiter,verb); 610 Dune::InverseOperatorResult stat; 611 solver.apply(z,r,stat); 612 res.converged = stat.converged; 613 res.iterations = stat.iterations; 614 res.elapsed = stat.elapsed; 615 res.reduction = stat.reduction; 616 res.conv_rate = stat.conv_rate; 617 } 618 619 /*! \brief Return access to result data */ result() const620 const Dune::PDELab::LinearSolverResult<double>& result() const 621 { 622 return res; 623 } 624 625 private: 626 const GFS& gfs; 627 PHELPER phelper; 628 Dune::PDELab::LinearSolverResult<double> res; 629 unsigned maxiter; 630 int verbose; 631 }; 632 633 //! Solver to be used for explicit time-steppers with (block-)diagonal mass matrix 634 template<typename GFS> 635 class ISTLBackend_NOVLP_ExplicitDiagonal 636 { 637 typedef ISTL::ParallelHelper<GFS> PHELPER; 638 639 const GFS& gfs; 640 PHELPER phelper; 641 Dune::PDELab::LinearSolverResult<double> res; 642 643 public: 644 /*! \brief make a linear solver object 645 646 \param[in] gfs_ GridFunctionSpace, used to identify DoFs for parallel 647 communication 648 */ ISTLBackend_NOVLP_ExplicitDiagonal(const GFS & gfs_)649 explicit ISTLBackend_NOVLP_ExplicitDiagonal(const GFS& gfs_) 650 : gfs(gfs_), phelper(gfs) 651 {} 652 653 /*! \brief compute global norm of a vector 654 655 \param[in] v the given vector 656 */ 657 template<class V> norm(const V & v) const658 typename V::ElementType norm (const V& v) const 659 { 660 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 661 V x(v); // make a copy because it has to be made consistent 662 PSP psp(gfs,phelper); 663 psp.make_consistent(x); 664 return psp.norm(x); 665 } 666 667 /*! \brief solve the given linear system 668 669 \param[in] A the given matrix 670 \param[out] z the solution vector to be computed 671 \param[in] r right hand side 672 \param[in] reduction to be achieved 673 */ 674 template<class M, class V, class W> apply(M & A,V & z,W & r,typename Dune::template FieldTraits<typename W::ElementType>::real_type reduction)675 void apply(M& A, V& z, W& r, typename Dune::template FieldTraits<typename W::ElementType >::real_type reduction) 676 { 677 Dune::SeqJac<M,V,W> jac(A,1,1.0); 678 jac.pre(z,r); 679 jac.apply(z,r); 680 jac.post(z); 681 if (gfs.gridView().comm().size()>1) 682 { 683 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,z); 684 gfs.gridView().communicate(adddh,Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication); 685 } 686 res.converged = true; 687 res.iterations = 1; 688 res.elapsed = 0.0; 689 res.reduction = reduction; 690 res.conv_rate = reduction; // pow(reduction,1.0/1) 691 } 692 693 /*! \brief Return access to result data */ result() const694 const Dune::PDELab::LinearSolverResult<double>& result() const 695 { 696 return res; 697 } 698 }; 699 //! \} Nonoverlapping Solvers 700 701 702 /*! \brief Utility base class for preconditioned novlp backends. 703 * \tparam GO The type of the grid operator for the spatial discretization. 704 * This class will be used to adjust the discretization matrix. 705 * and extract the trial grid function space. 706 * \tparam Preconditioner The type of preconditioner to use. 707 * \tparam Solver The type of solver to use. 708 */ 709 template<class GO, 710 template<class,class,class,int> class Preconditioner, 711 template<class> class Solver> 712 class ISTLBackend_NOVLP_BASE_PREC 713 { 714 typedef typename GO::Traits::TrialGridFunctionSpace GFS; 715 typedef ISTL::ParallelHelper<GFS> PHELPER; 716 717 public: 718 /*! \brief Constructor. 719 720 \param[in] gfs_ a grid function space 721 \param[in] maxiter_ maximum number of iterations to do 722 \param[in] steps_ number of preconditioner steps to apply as inner iteration 723 \param[in] verbose_ print messages if true 724 */ ISTLBackend_NOVLP_BASE_PREC(const GO & grid_operator,unsigned maxiter_=5000,unsigned steps_=5,int verbose_=1)725 explicit ISTLBackend_NOVLP_BASE_PREC (const GO& grid_operator, unsigned maxiter_ = 5000, unsigned steps_ = 5, int verbose_ = 1) 726 : _grid_operator(grid_operator) 727 , gfs(grid_operator.trialGridFunctionSpace()) 728 , phelper(gfs,verbose_) 729 , maxiter(maxiter_) 730 , steps(steps_) 731 , verbose(verbose_) 732 {} 733 734 /*! \brief Compute global norm of a vector. 735 736 \param[in] v the given vector 737 */ 738 template<class Vector> norm(const Vector & v) const739 typename Vector::ElementType norm (const Vector& v) const 740 { 741 Vector x(v); // make a copy because it has to be made consistent 742 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,Vector> PSP; 743 PSP psp(gfs,phelper); 744 psp.make_consistent(x); 745 return psp.norm(x); 746 } 747 748 /*! \brief Solve the given linear system. 749 750 \param[in] A the given matrix 751 \param[out] z the solution vector to be computed 752 \param[in] r right hand side 753 \param[in] reduction to be achieved 754 */ 755 template<class M, class V, class W> apply(M & A,V & z,W & r,typename V::ElementType reduction)756 void apply(M& A, V& z, W& r, typename V::ElementType reduction) 757 { 758 using MatrixType = Backend::Native<M>; 759 MatrixType& mat = Backend::native(A); 760 using VectorType = Backend::Native<W>; 761 #if HAVE_MPI 762 typedef typename ISTL::CommSelector<96,Dune::MPIHelper::isFake>::type Comm; 763 _grid_operator.make_consistent(A); 764 ISTL::assertParallelUG(gfs.gridView().comm()); 765 Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping); 766 phelper.createIndexSetAndProjectForAMG(mat, oocc); 767 typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother; 768 Smoother smoother(mat, steps, 1.0); 769 typedef Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> PSP; 770 PSP psp(oocc); 771 typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator; 772 Operator oop(mat,oocc); 773 typedef Dune::NonoverlappingBlockPreconditioner<Comm, Smoother> ParSmoother; 774 ParSmoother parsmoother(smoother, oocc); 775 #else 776 typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother; 777 ParSmoother parsmoother(mat, steps, 1.0); 778 typedef Dune::SeqScalarProduct<VectorType> PSP; 779 PSP psp; 780 typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator; 781 Operator oop(mat); 782 #endif 783 int verb=0; 784 if (gfs.gridView().comm().rank()==0) verb=verbose; 785 Solver<VectorType> solver(oop,psp,parsmoother,reduction,maxiter,verb); 786 Dune::InverseOperatorResult stat; 787 //make r consistent 788 if (gfs.gridView().comm().size()>1){ 789 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r); 790 gfs.gridView().communicate(adddh, 791 Dune::InteriorBorder_InteriorBorder_Interface, 792 Dune::ForwardCommunication); 793 } 794 795 solver.apply(z,r,stat); 796 res.converged = stat.converged; 797 res.iterations = stat.iterations; 798 res.elapsed = stat.elapsed; 799 res.reduction = stat.reduction; 800 res.conv_rate = stat.conv_rate; 801 } 802 803 /*! \brief Return access to result data. */ result() const804 const Dune::PDELab::LinearSolverResult<double>& result() const 805 { 806 return res; 807 } 808 809 private: 810 const GO& _grid_operator; 811 const GFS& gfs; 812 PHELPER phelper; 813 Dune::PDELab::LinearSolverResult<double> res; 814 unsigned maxiter; 815 unsigned steps; 816 int verbose; 817 }; 818 819 //! \addtogroup PDELab_novlpsolvers Nonoverlapping Solvers 820 //! \{ 821 822 /** 823 * @brief Nonoverlapping parallel BiCGSTAB solver preconditioned by block SSOR. 824 * @tparam GO The type of the grid operator used for the spatial discretization 825 * (or the fakeGOTraits class for the old grid operator space). It is used 826 * to adjust the discretization matrix and extract the trial grid function space. 827 * 828 * The solver uses a NonoverlappingBlockPreconditioner with underlying 829 * sequential SSOR preconditioner. The crucial step is to add up the matrix entries 830 * corresponding to the border vertices on each process. This is achieved by 831 * performing a VertexExchanger::sumEntries(Matrix&) before constructing the 832 * sequential SSOR. 833 */ 834 template<class GO> 835 class ISTLBackend_NOVLP_BCGS_SSORk 836 : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver> 837 { 838 839 public: 840 /*! \brief make a linear solver object 841 842 \param[in] gfs_ a grid function space 843 \param[in] maxiter_ maximum number of iterations to do 844 \param[in] steps_ number of SSOR steps to apply as inner iteration 845 \param[in] verbose_ print messages if true 846 */ ISTLBackend_NOVLP_BCGS_SSORk(const GO & grid_operator,unsigned maxiter_=5000,int steps_=5,int verbose_=1)847 explicit ISTLBackend_NOVLP_BCGS_SSORk (const GO& grid_operator, unsigned maxiter_=5000, 848 int steps_=5, int verbose_=1) 849 : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_, steps_, verbose_) 850 {} 851 }; 852 853 /** 854 * @brief Nonoverlapping parallel CG solver preconditioned by block SSOR. 855 * @tparam GO The type of the grid operator used for the spatial discretization. 856 * This class will be used to adjust the discretization matrix. 857 * and extract the trial grid function space. 858 */ 859 template<class GO> 860 class ISTLBackend_NOVLP_CG_SSORk 861 : public ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver> 862 { 863 864 public: 865 /*! \brief make a linear solver object 866 867 \param[in] gfs_ a grid function space 868 \param[in] maxiter_ maximum number of iterations to do 869 \param[in] steps_ number of SSOR steps to apply as inner iteration 870 \param[in] verbose_ print messages if true 871 */ ISTLBackend_NOVLP_CG_SSORk(const GO & grid_operator,unsigned maxiter_=5000,int steps_=5,int verbose_=1)872 explicit ISTLBackend_NOVLP_CG_SSORk (const GO& grid_operator, unsigned maxiter_=5000, 873 int steps_=5, int verbose_=1) 874 : ISTLBackend_NOVLP_BASE_PREC<GO,Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_, steps_, verbose_) 875 {} 876 }; 877 //! \} Nonoverlapping Solvers 878 //! \} group Backend 879 880 template<class GO,int s, template<class,class,class,int> class Preconditioner, 881 template<class> class Solver> 882 class ISTLBackend_AMG_NOVLP : public LinearResultStorage 883 { 884 typedef typename GO::Traits::TrialGridFunctionSpace GFS; 885 typedef typename ISTL::ParallelHelper<GFS> PHELPER; 886 typedef typename GO::Traits::Jacobian M; 887 typedef Backend::Native<M> MatrixType; 888 typedef typename GO::Traits::Domain V; 889 typedef Backend::Native<V> VectorType; 890 typedef typename ISTL::CommSelector<s,Dune::MPIHelper::isFake>::type Comm; 891 #if HAVE_MPI 892 typedef Preconditioner<MatrixType,VectorType,VectorType,1> Smoother; 893 typedef Dune::NonoverlappingBlockPreconditioner<Comm,Smoother> ParSmoother; 894 typedef Dune::NonoverlappingSchwarzOperator<MatrixType,VectorType,VectorType,Comm> Operator; 895 #else 896 typedef Preconditioner<MatrixType,VectorType,VectorType,1> ParSmoother; 897 typedef Dune::MatrixAdapter<MatrixType,VectorType,VectorType> Operator; 898 #endif 899 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs; 900 typedef Dune::Amg::AMG<Operator,VectorType,ParSmoother,Comm> AMG; 901 typedef Dune::Amg::Parameters Parameters; 902 903 public: ISTLBackend_AMG_NOVLP(const GO & grid_operator,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)904 ISTLBackend_AMG_NOVLP(const GO& grid_operator, unsigned maxiter_=5000, 905 int verbose_=1, bool reuse_=false, 906 bool usesuperlu_=true) 907 : _grid_operator(grid_operator) 908 , gfs(grid_operator.trialGridFunctionSpace()) 909 , phelper(gfs,verbose_) 910 , maxiter(maxiter_) 911 , params(15,2000,1.2,1.6,Dune::Amg::atOnceAccu) 912 , verbose(verbose_) 913 , reuse(reuse_) 914 , firstapply(true) 915 , usesuperlu(usesuperlu_) 916 { 917 params.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension); 918 params.setDebugLevel(verbose_); 919 #if !HAVE_SUPERLU 920 if (phelper.rank() == 0 && usesuperlu == true) 921 { 922 std::cout << "WARNING: You are using AMG without SuperLU!" 923 << " Please consider installing SuperLU," 924 << " or set the usesuperlu flag to false" 925 << " to suppress this warning." << std::endl; 926 } 927 #endif 928 } 929 930 /*! \brief set AMG parameters 931 932 \param[in] params_ a parameter object of Type Dune::Amg::Parameters 933 */ setParameters(const Parameters & params_)934 void setParameters(const Parameters& params_) 935 { 936 params = params_; 937 } 938 939 /** 940 * @brief Get the parameters describing the behaviuour of AMG. 941 * 942 * The returned object can be adjusted to ones needs and then can be 943 * reset using setParameters. 944 * @return The object holding the parameters of AMG. 945 */ parameters() const946 const Parameters& parameters() const 947 { 948 return params; 949 } 950 951 //! Set whether the AMG should be reused again during call to apply(). setReuse(bool reuse_)952 void setReuse(bool reuse_) 953 { 954 reuse = reuse_; 955 } 956 957 //! Return whether the AMG is reused during call to apply() getReuse() const958 bool getReuse() const 959 { 960 return reuse; 961 } 962 963 964 /*! \brief compute global norm of a vector 965 966 \param[in] v the given vector 967 */ norm(const V & v) const968 typename V::ElementType norm (const V& v) const 969 { 970 V x(v); // make a copy because it has to be made consistent 971 typedef Dune::PDELab::NonoverlappingScalarProduct<GFS,V> PSP; 972 PSP psp(gfs,phelper); 973 psp.make_consistent(x); 974 return psp.norm(x); 975 } 976 apply(M & A,V & z,V & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)977 void apply(M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 978 { 979 Timer watch; 980 MatrixType& mat = Backend::native(A); 981 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<MatrixType, 982 Dune::Amg::FirstDiagonal> > Criterion; 983 #if HAVE_MPI 984 Comm oocc(gfs.gridView().comm(),Dune::SolverCategory::nonoverlapping); 985 _grid_operator.make_consistent(A); 986 phelper.createIndexSetAndProjectForAMG(A, oocc); 987 Dune::NonoverlappingSchwarzScalarProduct<VectorType,Comm> sp(oocc); 988 Operator oop(mat, oocc); 989 #else 990 Comm oocc(gfs.gridView().comm()); 991 Operator oop(mat); 992 Dune::SeqScalarProduct<VectorType> sp; 993 #endif 994 SmootherArgs smootherArgs; 995 smootherArgs.iterations = 1; 996 smootherArgs.relaxationFactor = 1; 997 //use noAccu or atOnceAccu 998 Criterion criterion(params); 999 stats.tprepare=watch.elapsed(); 1000 watch.reset(); 1001 1002 int verb=0; 1003 if (gfs.gridView().comm().rank()==0) verb=verbose; 1004 //only construct a new AMG if the matrix changes 1005 if (reuse==false || firstapply==true){ 1006 amg.reset(new AMG(oop, criterion, smootherArgs, oocc)); 1007 firstapply = false; 1008 stats.tsetup = watch.elapsed(); 1009 stats.levels = amg->maxlevels(); 1010 stats.directCoarseLevelSolver=amg->usesDirectCoarseLevelSolver(); 1011 } 1012 1013 Dune::InverseOperatorResult stat; 1014 // make r consistent 1015 if (gfs.gridView().comm().size()>1) { 1016 Dune::PDELab::AddDataHandle<GFS,V> adddh(gfs,r); 1017 gfs.gridView().communicate(adddh, 1018 Dune::InteriorBorder_InteriorBorder_Interface, 1019 Dune::ForwardCommunication); 1020 } 1021 watch.reset(); 1022 Solver<VectorType> solver(oop,sp,*amg,reduction,maxiter,verb); 1023 solver.apply(Backend::native(z),Backend::native(r),stat); 1024 stats.tsolve= watch.elapsed(); 1025 res.converged = stat.converged; 1026 res.iterations = stat.iterations; 1027 res.elapsed = stat.elapsed; 1028 res.reduction = stat.reduction; 1029 res.conv_rate = stat.conv_rate; 1030 } 1031 1032 /** 1033 * @brief Get statistics of the AMG solver (no of levels, timings). 1034 * @return statistis of the AMG solver. 1035 */ statistics() const1036 const ISTLAMGStatistics& statistics() const 1037 { 1038 return stats; 1039 } 1040 1041 private: 1042 const GO& _grid_operator; 1043 const GFS& gfs; 1044 PHELPER phelper; 1045 unsigned maxiter; 1046 Parameters params; 1047 int verbose; 1048 bool reuse; 1049 bool firstapply; 1050 bool usesuperlu; 1051 std::shared_ptr<AMG> amg; 1052 ISTLAMGStatistics stats; 1053 }; 1054 1055 //! \addtogroup PDELab_novlpsolvers Nonoverlapping Solvers 1056 //! \{ 1057 1058 /** 1059 * @brief Nonoverlapping parallel CG solver preconditioned with AMG smoothed by SSOR. 1060 * @tparam GO The type of the grid operator 1061 * (or the fakeGOTraits class for the old grid operator space). 1062 * This class will be used to adjust the discretization matrix. 1063 * and extract the trial grid function space. 1064 * @tparam s The bits to use for the global index. 1065 * 1066 * The solver uses AMG with underlying 1067 * sequential SSOR preconditioner. The crucial step is to add up the matrix entries 1068 * corresponding to the border vertices on each process. This is achieved by 1069 * performing a VertexExchanger::sumEntries(Matrix&). 1070 */ 1071 1072 template<class GO, int s=96> 1073 class ISTLBackend_NOVLP_CG_AMG_SSOR 1074 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver> 1075 { 1076 1077 public: ISTLBackend_NOVLP_CG_AMG_SSOR(const GO & grid_operator,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)1078 ISTLBackend_NOVLP_CG_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000, 1079 int verbose_=1, bool reuse_=false, 1080 bool usesuperlu_=true) 1081 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::CGSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_) 1082 {} 1083 }; 1084 1085 /** 1086 * @brief Nonoverlapping parallel BiCGStab solver preconditioned with AMG smoothed by SSOR. 1087 * @tparam GO The type of the grid operator 1088 * (or the fakeGOTraits class for the old grid operator space). 1089 * This class will be used to adjust the discretization matrix. 1090 * and extract the trial grid function space. 1091 * @tparam s The bits to use for the global index. 1092 * 1093 * The solver uses AMG with underlying 1094 * sequential SSOR preconditioner. The crucial step is to add up the matrix entries 1095 * corresponding to the border vertices on each process. This is achieved by 1096 * performing a VertexExchanger::sumEntries(Matrix&). 1097 */ 1098 1099 template<class GO, int s=96> 1100 class ISTLBackend_NOVLP_BCGS_AMG_SSOR 1101 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver> 1102 { 1103 1104 public: ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO & grid_operator,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)1105 ISTLBackend_NOVLP_BCGS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000, 1106 int verbose_=1, bool reuse_=false, 1107 bool usesuperlu_=true) 1108 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::BiCGSTABSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_) 1109 {} 1110 }; 1111 1112 /** 1113 * @brief Nonoverlapping parallel LoopSolver preconditioned with AMG smoothed by SSOR. 1114 * @tparam GO The type of the grid operator used for the spatial discretization 1115 * (or the fakeGOTraits class for the old grid operator space). 1116 * This class will be used to adjust the discretization matrix. 1117 * and extract the trial grid function space. 1118 * @tparam s The bits to use for the global index. 1119 * 1120 * The solver uses AMG with underlying 1121 * sequential SSOR preconditioner. The crucial step is to add up the matrix entries 1122 * corresponding to the border vertices on each process. This is achieved by 1123 * performing a VertexExchanger::sumEntries(Matrix&). 1124 */ 1125 1126 template<class GO, int s=96> 1127 class ISTLBackend_NOVLP_LS_AMG_SSOR 1128 : public ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver> 1129 { 1130 1131 public: ISTLBackend_NOVLP_LS_AMG_SSOR(const GO & grid_operator,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)1132 ISTLBackend_NOVLP_LS_AMG_SSOR(const GO& grid_operator, unsigned maxiter_=5000, 1133 int verbose_=1, bool reuse_=false, 1134 bool usesuperlu_=true) 1135 : ISTLBackend_AMG_NOVLP<GO, s, Dune::SeqSSOR, Dune::LoopSolver>(grid_operator, maxiter_,verbose_,reuse_,usesuperlu_) 1136 {} 1137 }; 1138 //! \} Nonoverlapping Solvers 1139 //! \} group Backend 1140 1141 } // namespace PDELab 1142 } // namespace Dune 1143 1144 #endif // DUNE_PDELAB_BACKEND_ISTL_NOVLPISTLSOLVERBACKEND_HH 1145