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;
183 std::shared_ptr<CGOperator> cgop;
184 std::shared_ptr<AMG> amg;
191 std::size_t low_order_space_entries_per_row;
194 int smoother_relaxation = 1.0;
204 ISTLBackend_SEQ_AMG_4_DG(DGGO& dggo_, CGGFS& cggfs_,
unsigned maxiter_=5000,
int verbose_=1,
bool reuse_=
false,
bool usesuperlu_=
true)
207 , amg_parameters(15,2000)
212 , usesuperlu(usesuperlu_)
213 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
215 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,
MBE(low_order_space_entries_per_row))
219 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
220 amg_parameters.setDebugLevel(verbose_);
222 if (usesuperlu ==
true)
224 std::cout <<
"WARNING: You are using AMG without SuperLU!"
225 <<
" Please consider installing SuperLU,"
226 <<
" or set the usesuperlu flag to false"
227 <<
" to suppress this warning." << std::endl;
233 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
241 , amg_parameters(15,2000)
242 , maxiter(params.get<int>(
"max_iterations",5000))
243 , verbose(params.get<int>(
"verbose",1))
244 , reuse(params.get<bool>(
"reuse", false))
246 , usesuperlu(params.get<bool>(
"use_superlu",true))
247 , low_order_space_entries_per_row(params.get<
std::size_t>(
"low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
249 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,
MBE(low_order_space_entries_per_row))
252 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
253 amg_parameters.setDebugLevel(params.get<
int>(
"verbose",1));
255 if (usesuperlu ==
true)
257 std::cout <<
"WARNING: You are using AMG without SuperLU!"
258 <<
" Please consider installing SuperLU,"
259 <<
" or set the usesuperlu flag to false"
260 <<
" to suppress this warning." << std::endl;
267 if (verbose>0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
276 typename V::ElementType
norm (
const V& v)
const
287 amg_parameters = amg_parameters_;
299 return amg_parameters;
317 smoother_relaxation = relaxation_;
339 void apply (M& A, V& z, V& r,
typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
346 double triple_product_time = 0.0;
348 if(reuse ==
false || firstapply ==
true) {
350 Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A));
351 Dune::matMultMat(acg,ptadg,native(pmatrix));
352 triple_product_time = watch.elapsed();
354 std::cout <<
"=== triple matrix product " << triple_product_time <<
" s" << std::endl;
358 std::cout <<
"=== reuse CG matrix, SKIPPING triple matrix product " << std::endl;
361 typedef typename Dune::Amg::SmootherTraits<Smoother>::Arguments SmootherArgs;
362 SmootherArgs smootherArgs;
363 smootherArgs.iterations = 1;
364 smootherArgs.relaxationFactor = 1.0;
365 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
366 Criterion criterion(amg_parameters);
370 double amg_setup_time = 0.0;
371 if(reuse ==
false || firstapply ==
true) {
372 cgop.reset(
new CGOperator(acg));
373 amg.reset(
new AMG(*cgop,criterion,smootherArgs));
375 amg_setup_time = watch.elapsed();
377 std::cout <<
"=== AMG setup " <<amg_setup_time <<
" s" << std::endl;
380 std::cout <<
"=== reuse CG matrix, SKIPPING AMG setup " << std::endl;
383 Dune::MatrixAdapter<Matrix,Vector,Vector> op(native(A));
384 DGPrec<Matrix,Vector,Vector,1> dgprec(native(A),1,smoother_relaxation);
386 HybridPrec hybridprec(native(A),dgprec,*amg,native(pmatrix),n1,n2);
389 Solver<Vector> solver(op,hybridprec,reduction,maxiter,verbose);
392 Dune::InverseOperatorResult stat;
394 solver.apply(native(z),native(r),stat);
395 double amg_solve_time = watch.elapsed();
396 if (verbose>0) std::cout <<
"=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time <<
" s" << std::endl;
399 res.
elapsed = amg_solve_time+amg_setup_time+triple_product_time;
const P & p
Definition: constraints.hh:148
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:192
typename native_type< T >::type Native
Alias of the native container type associated with T or T itself if it is not a backend wrapper.
Definition: backend/interface.hh:176
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:188
Definition: seq_amg_dg_backend.hh:30
DGPrec::domain_type X
Definition: seq_amg_dg_backend.hh:39
CGPrec::range_type CGY
Definition: seq_amg_dg_backend.hh:42
CGPrec::domain_type CGX
Definition: seq_amg_dg_backend.hh:41
SeqDGAMGPrec(DGMatrix &dgmatrix_, DGPrec &dgprec_, CGPrec &cgprec_, P &p_, int n1_, int n2_)
Constructor.
Definition: seq_amg_dg_backend.hh:56
DGPrec::range_type Y
Definition: seq_amg_dg_backend.hh:40
virtual void pre(X &x, Y &b) override
Prepare the preconditioner.
Definition: seq_amg_dg_backend.hh:67
virtual void apply(X &x, const Y &b) override
Apply the precondioner.
Definition: seq_amg_dg_backend.hh:82
virtual void post(X &x) override
Clean up.
Definition: seq_amg_dg_backend.hh:126
SolverCategory::Category category() const override
Definition: seq_amg_dg_backend.hh:44
Definition: seq_amg_dg_backend.hh:147
void setNoDGPreSmoothSteps(int n1_)
set number of presmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:321
void setNoDGPostSmoothSteps(int n2_)
set number of postsmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:327
void setParameters(const Parameters &amg_parameters_)
set AMG parameters
Definition: seq_amg_dg_backend.hh:285
void setDGSmootherRelaxation(double relaxation_)
set number of presmoothing steps on the DG level
Definition: seq_amg_dg_backend.hh:315
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: seq_amg_dg_backend.hh:204
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: seq_amg_dg_backend.hh:303
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: seq_amg_dg_backend.hh:297
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: seq_amg_dg_backend.hh:309
ISTLBackend_SEQ_AMG_4_DG(DGGO &dggo_, CGGFS &cggfs_, const ParameterTree ¶ms)
Definition: seq_amg_dg_backend.hh:238
V::ElementType norm(const V &v) const
compute global norm of a vector
Definition: seq_amg_dg_backend.hh:276
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: seq_amg_dg_backend.hh:339
RFType conv_rate
Definition: solver.hh:36
bool converged
Definition: solver.hh:32
RFType reduction
Definition: solver.hh:35
unsigned int iterations
Definition: solver.hh:33
double elapsed
Definition: solver.hh:34
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:63
Definition: constraintstransformation.hh:112
Standard grid operator implementation.
Definition: gridoperator.hh:36
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:180
Dune::PDELab::Backend::Matrix< MBE, Domain, Range, field_type > Jacobian
The type of the jacobian.
Definition: gridoperator.hh:47
Definition: recipe-operator-splitting.cc:108