dune-pdelab  2.0.0
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/istl/matrixmatrix.hh>
5 
6 #include <dune/grid/common/datahandleif.hh>
7 
12 
13 namespace Dune {
14  namespace PDELab {
15 
24  template<class DGMatrix, class DGPrec, class CGPrec, class P>
25  class SeqDGAMGPrec : public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
26  {
27  DGMatrix& dgmatrix;
28  DGPrec& dgprec;
29  CGPrec& cgprec;
30  P& p;
31  int n1,n2;
32  bool firstapply;
33 
34  public:
35  typedef typename DGPrec::domain_type X;
36  typedef typename DGPrec::range_type Y;
37  typedef typename CGPrec::domain_type CGX;
38  typedef typename CGPrec::range_type CGY;
39 
40  // define the category
41  enum {
43  category=Dune::SolverCategory::sequential
44  };
45 
53  SeqDGAMGPrec (DGMatrix& dgmatrix_, DGPrec& dgprec_, CGPrec& cgprec_, P& p_, int n1_, int n2_)
54  : dgmatrix(dgmatrix_), dgprec(dgprec_), cgprec(cgprec_), p(p_), n1(n1_), n2(n2_),
55  firstapply(true)
56  {
57  }
58 
64  virtual void pre (X& x, Y& b)
65  {
66  dgprec.pre(x,b);
67  CGY cgd(p.M());
68  cgd = 0.0;
69  CGX cgv(p.M());
70  cgv = 0.0;
71  cgprec.pre(cgv,cgd);
72  }
73 
79  virtual void apply (X& x, const Y& b)
80  {
81  // need local copies to store defect and solution
82  Y d(b);
83  X v(x);
84 
85  // pre-smoothing on DG matrix
86  for (int i=0; i<n1; i++)
87  {
88  v = 0.0;
89  dgprec.apply(v,d);
90  dgmatrix.mmv(v,d);
91  x += v;
92  }
93 
94  // restrict defect to CG subspace
95  CGY cgd(p.M());
96  p.mtv(d,cgd);
97  CGX cgv(p.M());
98  cgv = 0.0;
99 
100  // apply AMG
101  cgprec.apply(cgv,cgd);
102 
103  // prolongate correction
104  p.mv(cgv,v);
105  dgmatrix.mmv(v,d);
106  x += v;
107 
108  // pre-smoothing on DG matrix
109  for (int i=0; i<n2; i++)
110  {
111  v = 0.0;
112  dgprec.apply(v,d);
113  dgmatrix.mmv(v,d);
114  x += v;
115  }
116  }
117 
123  virtual void post (X& x)
124  {
125  dgprec.post(x);
126  CGX cgv(p.M());
127  cgv = 0.0;
128  cgprec.post(cgv);
129  }
130  };
131 
132 
142  template<class DGGO, class CGGFS, class TransferLOP, template<class,class,class,int> class DGPrec, template<class> class Solver>
144  {
145  // DG grid function space
146  typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
147 
148  // vectors and matrices on DG level
149  typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix
150  typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector
151  typedef typename M::BaseT Matrix; // istl DG matrix
152  typedef typename V::BaseT Vector; // istl DG vector
153  typedef typename Vector::field_type field_type;
154 
155  // vectors and matrices on CG level
156  typedef typename Dune::PDELab::BackendVectorSelector<CGGFS,field_type>::Type CGV; // wrapped istl CG vector
157  typedef typename CGV::BaseT CGVector; // istl CG vector
158 
159  // prolongation matrix
162  typedef TransferLOP CGTODGLOP; // local operator
164  typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix
165  typedef typename PMatrix::BaseT P; // ISTL prolongation matrix
166 
167  // CG subspace matrix
168  typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
169  typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG; // istl coarse space matrix
170  typedef ACG CGMatrix; // another name
171 
172  DGGO& dggo;
173  CGGFS& cggfs;
174  unsigned maxiter;
175  int verbose;
176  bool usesuperlu;
177 
178  CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
179  PGO pgo; // grid operator to assemble prolongation matrix
180  PMatrix pmatrix; // wrapped prolongation matrix
181 
182  public:
183  ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
184  : dggo(dggo_), cggfs(cggfs_), maxiter(maxiter_), verbose(verbose_), usesuperlu(usesuperlu_),
185  cgtodglop(), pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop), pmatrix(pgo)
186  {
187 #if !HAVE_SUPERLU
188  if (usesuperlu == true)
189  {
190  std::cout << "WARNING: You are using AMG without SuperLU!"
191  << " Please consider installing SuperLU,"
192  << " or set the usesuperlu flag to false"
193  << " to suppress this warning." << std::endl;
194  }
195 #endif
196 
197  // assemble prolongation matrix; this will not change from one apply to the next
198  pmatrix = 0.0;
199  if (verbose>0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
200  CGV cgx(cggfs,0.0); // need vector to call jacobian
201  pgo.jacobian(cgx,pmatrix);
202  }
203 
208  typename V::ElementType norm (const V& v) const
209  {
210  return Dune::PDELab::istl::raw(v).two_norm();
211  }
212 
220  void apply (M& A, V& z, V& r, typename V::ElementType reduction)
221  {
222  // do triple matrix product ACG = P^T ADG P
223  Dune::Timer watch;
224  watch.reset();
225  tags::attached_container attached_container;
226  ACG acg(attached_container);
227  {
228  PTADG ptadg;
229  Dune::transposeMatMultMat(ptadg,Dune::PDELab::istl::raw(pmatrix),Dune::PDELab::istl::raw(A));
230  Dune::matMultMat(acg,ptadg,Dune::PDELab::istl::raw(pmatrix));
231  }
232  double triple_product_time = watch.elapsed();
233  if (verbose>0) std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
234  //Dune::printmatrix(std::cout,acg,"triple product matrix","row",10,2);
235 
236  // set up AMG solver for the CG subspace
237  typedef ACG CGMatrix;
238  typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
239  CGOperator cgop(acg);
240  typedef Dune::Amg::Parameters Parameters; // AMG parameters (might be nice to change from outside)
241  Parameters params(15,2000);
242  params.setDefaultValuesIsotropic(CGGFS::Traits::GridViewType::Traits::Grid::dimension);
243  params.setDebugLevel(verbose);
244  params.setCoarsenTarget(1000);
245  params.setMaxLevel(20);
246  params.setProlongationDampingFactor(1.8);
247  params.setNoPreSmoothSteps(2);
248  params.setNoPostSmoothSteps(2);
249  params.setGamma(1);
250  params.setAdditive(false);
251  typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
252  typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
253  SmootherArgs smootherArgs;
254  smootherArgs.iterations = 2;
255  smootherArgs.relaxationFactor = 1.0;
256  typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
257  Criterion criterion(params);
258  typedef Dune::Amg::AMG<CGOperator,CGVector,Smoother> AMG;
259  watch.reset();
260  AMG amg(cgop,criterion,smootherArgs);
261  double amg_setup_time = watch.elapsed();
262  if (verbose>0) std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
263 
264  // set up hybrid DG/CG preconditioner
265  Dune::MatrixAdapter<Matrix,Vector,Vector> op(Dune::PDELab::istl::raw(A));
266  DGPrec<Matrix,Vector,Vector,1> dgprec(Dune::PDELab::istl::raw(A),1,1);
267  typedef SeqDGAMGPrec<Matrix,DGPrec<Matrix,Vector,Vector,1>,AMG,P> HybridPrec;
268  HybridPrec hybridprec(Dune::PDELab::istl::raw(A),dgprec,amg,Dune::PDELab::istl::raw(pmatrix),2,2);
269 
270  // set up solver
271  Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
272 
273  // solve
274  Dune::InverseOperatorResult stat;
275  watch.reset();
276  solver.apply(Dune::PDELab::istl::raw(z),Dune::PDELab::istl::raw(r),stat);
277  double amg_solve_time = watch.elapsed();
278  if (verbose>0) std::cout << "=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time << " s" << std::endl;
279  res.converged = stat.converged;
280  res.iterations = stat.iterations;
281  res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
282  res.reduction = stat.reduction;
283  res.conv_rate = stat.conv_rate;
284  }
285 
286  };
287  }
288 }
289 #endif
virtual void pre(X &x, Y &b)
Prepare the preconditioner.
Definition: seq_amg_dg_backend.hh:64
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:172
Definition: seq_amg_dg_backend.hh:25
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
Definition: seq_amg_dg_backend.hh:183
bool converged
Definition: solver.hh:32
Tag for requesting a vector or matrix container with a pre-attached underlying object.
Definition: backend/tags.hh:27
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: seq_amg_dg_backend.hh:208
CGPrec::range_type CGY
Definition: seq_amg_dg_backend.hh:38
double elapsed
Definition: solver.hh:34
virtual void apply(X &x, const Y &b)
Apply the precondioner.
Definition: seq_amg_dg_backend.hh:79
SeqDGAMGPrec(DGMatrix &dgmatrix_, DGPrec &dgprec_, CGPrec &cgprec_, P &p_, int n1_, int n2_)
Constructor.
Definition: seq_amg_dg_backend.hh:53
Definition: common/constraintstransformation.hh:127
DGPrec::range_type Y
Definition: seq_amg_dg_backend.hh:36
Backend using ISTL matrices.
Definition: istl/descriptors.hh:69
CGPrec::domain_type CGX
Definition: seq_amg_dg_backend.hh:37
DGPrec::domain_type X
Definition: seq_amg_dg_backend.hh:35
The category the preconditioner is part of.
Definition: seq_amg_dg_backend.hh:43
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
solve the given linear system
Definition: seq_amg_dg_backend.hh:220
V & raw(V &v)
Returns the raw ISTL object associated with v, or v itself it is already an ISTL object.
Definition: backend/istl/utility.hh:26
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52
Definition: seq_amg_dg_backend.hh:143
virtual void post(X &x)
Clean up.
Definition: seq_amg_dg_backend.hh:123
RFType reduction
Definition: solver.hh:35
BackendVectorSelectorHelper< Backend, GridFunctionSpace, FieldType >::Type Type
Definition: backendselector.hh:14
Definition: solver.hh:42
RFType conv_rate
Definition: solver.hh:36
Dune::PDELab::BackendMatrixSelector< MBE, Domain, Range, field_type >::Type Jacobian
The type of the jacobian.
Definition: gridoperator.hh:46
unsigned int iterations
Definition: solver.hh:33