dune-pdelab  2.5-dev
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 // this is here for backwards compatibility and deprecation warnings, remove after 2.5.0
5 #include "ensureistlinclude.hh"
6 
7 #include <dune/common/power.hh>
8 #include <dune/common/parametertree.hh>
9 
10 #include <dune/istl/matrixmatrix.hh>
11 
12 #include <dune/grid/common/datahandleif.hh>
13 
19 
20 namespace Dune {
21  namespace PDELab {
22 
31  template<class DGMatrix, class DGPrec, class CGPrec, class P>
32  class SeqDGAMGPrec : public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
33  {
34  DGMatrix& dgmatrix;
35  DGPrec& dgprec;
36  CGPrec& cgprec;
37  P& p;
38  int n1,n2;
39  bool firstapply;
40 
41  public:
42  typedef typename DGPrec::domain_type X;
43  typedef typename DGPrec::range_type Y;
44  typedef typename CGPrec::domain_type CGX;
45  typedef typename CGPrec::range_type CGY;
46 
47  // define the category
48  enum {
50  category=Dune::SolverCategory::sequential
51  };
52 
60  SeqDGAMGPrec (DGMatrix& dgmatrix_, DGPrec& dgprec_, CGPrec& cgprec_, P& p_, int n1_, int n2_)
61  : dgmatrix(dgmatrix_), dgprec(dgprec_), cgprec(cgprec_), p(p_), n1(n1_), n2(n2_),
62  firstapply(true)
63  {
64  }
65 
71  virtual void pre (X& x, Y& b)
72  {
73  dgprec.pre(x,b);
74  CGY cgd(p.M());
75  cgd = 0.0;
76  CGX cgv(p.M());
77  cgv = 0.0;
78  cgprec.pre(cgv,cgd);
79  }
80 
86  virtual void apply (X& x, const Y& b)
87  {
88  // need local copies to store defect and solution
89  Y d(b);
90  X v(x);
91 
92  // pre-smoothing on DG matrix
93  for (int i=0; i<n1; i++)
94  {
95  v = 0.0;
96  dgprec.apply(v,d);
97  dgmatrix.mmv(v,d);
98  x += v;
99  }
100 
101  // restrict defect to CG subspace
102  CGY cgd(p.M());
103  p.mtv(d,cgd);
104  CGX cgv(p.M());
105  cgv = 0.0;
106 
107  // apply AMG
108  cgprec.apply(cgv,cgd);
109 
110  // prolongate correction
111  p.mv(cgv,v);
112  dgmatrix.mmv(v,d);
113  x += v;
114 
115  // pre-smoothing on DG matrix
116  for (int i=0; i<n2; i++)
117  {
118  v = 0.0;
119  dgprec.apply(v,d);
120  dgmatrix.mmv(v,d);
121  x += v;
122  }
123  }
124 
130  virtual void post (X& x)
131  {
132  dgprec.post(x);
133  CGX cgv(p.M());
134  cgv = 0.0;
135  cgprec.post(cgv);
136  }
137  };
138 
139 
149  template<class DGGO, class CGGFS, class TransferLOP, template<class,class,class,int> class DGPrec, template<class> class Solver>
151  {
152  // DG grid function space
153  typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
154 
155  // vectors and matrices on DG level
156  typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix
157  typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector
158  typedef Backend::Native<M> Matrix; // istl DG matrix
159  typedef Backend::Native<V> Vector; // istl DG vector
160  typedef typename Vector::field_type field_type;
161 
162  // vectors and matrices on CG level
163  using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector
164  typedef Backend::Native<CGV> CGVector; // istl CG vector
165 
166  // prolongation matrix
169  typedef TransferLOP CGTODGLOP; // local operator
171  typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix
172  typedef Backend::Native<PMatrix> P; // ISTL prolongation matrix
173 
174  // CG subspace matrix
175  typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
176  typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG; // istl coarse space matrix
177  typedef ACG CGMatrix; // another name
178 
179  // AMG in CG-subspace
180  typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
181  typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
182  typedef Dune::Amg::AMG<CGOperator,CGVector,Smoother> AMG;
183  typedef Dune::Amg::Parameters Parameters;
184 
185  DGGO& dggo;
186  CGGFS& cggfs;
187  std::shared_ptr<CGOperator> cgop;
188  std::shared_ptr<AMG> amg;
189  Parameters amg_parameters;
190  unsigned maxiter;
191  int verbose;
192  bool reuse;
193  bool firstapply;
194  bool usesuperlu;
195  std::size_t low_order_space_entries_per_row;
196 
197  CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
198  PGO pgo; // grid operator to assemble prolongation matrix
199  PMatrix pmatrix; // wrapped prolongation matrix
200  ACG acg; // CG-subspace matrix
201 
202  public:
203  ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
204  : dggo(dggo_)
205  , cggfs(cggfs_)
206  , amg_parameters(15,2000)
207  , maxiter(maxiter_)
208  , verbose(verbose_)
209  , reuse(reuse_)
210  , firstapply(true)
211  , usesuperlu(usesuperlu_)
212  , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
213  , cgtodglop()
214  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
215  , pmatrix(pgo)
216  , acg()
217  {
218  amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
219  amg_parameters.setDebugLevel(verbose_);
220 #if !HAVE_SUPERLU
221  if (usesuperlu == true)
222  {
223  std::cout << "WARNING: You are using AMG without SuperLU!"
224  << " Please consider installing SuperLU,"
225  << " or set the usesuperlu flag to false"
226  << " to suppress this warning." << std::endl;
227  }
228 #endif
229 
230  // assemble prolongation matrix; this will not change from one apply to the next
231  pmatrix = 0.0;
232  if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
233  CGV cgx(cggfs,0.0); // need vector to call jacobian
234  pgo.jacobian(cgx,pmatrix);
235  }
236 
237  ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, const ParameterTree& params)//unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
238  : dggo(dggo_)
239  , cggfs(cggfs_)
240  , amg_parameters(15,2000)
241  , maxiter(params.get<int>("max_iterations",5000))
242  , verbose(params.get<int>("verbose",1))
243  , usesuperlu(params.get<bool>("use_superlu",true))
244  , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
245  , cgtodglop()
246  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
247  , pmatrix(pgo)
248  {
249  amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
250  amg_parameters.setDebugLevel(params.get<int>("verbose",1));
251 #if !HAVE_SUPERLU
252  if (usesuperlu == true)
253  {
254  std::cout << "WARNING: You are using AMG without SuperLU!"
255  << " Please consider installing SuperLU,"
256  << " or set the usesuperlu flag to false"
257  << " to suppress this warning." << std::endl;
258  }
259 #endif
260 
261 
262  // assemble prolongation matrix; this will not change from one apply to the next
263  pmatrix = 0.0;
264  if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
265  CGV cgx(cggfs,0.0); // need vector to call jacobian
266  pgo.jacobian(cgx,pmatrix);
267  }
268 
273  typename V::ElementType norm (const V& v) const
274  {
275  return Backend::native(v).two_norm();
276  }
277 
282  void setParameters(const Parameters& amg_parameters_)
283  {
284  amg_parameters = amg_parameters_;
285  }
286 
294  const Parameters& parameters() const
295  {
296  return amg_parameters;
297  }
298 
300  void setReuse(bool reuse_)
301  {
302  reuse = reuse_;
303  }
304 
306  bool getReuse() const
307  {
308  return reuse;
309  }
310 
318  void apply (M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
319  {
320  using Backend::native;
321  // do triple matrix product ACG = P^T ADG P
322  Dune::Timer watch;
323  watch.reset();
324  // only do triple matrix product if the matrix changes
325  double triple_product_time = 0.0;
326  // no need to set acg here back to zero, this is done in matMultmat
327  if(reuse == false || firstapply == true) {
328  PTADG ptadg;
329  Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A));
330  Dune::matMultMat(acg,ptadg,native(pmatrix));
331  triple_product_time = watch.elapsed();
332  if (verbose>0)
333  std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
334  //Dune::printmatrix(std::cout,acg,"triple product matrix","row",10,2);
335  }
336  else if (verbose>0)
337  std::cout << "=== reuse CG matrix, SKIPPING triple matrix product " << std::endl;
338 
339  // set up AMG solver for the CG subspace
340  typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
341  SmootherArgs smootherArgs;
342  smootherArgs.iterations = 1;
343  smootherArgs.relaxationFactor = 1.0;
344  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
345  Criterion criterion(amg_parameters);
346  watch.reset();
347 
348  // only construct a new AMG for the CG-subspace if the matrix changes
349  double amg_setup_time = 0.0;
350  if(reuse == false || firstapply == true) {
351  cgop.reset(new CGOperator(acg));
352  amg.reset(new AMG(*cgop,criterion,smootherArgs));
353  firstapply = false;
354  amg_setup_time = watch.elapsed();
355  if (verbose>0)
356  std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
357  }
358  else if (verbose>0)
359  std::cout << "=== reuse CG matrix, SKIPPING AMG setup " << std::endl;
360 
361  // set up hybrid DG/CG preconditioner
362  Dune::MatrixAdapter<Matrix,Vector,Vector> op(native(A));
363  DGPrec<Matrix,Vector,Vector,1> dgprec(native(A),1,1);
364  typedef SeqDGAMGPrec<Matrix,DGPrec<Matrix,Vector,Vector,1>,AMG,P> HybridPrec;
365  HybridPrec hybridprec(native(A),dgprec,*amg,native(pmatrix),2,2);
366 
367  // set up solver
368  Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
369 
370  // solve
371  Dune::InverseOperatorResult stat;
372  watch.reset();
373  solver.apply(native(z),native(r),stat);
374  double amg_solve_time = watch.elapsed();
375  if (verbose>0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl;
376  res.converged = stat.converged;
377  res.iterations = stat.iterations;
378  res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
379  res.reduction = stat.reduction;
380  res.conv_rate = stat.conv_rate;
381  }
382 
383  };
384  }
385 }
386 #endif // DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: seq_amg_dg_backend.hh:71
CGPrec::range_type CGY
Definition: seq_amg_dg_backend.hh:45
STL namespace.
Definition: constraintstransformation.hh:111
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: seq_amg_dg_backend.hh:273
Definition: seq_amg_dg_backend.hh:150
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: seq_amg_dg_backend.hh:300
Definition: seq_amg_dg_backend.hh:32
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: seq_amg_dg_backend.hh:294
typename impl::BackendMatrixSelector< Backend, VU, VV, E >::Type Matrix
alias of the return type of BackendMatrixSelector
Definition: backend/interface.hh:127
The category the preconditioner is part of.
Definition: seq_amg_dg_backend.hh:50
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, const ParameterTree &params)
Definition: seq_amg_dg_backend.hh:237
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:203
virtual void apply(X &x, const Y &b)
Apply the precondioner.
Definition: seq_amg_dg_backend.hh:86
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:190
DGPrec::domain_type X
Definition: seq_amg_dg_backend.hh:42
DGPrec::range_type Y
Definition: seq_amg_dg_backend.hh:43
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: seq_amg_dg_backend.hh:306
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > &>::type native(T &t)
Definition: backend/interface.hh:192
SeqDGAMGPrec(DGMatrix &dgmatrix_, DGPrec &dgprec_, CGPrec &cgprec_, P &p_, int n1_, int n2_)
Constructor.
Definition: seq_amg_dg_backend.hh:60
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
virtual void post(X &x)
Clean up.
Definition: seq_amg_dg_backend.hh:130
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:205
CGPrec::domain_type CGX
Definition: seq_amg_dg_backend.hh:44
Definition: solver.hh:42
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:318
void setParameters(const Parameters &amg_parameters_)
set AMG parameters
Definition: seq_amg_dg_backend.hh:282