dune-pdelab  2.7-git
seq_amg_dg_backend.hh
Go to the documentation of this file.
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 
16 
17 namespace Dune {
18  namespace PDELab {
19 
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 
44  SolverCategory::Category category() const override
45  {
46  return SolverCategory::sequential;
47  }
48 
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 
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 
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 
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 
145  template<class DGGO, class CGGFS, class TransferLOP, template<class,class,class,int> class DGPrec, template<class> class Solver>
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
165  typedef TransferLOP CGTODGLOP; // local operator
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:
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 
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 
276  typename V::ElementType norm (const V& v) const
277  {
278  return Backend::native(v).two_norm();
279  }
280 
285  void setParameters(const Parameters& amg_parameters_)
286  {
287  amg_parameters = amg_parameters_;
288  }
289 
297  const Parameters& parameters() const
298  {
299  return amg_parameters;
300  }
301 
303  void setReuse(bool reuse_)
304  {
305  reuse = reuse_;
306  }
307 
309  bool getReuse() const
310  {
311  return reuse;
312  }
313 
315  void setDGSmootherRelaxation(double relaxation_)
316  {
317  smoother_relaxation = relaxation_;
318  }
319 
321  void setNoDGPreSmoothSteps(int n1_)
322  {
323  n1 = n1_;
324  }
325 
328  {
329  n2 = n2_;
330  }
331 
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
const P & p
Definition: constraints.hh:148
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper.
Definition: backend/interface.hh:176
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:188
Definition: seq_amg_dg_backend.hh:30
DGPrec::domain_type X
Definition: seq_amg_dg_backend.hh:39
CGPrec::range_type CGY
Definition: seq_amg_dg_backend.hh:42
CGPrec::domain_type CGX
Definition: seq_amg_dg_backend.hh:41
SeqDGAMGPrec(DGMatrix &dgmatrix_, DGPrec &dgprec_, CGPrec &cgprec_, P &p_, int n1_, int n2_)
Constructor.
Definition: seq_amg_dg_backend.hh:56
DGPrec::range_type Y
Definition: seq_amg_dg_backend.hh:40
virtual void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: seq_amg_dg_backend.hh:67
virtual void apply(X &x, const Y &b) override
Apply the precondioner.
Definition: seq_amg_dg_backend.hh:82
virtual void post(X &x) override
Clean up.
Definition: seq_amg_dg_backend.hh:126
SolverCategory::Category category() const override
Definition: seq_amg_dg_backend.hh:44
Definition: seq_amg_dg_backend.hh:147
void setNoDGPreSmoothSteps(int n1_)
set number of presmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:321
void setNoDGPostSmoothSteps(int n2_)
set number of postsmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:327
void setParameters(const Parameters &amg_parameters_)
set AMG parameters
Definition: seq_amg_dg_backend.hh:285
void setDGSmootherRelaxation(double relaxation_)
set number of presmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:315
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: seq_amg_dg_backend.hh:204
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: seq_amg_dg_backend.hh:303
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: seq_amg_dg_backend.hh:309
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: seq_amg_dg_backend.hh:297
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, const ParameterTree &params)
Definition: seq_amg_dg_backend.hh:238
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: seq_amg_dg_backend.hh:276
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: seq_amg_dg_backend.hh:339
RFType conv_rate
Definition: solver.hh:36
bool converged
Definition: solver.hh:32
RFType reduction
Definition: solver.hh:35
unsigned int iterations
Definition: solver.hh:33
double elapsed
Definition: solver.hh:34
Definition: solver.hh:54
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:63
Definition: constraintstransformation.hh:112
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:180
Dune::PDELab::Backend::Matrix< MBE, Domain, Range, field_type > Jacobian
The type of the jacobian.
Definition: gridoperator.hh:47
Definition: recipe-operator-splitting.cc:108