dune-pdelab  2.4-dev
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 
15 
16 namespace Dune {
17  namespace PDELab {
18 
27  template<class DGMatrix, class DGPrec, class CGPrec, class P>
28  class SeqDGAMGPrec : public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
29  {
30  DGMatrix& dgmatrix;
31  DGPrec& dgprec;
32  CGPrec& cgprec;
33  P& p;
34  int n1,n2;
35  bool firstapply;
36 
37  public:
38  typedef typename DGPrec::domain_type X;
39  typedef typename DGPrec::range_type Y;
40  typedef typename CGPrec::domain_type CGX;
41  typedef typename CGPrec::range_type CGY;
42 
43  // define the category
44  enum {
46  category=Dune::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)
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)
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)
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 typename M::BaseT Matrix; // istl DG matrix
155  typedef typename V::BaseT 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 typename CGV::BaseT 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 typename PMatrix::BaseT 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  DGGO& dggo;
176  CGGFS& cggfs;
177  unsigned maxiter;
178  int verbose;
179  bool usesuperlu;
180  std::size_t low_order_space_entries_per_row;
181 
182  CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
183  PGO pgo; // grid operator to assemble prolongation matrix
184  PMatrix pmatrix; // wrapped prolongation matrix
185 
186  public:
187  ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
188  : dggo(dggo_)
189  , cggfs(cggfs_)
190  , maxiter(maxiter_)
191  , verbose(verbose_)
192  , usesuperlu(usesuperlu_)
193  , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
194  , cgtodglop()
195  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
196  , pmatrix(pgo)
197  {
198 #if !HAVE_SUPERLU
199  if (usesuperlu == true)
200  {
201  std::cout << "WARNING: You are using AMG without SuperLU!"
202  << " Please consider installing SuperLU,"
203  << " or set the usesuperlu flag to false"
204  << " to suppress this warning." << std::endl;
205  }
206 #endif
207 
208  // assemble prolongation matrix; this will not change from one apply to the next
209  pmatrix = 0.0;
210  if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
211  CGV cgx(cggfs,0.0); // need vector to call jacobian
212  pgo.jacobian(cgx,pmatrix);
213  }
214 
215  ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, const ParameterTree& params)//unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
216  : dggo(dggo_)
217  , cggfs(cggfs_)
218  , maxiter(params.get<int>("max_iterations",5000))
219  , verbose(params.get<int>("verbose",1))
220  , usesuperlu(params.get<bool>("use_superlu",true))
221  , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
222  , cgtodglop()
223  , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
224  , pmatrix(pgo)
225  {
226 #if !HAVE_SUPERLU
227  if (usesuperlu == true)
228  {
229  std::cout << "WARNING: You are using AMG without SuperLU!"
230  << " Please consider installing SuperLU,"
231  << " or set the usesuperlu flag to false"
232  << " to suppress this warning." << std::endl;
233  }
234 #endif
235 
236 
237  // assemble prolongation matrix; this will not change from one apply to the next
238  pmatrix = 0.0;
239  if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
240  CGV cgx(cggfs,0.0); // need vector to call jacobian
241  pgo.jacobian(cgx,pmatrix);
242  }
243 
248  typename V::ElementType norm (const V& v) const
249  {
250  return Backend::native(v).two_norm();
251  }
252 
260  void apply (M& A, V& z, V& r, typename V::ElementType reduction)
261  {
262  using Backend::native;
263  // do triple matrix product ACG = P^T ADG P
264  Dune::Timer watch;
265  watch.reset();
266  ACG acg;
267  {
268  PTADG ptadg;
269  Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A));
270  Dune::matMultMat(acg,ptadg,native(pmatrix));
271  }
272  double triple_product_time = watch.elapsed();
273  if (verbose>0) std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
274  //Dune::printmatrix(std::cout,acg,"triple product matrix","row",10,2);
275 
276  // set up AMG solver for the CG subspace
277  typedef ACG CGMatrix;
278  typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
279  CGOperator cgop(acg);
280  typedef Dune::Amg::Parameters Parameters; // AMG parameters (might be nice to change from outside)
281  Parameters params(15,2000);
282  params.setDefaultValuesIsotropic(CGGFS::Traits::GridViewType::Traits::Grid::dimension);
283  params.setDebugLevel(verbose);
284  params.setCoarsenTarget(1000);
285  params.setMaxLevel(20);
286  params.setProlongationDampingFactor(1.8);
287  params.setNoPreSmoothSteps(2);
288  params.setNoPostSmoothSteps(2);
289  params.setGamma(1);
290  params.setAdditive(false);
291  typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
292  typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
293  SmootherArgs smootherArgs;
294  smootherArgs.iterations = 2;
295  smootherArgs.relaxationFactor = 1.0;
296  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
297  Criterion criterion(params);
298  typedef Dune::Amg::AMG<CGOperator,CGVector,Smoother> AMG;
299  watch.reset();
300  AMG amg(cgop,criterion,smootherArgs);
301  double amg_setup_time = watch.elapsed();
302  if (verbose>0) std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
303 
304  // set up hybrid DG/CG preconditioner
305  Dune::MatrixAdapter<Matrix,Vector,Vector> op(native(A));
306  DGPrec<Matrix,Vector,Vector,1> dgprec(native(A),1,1);
307  typedef SeqDGAMGPrec<Matrix,DGPrec<Matrix,Vector,Vector,1>,AMG,P> HybridPrec;
308  HybridPrec hybridprec(native(A),dgprec,amg,native(pmatrix),2,2);
309 
310  // set up solver
311  Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
312 
313  // solve
314  Dune::InverseOperatorResult stat;
315  watch.reset();
316  solver.apply(native(z),native(r),stat);
317  double amg_solve_time = watch.elapsed();
318  if (verbose>0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl;
319  res.converged = stat.converged;
320  res.iterations = stat.iterations;
321  res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
322  res.reduction = stat.reduction;
323  res.conv_rate = stat.conv_rate;
324  }
325 
326  };
327  }
328 }
329 #endif
The category the preconditioner is part of.
Definition: seq_amg_dg_backend.hh:46
virtual void apply(X &x, const Y &b)
Apply the precondioner.
Definition: seq_amg_dg_backend.hh:82
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:172
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: seq_amg_dg_backend.hh:67
Definition: seq_amg_dg_backend.hh:28
DGPrec::domain_type X
Definition: seq_amg_dg_backend.hh:38
Definition: solver.hh:42
SeqDGAMGPrec(DGMatrix &dgmatrix_, DGPrec &dgprec_, CGPrec &cgprec_, P &p_, int n1_, int n2_)
Constructor.
Definition: seq_amg_dg_backend.hh:56
CGPrec::range_type CGY
Definition: seq_amg_dg_backend.hh:41
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
Definition: seq_amg_dg_backend.hh:187
Definition: seq_amg_dg_backend.hh:146
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:112
RFType reduction
Definition: solver.hh:35
STL namespace.
virtual void post(X &x)
Clean up.
Definition: seq_amg_dg_backend.hh:126
CGPrec::domain_type CGX
Definition: seq_amg_dg_backend.hh:40
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
solve the given linear system
Definition: seq_amg_dg_backend.hh:260
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, const ParameterTree &params)
Definition: seq_amg_dg_backend.hh:215
unsigned int iterations
Definition: solver.hh:33
DGPrec::range_type Y
Definition: seq_amg_dg_backend.hh:39
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:198
bool converged
Definition: solver.hh:32
Dune::PDELab::Backend::Matrix< MBE, Domain, Range, field_type > Jacobian
The type of the jacobian.
Definition: gridoperator.hh:46
RFType conv_rate
Definition: solver.hh:36
double elapsed
Definition: solver.hh:34
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: seq_amg_dg_backend.hh:248
Definition: adaptivity.hh:27
Definition: constraintstransformation.hh:111
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:187
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52