1 #ifndef DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH 2 #define DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH 3 4 #include <dune/common/power.hh> 5 #include <dune/common/parametertree.hh> 6 7 #include <dune/istl/matrixmatrix.hh> 8 9 #include <dune/grid/common/datahandleif.hh> 10 11 #include <dune/pdelab/backend/istl/vector.hh> 12 #include <dune/pdelab/backend/istl/bcrsmatrix.hh> 13 #include <dune/pdelab/backend/istl/bcrsmatrixbackend.hh> 14 #include <dune/pdelab/backend/istl/ovlpistlsolverbackend.hh> 15 #include <dune/pdelab/gridoperator/gridoperator.hh> 16 17 namespace Dune { 18 namespace PDELab { 19 20 /** An ISTL preconditioner for DG based on AMG applied to CG subspace 21 22 The template parameters are: 23 DGMatrix BCRSMatrix assembled with DG 24 DGPrec preconditioner to be used for DG 25 CGPrec preconditioner to be used on CG subspace 26 P BCRSMatrix for grid transfer 27 */ 28 template<class DGMatrix, class DGPrec, class CGPrec, class P> 29 class SeqDGAMGPrec : public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type> 30 { 31 DGMatrix& dgmatrix; 32 DGPrec& dgprec; 33 CGPrec& cgprec; 34 P& p; 35 int n1,n2; 36 bool firstapply; 37 38 public: 39 typedef typename DGPrec::domain_type X; 40 typedef typename DGPrec::range_type Y; 41 typedef typename CGPrec::domain_type CGX; 42 typedef typename CGPrec::range_type CGY; 43 category() const44 SolverCategory::Category category() const override 45 { 46 return SolverCategory::sequential; 47 } 48 49 /*! \brief Constructor. 50 51 Constructor gets all parameters to operate the prec. 52 \param A The matrix to operate on. 53 \param n The number of iterations to perform. 54 \param w The relaxation factor. 55 */ SeqDGAMGPrec(DGMatrix & dgmatrix_,DGPrec & dgprec_,CGPrec & cgprec_,P & p_,int n1_,int n2_)56 SeqDGAMGPrec (DGMatrix& dgmatrix_, DGPrec& dgprec_, CGPrec& cgprec_, P& p_, int n1_, int n2_) 57 : dgmatrix(dgmatrix_), dgprec(dgprec_), cgprec(cgprec_), p(p_), n1(n1_), n2(n2_), 58 firstapply(true) 59 { 60 } 61 62 /*! 63 \brief Prepare the preconditioner. 64 65 \copydoc Preconditioner::pre(X&,Y&) 66 */ pre(X & x,Y & b)67 virtual void pre (X& x, Y& b) override 68 { 69 dgprec.pre(x,b); 70 CGY cgd(p.M()); 71 cgd = 0.0; 72 CGX cgv(p.M()); 73 cgv = 0.0; 74 cgprec.pre(cgv,cgd); 75 } 76 77 /*! 78 \brief Apply the precondioner. 79 80 \copydoc Preconditioner::apply(X&,const Y&) 81 */ apply(X & x,const Y & b)82 virtual void apply (X& x, const Y& b) override 83 { 84 // need local copies to store defect and solution 85 Y d(b); 86 X v(x); 87 88 // pre-smoothing on DG matrix 89 for (int i=0; i<n1; i++) 90 { 91 v = 0.0; 92 dgprec.apply(v,d); 93 dgmatrix.mmv(v,d); 94 x += v; 95 } 96 97 // restrict defect to CG subspace 98 CGY cgd(p.M()); 99 p.mtv(d,cgd); 100 CGX cgv(p.M()); 101 cgv = 0.0; 102 103 // apply AMG 104 cgprec.apply(cgv,cgd); 105 106 // prolongate correction 107 p.mv(cgv,v); 108 dgmatrix.mmv(v,d); 109 x += v; 110 111 // pre-smoothing on DG matrix 112 for (int i=0; i<n2; i++) 113 { 114 v = 0.0; 115 dgprec.apply(v,d); 116 dgmatrix.mmv(v,d); 117 x += v; 118 } 119 } 120 121 /*! 122 \brief Clean up. 123 124 \copydoc Preconditioner::post(X&) 125 */ post(X & x)126 virtual void post (X& x) override 127 { 128 dgprec.post(x); 129 CGX cgv(p.M()); 130 cgv = 0.0; 131 cgprec.post(cgv); 132 } 133 }; 134 135 136 /** Sequential solver backend for using AMG for DG in PDELab 137 138 The template parameters are: 139 DGGO GridOperator for DG discretization, allows access to matrix, vector and grid function space 140 CGGFS grid function space for CG subspace 141 DGPrec preconditioner for DG problem 142 Solver solver to be used on the complete problem 143 144 */ 145 template<class DGGO, class CGGFS, class TransferLOP, template<class,class,class,int> class DGPrec, template<class> class Solver> 146 class ISTLBackend_SEQ_AMG_4_DG : public Dune::PDELab::LinearResultStorage 147 { 148 // DG grid function space 149 typedef typename DGGO::Traits::TrialGridFunctionSpace GFS; 150 151 // vectors and matrices on DG level 152 typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix 153 typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector 154 typedef Backend::Native<M> Matrix; // istl DG matrix 155 typedef Backend::Native<V> Vector; // istl DG vector 156 typedef typename Vector::field_type field_type; 157 158 // vectors and matrices on CG level 159 using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector 160 typedef Backend::Native<CGV> CGVector; // istl CG vector 161 162 // prolongation matrix 163 typedef Dune::PDELab::ISTL::BCRSMatrixBackend<> MBE; 164 typedef Dune::PDELab::EmptyTransformation CC; 165 typedef TransferLOP CGTODGLOP; // local operator 166 typedef Dune::PDELab::GridOperator<CGGFS,GFS,CGTODGLOP,MBE,field_type,field_type,field_type,CC,CC> PGO; 167 typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix 168 typedef Backend::Native<PMatrix> P; // ISTL prolongation matrix 169 170 // CG subspace matrix 171 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG; 172 typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG; // istl coarse space matrix 173 typedef ACG CGMatrix; // another name 174 175 // AMG in CG-subspace 176 typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator; 177 typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother; 178 typedef Dune::Amg::AMG<CGOperator,CGVector,Smoother> AMG; 179 typedef Dune::Amg::Parameters Parameters; 180 181 DGGO& dggo; 182 CGGFS& cggfs; 183 std::shared_ptr<CGOperator> cgop; 184 std::shared_ptr<AMG> amg; 185 Parameters amg_parameters; 186 unsigned maxiter; 187 int verbose; 188 bool reuse; 189 bool firstapply; 190 bool usesuperlu; 191 std::size_t low_order_space_entries_per_row; 192 193 // Parameters for DG smoother 194 int smoother_relaxation = 1.0; // Relaxation parameter of DG smoother 195 int n1=2; // Number of DG pre-smoothing steps 196 int n2=2; // Number of DG post-smoothing steps 197 198 CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix 199 PGO pgo; // grid operator to assemble prolongation matrix 200 PMatrix pmatrix; // wrapped prolongation matrix 201 ACG acg; // CG-subspace matrix 202 203 public: ISTLBackend_SEQ_AMG_4_DG(DGGO & dggo_,CGGFS & cggfs_,unsigned maxiter_=5000,int verbose_=1,bool reuse_=false,bool usesuperlu_=true)204 ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true) 205 : dggo(dggo_) 206 , cggfs(cggfs_) 207 , amg_parameters(15,2000) 208 , maxiter(maxiter_) 209 , verbose(verbose_) 210 , reuse(reuse_) 211 , firstapply(true) 212 , usesuperlu(usesuperlu_) 213 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power) 214 , cgtodglop() 215 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row)) 216 , pmatrix(pgo) 217 , acg() 218 { 219 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension); 220 amg_parameters.setDebugLevel(verbose_); 221 #if !HAVE_SUPERLU 222 if (usesuperlu == true) 223 { 224 std::cout << "WARNING: You are using AMG without SuperLU!" 225 << " Please consider installing SuperLU," 226 << " or set the usesuperlu flag to false" 227 << " to suppress this warning." << std::endl; 228 } 229 #endif 230 231 // assemble prolongation matrix; this will not change from one apply to the next 232 pmatrix = 0.0; 233 if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl; 234 CGV cgx(cggfs,0.0); // need vector to call jacobian 235 pgo.jacobian(cgx,pmatrix); 236 } 237 ISTLBackend_SEQ_AMG_4_DG(DGGO & dggo_,CGGFS & cggfs_,const ParameterTree & params)238 ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, const ParameterTree& params)//unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true) 239 : dggo(dggo_) 240 , cggfs(cggfs_) 241 , amg_parameters(15,2000) 242 , maxiter(params.get<int>("max_iterations",5000)) 243 , verbose(params.get<int>("verbose",1)) 244 , reuse(params.get<bool>("reuse", false)) 245 , firstapply(true) 246 , usesuperlu(params.get<bool>("use_superlu",true)) 247 , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power)) 248 , cgtodglop() 249 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row)) 250 , pmatrix(pgo) 251 { 252 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension); 253 amg_parameters.setDebugLevel(params.get<int>("verbose",1)); 254 #if !HAVE_SUPERLU 255 if (usesuperlu == true) 256 { 257 std::cout << "WARNING: You are using AMG without SuperLU!" 258 << " Please consider installing SuperLU," 259 << " or set the usesuperlu flag to false" 260 << " to suppress this warning." << std::endl; 261 } 262 #endif 263 264 265 // assemble prolongation matrix; this will not change from one apply to the next 266 pmatrix = 0.0; 267 if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl; 268 CGV cgx(cggfs,0.0); // need vector to call jacobian 269 pgo.jacobian(cgx,pmatrix); 270 } 271 272 /*! \brief compute global norm of a vector 273 274 \param[in] v the given vector 275 */ norm(const V & v) const276 typename V::ElementType norm (const V& v) const 277 { 278 return Backend::native(v).two_norm(); 279 } 280 281 /*! \brief set AMG parameters 282 283 \param[in] amg_parameters_ a parameter object of Type Dune::Amg::Parameters 284 */ setParameters(const Parameters & amg_parameters_)285 void setParameters(const Parameters& amg_parameters_) 286 { 287 amg_parameters = amg_parameters_; 288 } 289 290 /** 291 * @brief Get the parameters describing the behaviuour of AMG. 292 * 293 * The returned object can be adjusted to ones needs and then can be 294 * reset using setParameters. 295 * @return The object holding the parameters of AMG. 296 */ parameters() const297 const Parameters& parameters() const 298 { 299 return amg_parameters; 300 } 301 302 //! Set whether the AMG should be reused again during call to apply(). setReuse(bool reuse_)303 void setReuse(bool reuse_) 304 { 305 reuse = reuse_; 306 } 307 308 //! Return whether the AMG is reused during call to apply() getReuse() const309 bool getReuse() const 310 { 311 return reuse; 312 } 313 314 //! set number of presmoothing steps on the DG level setDGSmootherRelaxation(double relaxation_)315 void setDGSmootherRelaxation(double relaxation_) 316 { 317 smoother_relaxation = relaxation_; 318 } 319 320 //! set number of presmoothing steps on the DG level setNoDGPreSmoothSteps(int n1_)321 void setNoDGPreSmoothSteps(int n1_) 322 { 323 n1 = n1_; 324 } 325 326 //! set number of postsmoothing steps on the DG level setNoDGPostSmoothSteps(int n2_)327 void setNoDGPostSmoothSteps(int n2_) 328 { 329 n2 = n2_; 330 } 331 332 /*! \brief solve the given linear system 333 334 \param[in] A the given matrix 335 \param[out] z the solution vector to be computed 336 \param[in] r right hand side 337 \param[in] reduction to be achieved 338 */ apply(M & A,V & z,V & r,typename Dune::template FieldTraits<typename V::ElementType>::real_type reduction)339 void apply (M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction) 340 { 341 using Backend::native; 342 // do triple matrix product ACG = P^T ADG P 343 Dune::Timer watch; 344 watch.reset(); 345 // only do triple matrix product if the matrix changes 346 double triple_product_time = 0.0; 347 // no need to set acg here back to zero, this is done in matMultmat 348 if(reuse == false || firstapply == true) { 349 PTADG ptadg; 350 Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A)); 351 Dune::matMultMat(acg,ptadg,native(pmatrix)); 352 triple_product_time = watch.elapsed(); 353 if (verbose>0) 354 std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl; 355 //Dune::printmatrix(std::cout,acg,"triple product matrix","row",10,2); 356 } 357 else if (verbose>0) 358 std::cout << "=== reuse CG matrix, SKIPPING triple matrix product " << std::endl; 359 360 // set up AMG solver for the CG subspace 361 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs; 362 SmootherArgs smootherArgs; 363 smootherArgs.iterations = 1; 364 smootherArgs.relaxationFactor = 1.0; 365 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion; 366 Criterion criterion(amg_parameters); 367 watch.reset(); 368 369 // only construct a new AMG for the CG-subspace if the matrix changes 370 double amg_setup_time = 0.0; 371 if(reuse == false || firstapply == true) { 372 cgop.reset(new CGOperator(acg)); 373 amg.reset(new AMG(*cgop,criterion,smootherArgs)); 374 firstapply = false; 375 amg_setup_time = watch.elapsed(); 376 if (verbose>0) 377 std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl; 378 } 379 else if (verbose>0) 380 std::cout << "=== reuse CG matrix, SKIPPING AMG setup " << std::endl; 381 382 // set up hybrid DG/CG preconditioner 383 Dune::MatrixAdapter<Matrix,Vector,Vector> op(native(A)); 384 DGPrec<Matrix,Vector,Vector,1> dgprec(native(A),1,smoother_relaxation); 385 typedef SeqDGAMGPrec<Matrix,DGPrec<Matrix,Vector,Vector,1>,AMG,P> HybridPrec; 386 HybridPrec hybridprec(native(A),dgprec,*amg,native(pmatrix),n1,n2); 387 388 // set up solver 389 Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose); 390 391 // solve 392 Dune::InverseOperatorResult stat; 393 watch.reset(); 394 solver.apply(native(z),native(r),stat); 395 double amg_solve_time = watch.elapsed(); 396 if (verbose>0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl; 397 res.converged = stat.converged; 398 res.iterations = stat.iterations; 399 res.elapsed = amg_solve_time+amg_setup_time+triple_product_time; 400 res.reduction = stat.reduction; 401 res.conv_rate = stat.conv_rate; 402 } 403 404 }; 405 } 406 } 407 #endif // DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH 408