1 #ifndef DUNE_PDELAB_SEQ_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_SEQ_AMG_DG_BACKEND_HH
4 #include <dune/istl/matrixmatrix.hh>
6 #include <dune/grid/common/datahandleif.hh>
24 template<
class DGMatrix,
class DGPrec,
class CGPrec,
class P>
25 class SeqDGAMGPrec :
public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
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;
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_),
86 for (
int i=0; i<n1; i++)
101 cgprec.apply(cgv,cgd);
109 for (
int i=0; i<n2; i++)
142 template<
class DGGO,
class CGGFS,
class TransferLOP,
template<
class,
class,
class,
int>
class DGPrec,
template<
class>
class Solver>
146 typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
149 typedef typename DGGO::Traits::Jacobian M;
150 typedef typename DGGO::Traits::Domain V;
151 typedef typename M::BaseT Matrix;
152 typedef typename V::BaseT Vector;
153 typedef typename Vector::field_type field_type;
157 typedef typename CGV::BaseT CGVector;
162 typedef TransferLOP CGTODGLOP;
165 typedef typename PMatrix::BaseT P;
168 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
169 typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG;
170 typedef ACG CGMatrix;
184 : dggo(dggo_), cggfs(cggfs_), maxiter(maxiter_), verbose(verbose_), usesuperlu(usesuperlu_),
185 cgtodglop(), pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop), pmatrix(pgo)
188 if (usesuperlu ==
true)
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;
199 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
208 typename V::ElementType
norm (
const V& v)
const
220 void apply (M& A, V& z, V& r,
typename V::ElementType reduction)
226 ACG acg(attached_container);
232 double triple_product_time = watch.elapsed();
233 if (verbose>0) std::cout <<
"=== triple matrix product " << triple_product_time <<
" s" << std::endl;
237 typedef ACG CGMatrix;
238 typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
239 CGOperator cgop(acg);
240 typedef Dune::Amg::Parameters Parameters;
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);
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;
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;
271 Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
274 Dune::InverseOperatorResult 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;
281 res.
elapsed = amg_solve_time+amg_setup_time+triple_product_time;
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
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