1 #ifndef DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_BACKEND_ISTL_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>
28 template<
class DGMatrix,
class DGPrec,
class CGPrec,
class P>
29 class SeqDGAMGPrec :
public Dune::Preconditioner<typename DGPrec::domain_type,typename DGPrec::range_type>
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;
44 SolverCategory::Category
category()
const override
46 return SolverCategory::sequential;
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)
override
82 virtual void apply (
X& x,
const Y& b)
override
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;
156 typedef typename Vector::field_type field_type;
165 typedef TransferLOP CGTODGLOP;
171 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
172 typedef typename Dune::MatMultMatResult<PTADG,P>::type ACG;
173 typedef ACG CGMatrix;
176 typedef Dune::MatrixAdapter<CGMatrix,CGVector,CGVector> CGOperator;
177 typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
178 typedef Dune::Amg::AMG<CGOperator,CGVector,Smoother> AMG;
179 typedef Dune::Amg::Parameters Parameters;
183 std::shared_ptr<CGOperator> cgop;
184 std::shared_ptr<AMG> amg;
185 Parameters amg_parameters;
191 std::size_t low_order_space_entries_per_row;
199 ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_,
unsigned maxiter_=5000,
int verbose_=1,
bool reuse_=
false,
bool usesuperlu_=
true)
202 , amg_parameters(15,2000)
207 , usesuperlu(usesuperlu_)
208 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
210 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,
MBE(low_order_space_entries_per_row))
214 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
215 amg_parameters.setDebugLevel(verbose_);
217 if (usesuperlu ==
true)
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;
228 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
236 , amg_parameters(15,2000)
237 , maxiter(params.get<int>(
"max_iterations",5000))
238 , verbose(params.get<int>(
"verbose",1))
239 , reuse(params.get<bool>(
"reuse", false))
241 , usesuperlu(params.get<bool>(
"use_superlu",true))
242 , low_order_space_entries_per_row(params.get<std::size_t>(
"low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
244 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,
MBE(low_order_space_entries_per_row))
247 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
248 amg_parameters.setDebugLevel(params.get<
int>(
"verbose",1));
250 if (usesuperlu ==
true)
252 std::cout <<
"WARNING: You are using AMG without SuperLU!"
253 <<
" Please consider installing SuperLU,"
254 <<
" or set the usesuperlu flag to false"
255 <<
" to suppress this warning." << std::endl;
262 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
271 typename V::ElementType
norm (
const V& v)
const
282 amg_parameters = amg_parameters_;
294 return amg_parameters;
316 void apply (M& A, V& z, V& r,
typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
323 double triple_product_time = 0.0;
325 if(reuse ==
false || firstapply ==
true) {
327 Dune::transposeMatMultMat(ptadg,
native(pmatrix),
native(A));
328 Dune::matMultMat(acg,ptadg,
native(pmatrix));
329 triple_product_time = watch.elapsed();
331 std::cout <<
"=== triple matrix product " << triple_product_time <<
" s" << std::endl;
335 std::cout <<
"=== reuse CG matrix, SKIPPING triple matrix product " << std::endl;
338 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
339 SmootherArgs smootherArgs;
340 smootherArgs.iterations = 1;
341 smootherArgs.relaxationFactor = 1.0;
342 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
343 Criterion criterion(amg_parameters);
347 double amg_setup_time = 0.0;
348 if(reuse ==
false || firstapply ==
true) {
349 cgop.reset(
new CGOperator(acg));
350 amg.reset(
new AMG(*cgop,criterion,smootherArgs));
352 amg_setup_time = watch.elapsed();
354 std::cout <<
"=== AMG setup " <<amg_setup_time <<
" s" << std::endl;
357 std::cout <<
"=== reuse CG matrix, SKIPPING AMG setup " << std::endl;
360 Dune::MatrixAdapter<Matrix,Vector,Vector> op(
native(A));
361 DGPrec<Matrix,Vector,Vector,1> dgprec(
native(A),1,1);
363 HybridPrec hybridprec(
native(A),dgprec,*amg,
native(pmatrix),2,2);
366 Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
369 Dune::InverseOperatorResult stat;
372 double amg_solve_time = watch.elapsed();
373 if (verbose>0) std::cout <<
"=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time <<
" s" << std::endl;
376 res.
elapsed = amg_solve_time+amg_setup_time+triple_product_time;
384 #endif // DUNE_PDELAB_BACKEND_ISTL_SEQ_AMG_DG_BACKEND_HH