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