1 #ifndef DUNE_PDELAB_SEQ_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_SEQ_AMG_DG_BACKEND_HH
4 #include <dune/common/power.hh>
5 #include <dune/common/parametertree.hh>
7 #include <dune/istl/matrixmatrix.hh>
9 #include <dune/grid/common/datahandleif.hh>
27 template<
class DGMatrix,
class DGPrec,
class CGPrec,
class P>
28 class SeqDGAMGPrec :
public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
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;
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_),
67 virtual void pre (X& x, Y& b)
82 virtual void apply (X& x,
const Y& b)
89 for (
int i=0; i<n1; i++)
104 cgprec.apply(cgv,cgd);
112 for (
int i=0; i<n2; i++)
145 template<
class DGGO,
class CGGFS,
class TransferLOP,
template<
class,
class,
class,
int>
class DGPrec,
template<
class>
class Solver>
149 typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
152 typedef typename DGGO::Traits::Jacobian M;
153 typedef typename DGGO::Traits::Domain V;
154 typedef typename M::BaseT Matrix;
155 typedef typename V::BaseT Vector;
156 typedef typename Vector::field_type field_type;
160 typedef typename CGV::BaseT CGVector;
165 typedef TransferLOP CGTODGLOP;
168 typedef typename PMatrix::BaseT P;
171 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
172 typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG;
173 typedef ACG CGMatrix;
180 std::size_t low_order_space_entries_per_row;
192 , usesuperlu(usesuperlu_)
193 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
195 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
199 if (usesuperlu ==
true)
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;
210 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
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))
223 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
227 if (usesuperlu ==
true)
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;
239 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
248 typename V::ElementType
norm (
const V& v)
const
260 void apply (M& A, V& z, V& r,
typename V::ElementType reduction)
269 Dune::transposeMatMultMat(ptadg,
native(pmatrix),
native(A));
270 Dune::matMultMat(acg,ptadg,
native(pmatrix));
272 double triple_product_time = watch.elapsed();
273 if (verbose>0) std::cout <<
"=== triple matrix product " << triple_product_time <<
" s" << std::endl;
277 typedef ACG CGMatrix;
278 typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
279 CGOperator cgop(acg);
280 typedef Dune::Amg::Parameters Parameters;
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);
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;
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;
305 Dune::MatrixAdapter<Matrix,Vector,Vector> op(
native(A));
306 DGPrec<Matrix,Vector,Vector,1> dgprec(
native(A),1,1);
308 HybridPrec hybridprec(
native(A),dgprec,amg,
native(pmatrix),2,2);
311 Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
314 Dune::InverseOperatorResult 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;
321 res.
elapsed = amg_solve_time+amg_setup_time+triple_product_time;
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
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
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 ¶ms)
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