1 #ifndef DUNE_PDELAB_OVLP_AMG_DG_BACKEND_HH
2 #define DUNE_PDELAB_OVLP_AMG_DG_BACKEND_HH
4 #include <dune/istl/matrixmatrix.hh>
6 #include <dune/grid/common/datahandleif.hh>
29 :
public Dune::CommDataHandleIF<LocalGlobalMapDataHandle<GFS>,int>
33 typedef typename GFS::Traits::GridView
GV;
36 typedef typename GV::Grid
Grid;
38 typedef typename GlobalIdSet::IdType
IdType;
63 template<
class EntityType>
64 size_t size (EntityType&
e)
const
70 template<
class MessageBuffer,
class EntityType>
71 void gather (MessageBuffer& buff,
const EntityType&
e)
const
74 IndexType myindex = indexSet.index(e);
75 IdType myid = globalIdSet.id(e);
87 template<
class MessageBuffer,
class EntityType>
88 void scatter (MessageBuffer& buff,
const EntityType&
e,
size_t n)
93 IndexType myindex = indexSet.index(e);
94 IdType myid = globalIdSet.id(e);
102 : gfs(gfs_), gv(gfs.gridView()), indexSet(gv.indexSet()),
103 grid(gv.grid()), globalIdSet(grid.globalIdSet()),
111 const IndexSet& indexSet;
113 const GlobalIdSet& globalIdSet;
114 LocalToGlobalMap& l2g;
115 GlobalToLocalMap& g2l;
120 template<
class GFS,
class M>
122 :
public Dune::CommDataHandleIF<MatrixExchangeDataHandle<GFS,M>,
123 std::pair<typename GFS::Traits::GridView::Grid::Traits::GlobalIdSet::IdType,
124 typename M::block_type> >
128 typedef typename GFS::Traits::GridView
GV;
131 typedef typename GV::Grid
Grid;
133 typedef typename GlobalIdSet::IdType
IdType;
136 typedef typename M::block_type
B;
161 template<
class EntityType>
164 IndexType myindex = indexSet.index(e);
165 typename M::row_type myrow = m[myindex];
166 typename M::row_type::iterator endit=myrow.end();
168 for (
typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
170 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
181 template<
class MessageBuffer,
class EntityType>
182 void gather (MessageBuffer& buff,
const EntityType&
e)
const
184 IndexType myindex = indexSet.index(e);
186 typename M::row_type myrow = m[myindex];
187 typename M::row_type::iterator endit=myrow.end();
189 for (
typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
191 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
194 buff.write(std::make_pair(find->second,*it));
206 template<
class MessageBuffer,
class EntityType>
207 void scatter (MessageBuffer& buff,
const EntityType&
e,
size_t n)
209 IndexType myindex = indexSet.index(e);
210 std::cout << gv.comm().rank() <<
": begin scatter local=" << myindex <<
" size=" << n << std::endl;
213 for (
size_t i=0; i<n; ++i)
216 std::cout << gv.comm().rank() <<
": --> received global=" << x.first << std::endl;
217 typename GlobalToLocalMap::const_iterator find=g2l.find(x.first);
220 IndexType nbindex=find->second;
221 if (m.exists(myindex,nbindex))
223 m[myindex][nbindex] = x.second;
225 block -= m2[myindex][nbindex];
226 std::cout << gv.comm().rank() <<
": compare i=" << myindex <<
" j=" << nbindex
227 <<
" norm=" << block.infinity_norm() << std::endl;
239 : gfs(gfs_), m(m_), gv(gfs.gridView()), indexSet(gv.indexSet()),
240 grid(gv.grid()), globalIdSet(grid.globalIdSet()),
241 l2g(l2g_), g2l(g2l_), m2(m2_)
249 const IndexSet& indexSet;
251 const GlobalIdSet& globalIdSet;
252 const LocalToGlobalMap& l2g;
253 const GlobalToLocalMap& g2l;
260 template<
class GFS,
class T,
class A,
int n,
int m>
262 Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix2)
264 typedef Dune::FieldMatrix<T,n,m> B;
265 typedef Dune::BCRSMatrix<B,A> M;
270 LocalToGlobalMap l2g;
271 GlobalToLocalMap g2l;
273 if (gfs.gridView().comm().size()>1)
274 gfs.gridView().communicate(lgdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
278 if (gfs.gridView().comm().size()>1)
279 gfs.gridView().communicate(mexdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
295 template<
class DGGFS,
class DGMatrix,
class DGPrec,
class DGCC,
296 class CGGFS,
class CGPrec,
class CGCC,
297 class P,
class DGHelper,
class Comm>
299 :
public Dune::Preconditioner<Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>,
300 Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>>
310 const DGHelper& dghelper;
317 typedef typename V::BaseT
X;
318 typedef typename W::BaseT
Y;
335 OvlpDGAMGPrec (
const DGGFS& dggfs_, DGMatrix& dgmatrix_, DGPrec& dgprec_,
const DGCC& dgcc_,
336 const CGGFS& cggfs_, CGPrec& cgprec_,
const CGCC& cgcc_, P& p_,
337 const DGHelper& dghelper_,
const Comm& comm_,
int n1_,
int n2_)
338 : dggfs(dggfs_), dgmatrix(dgmatrix_), dgprec(dgprec_), dgcc(dgcc_),
339 cggfs(cggfs_), cgprec(cgprec_), cgcc(cgcc_), p(p_), dghelper(dghelper_),
340 comm(comm_), n1(n1_), n2(n2_)
372 for (
int i=0; i<n1; i++)
378 if (dggfs.gridView().comm().size()>1)
379 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
386 dghelper.maskForeignDOFs(d);
390 if (cggfs.gridView().comm().size()>1)
391 cggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
393 comm.project(
native(cgd));
407 for (
int i=0; i<n2; i++)
412 if (dggfs.gridView().comm().size()>1)
413 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
449 template<
class DGGO,
class DGCC,
class CGGFS,
class CGCC,
class TransferLOP,
450 template<
class,
class,
class,
int>
class DGPrec,
template<
class>
class Solver,
int s=96>
457 typedef typename DGGO::Traits::TrialGridFunctionSpace
GFS;
460 typedef typename DGGO::Traits::Jacobian
M;
461 typedef typename DGGO::Traits::Domain
V;
476 typedef typename PMatrix::BaseT
P;
479 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type
PTADG;
480 typedef typename Dune::MatMultMatResult<PTADG,P>::type
CGMatrix;
492 std::size_t low_order_space_entries_per_row;
518 unsigned maxiter_=5000,
int verbose_=1,
bool usesuperlu_=
true)
520 , gfs(dggo_.trialGridFunctionSpace())
527 , usesuperlu(usesuperlu_)
528 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
530 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
534 if (usesuperlu ==
true)
536 if (gfs.gridView().comm().rank()==0)
537 std::cout <<
"WARNING: You are using AMG without SuperLU!"
538 <<
" Please consider installing SuperLU,"
539 <<
" or set the usesuperlu flag to false"
540 <<
" to suppress this warning." << std::endl;
546 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
555 const ParameterTree& params)
557 , gfs(dggo_.trialGridFunctionSpace())
562 , maxiter(params.get<int>(
"max_iterations",5000))
563 , verbose(params.get<int>(
"verbose",1))
564 , usesuperlu(params.get<bool>(
"use_superlu",true))
565 , low_order_space_entries_per_row(params.get<
std::size_t>(
"low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
567 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
571 if (usesuperlu ==
true)
573 if (gfs.gridView().comm().rank()==0)
574 std::cout <<
"WARNING: You are using AMG without SuperLU!"
575 <<
" Please consider installing SuperLU,"
576 <<
" or set the usesuperlu flag to false"
577 <<
" to suppress this warning." << std::endl;
583 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout <<
"allocated prolongation matrix of size " << pmatrix.N() <<
" x " << pmatrix.M() << std::endl;
596 void apply (M& A, V& z, V& r,
typename V::ElementType reduction)
609 CGGO cggo(cggfs,cgcc,cggfs,cgcc,emptylop,
MBE(low_order_space_entries_per_row));
610 typedef typename CGGO::Jacobian CGM;
616 CGM acg(attached_container);
619 Dune::transposeMatMultMat(ptadg,
native(pmatrix),
native(A));
623 double triple_product_time = watch.elapsed();
624 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout <<
"=== triple matrix product " << triple_product_time <<
" s" << std::endl;
627 cggo.jacobian(cgx,acg);
632 DGGOEmpty dggoempty(gfs,dgcc,gfs,dgcc,emptylop,
MBE(1 << GFS::Traits::GridView::dimension));
633 dggoempty.jacobian(z,A);
640 Comm oocc(gfs.gridView().comm());
642 CGHELPER cghelper(cggfs,2);
643 cghelper.createIndexSetAndProjectForAMG(acg,oocc);
644 typedef Dune::OverlappingSchwarzOperator<CGMatrix,CGVector,CGVector,Comm> ParCGOperator;
645 ParCGOperator paroop(
native(acg),oocc);
646 Dune::OverlappingSchwarzScalarProduct<CGVector,Comm> sp(oocc);
647 typedef Dune::Amg::Parameters Parameters;
648 Parameters params(15,2000);
649 params.setDefaultValuesIsotropic(CGGFS::Traits::GridViewType::Traits::Grid::dimension);
650 params.setDebugLevel(verbose);
651 params.setCoarsenTarget(2000);
652 params.setMaxLevel(20);
653 params.setProlongationDampingFactor(1.6);
654 params.setNoPreSmoothSteps(3);
655 params.setNoPostSmoothSteps(3);
657 params.setAdditive(
false);
659 typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
660 typedef Dune::BlockPreconditioner<CGVector,CGVector,Comm,Smoother> ParSmoother;
661 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
662 SmootherArgs smootherArgs;
663 smootherArgs.iterations = 2;
664 smootherArgs.relaxationFactor = 0.92;
665 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
666 Criterion criterion(params);
667 typedef Dune::Amg::AMG<ParCGOperator,CGVector,ParSmoother,Comm> AMG;
669 AMG amg(paroop,criterion,smootherArgs,oocc);
670 double amg_setup_time = watch.elapsed();
671 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout <<
"=== AMG setup " <<amg_setup_time <<
" s" << std::endl;
674 typedef DGPrec<Matrix,Vector,Vector,1> DGPrecType;
675 DGPrecType dgprec(
native(A),1,0.92);
679 HybridPrec hybridprec(gfs,
native(A),dgprec,dgcc,cggfs,amg,cgcc,
native(pmatrix),
694 if (gfs.gridView().comm().rank()>0) verb=0;
695 Solver<V> solver(pop,psp,hybridprec,reduction,maxiter,verb);
698 Dune::InverseOperatorResult stat;
700 solver.apply(z,r,stat);
701 double amg_solve_time = watch.elapsed();
702 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout <<
"=== Hybrid total solve time " << amg_solve_time+amg_setup_time+triple_product_time <<
" s" << std::endl;
705 res.
elapsed = amg_solve_time+amg_setup_time+triple_product_time;
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:146
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: ovlp_amg_dg_backend.hh:152
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:38
void jacobian(const Domain &x, Jacobian &a) const
Assembler jacobian.
Definition: gridoperator.hh:172
Definition: istl/ovlpistlsolverbackend.hh:323
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:88
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:182
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:128
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::range_type::field_type > CGW
Definition: ovlp_amg_dg_backend.hh:320
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:803
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:130
DGGO::Traits::TrialGridFunctionSpace GFS
Definition: ovlp_amg_dg_backend.hh:457
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:48
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:64
virtual void post(V &x)
Clean up.
Definition: ovlp_amg_dg_backend.hh:425
Vector::field_type field_type
Definition: ovlp_amg_dg_backend.hh:464
MatrixExchangeDataHandle(const GFS &gfs_, M &m_, const LocalToGlobalMap &l2g_, const GlobalToLocalMap &g2l_, M &m2_)
constructor
Definition: ovlp_amg_dg_backend.hh:238
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:34
OvlpDGAMGPrec(const DGGFS &dggfs_, DGMatrix &dgmatrix_, DGPrec &dgprec_, const DGCC &dgcc_, const CGGFS &cggfs_, CGPrec &cgprec_, const CGCC &cgcc_, P &p_, const DGHelper &dghelper_, const Comm &comm_, int n1_, int n2_)
Constructor.
Definition: ovlp_amg_dg_backend.hh:335
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:112
const E & e
Definition: interpolate.hh:172
Implement jacobian_apply_volume() based on alpha_volume()
Definition: defaultimp.hh:321
RFType reduction
Definition: solver.hh:35
PGO::Jacobian PMatrix
Definition: ovlp_amg_dg_backend.hh:475
PMatrix::BaseT P
Definition: ovlp_amg_dg_backend.hh:476
Definition: parallelhelper.hh:45
void restore_overlap_entries(const GFS &gfs, Dune::BCRSMatrix< Dune::FieldMatrix< T, n, m >, A > &matrix, Dune::BCRSMatrix< Dune::FieldMatrix< T, n, m >, A > &matrix2)
Definition: ovlp_amg_dg_backend.hh:261
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:71
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, const ParameterTree ¶ms)
Definition: ovlp_amg_dg_backend.hh:554
virtual void pre(V &x, W &b)
Prepare the preconditioner.
Definition: ovlp_amg_dg_backend.hh:349
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:129
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:45
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:133
unsigned int iterations
Definition: solver.hh:33
Backend::attached_container attached_container
Definition: backend/common/tags.hh:37
The category the preconditioner is part of.
Definition: ovlp_amg_dg_backend.hh:325
Default flags for all local operators.
Definition: flags.hh:18
sparsity pattern generator
Definition: pattern.hh:13
int DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:41
Definition: istl/ovlpistlsolverbackend.hh:371
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, unsigned maxiter_=5000, int verbose_=1, bool usesuperlu_=true)
Definition: ovlp_amg_dg_backend.hh:517
std::enable_if< std::is_base_of< impl::WrapperBase, T >::value, Native< T > & >::type native(T &t)
Definition: backend/interface.hh:198
TransferLOP CGTODGLOP
Definition: ovlp_amg_dg_backend.hh:473
bool converged
Definition: solver.hh:32
Tag for requesting a vector or matrix container with a pre-attached underlying object.
Definition: backend/common/tags.hh:27
CGV::BaseT CGVector
Definition: ovlp_amg_dg_backend.hh:468
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:36
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:33
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:162
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
Dune::MatMultMatResult< PTADG, P >::type CGMatrix
Definition: ovlp_amg_dg_backend.hh:480
Definition: genericdatahandle.hh:622
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:37
LocalGlobalMapDataHandle(const GFS &gfs_, LocalToGlobalMap &l2g_, GlobalToLocalMap &g2l_)
constructor
Definition: ovlp_amg_dg_backend.hh:101
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:132
W::BaseT Y
Definition: ovlp_amg_dg_backend.hh:318
Definition: adaptivity.hh:27
V::BaseT X
Definition: ovlp_amg_dg_backend.hh:317
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:421
Definition: constraintstransformation.hh:111
Dune::PDELab::GridOperator< CGGFS, GFS, CGTODGLOP, MBE, field_type, field_type, field_type, CC, CC > PGO
Definition: ovlp_amg_dg_backend.hh:474
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:131
DGGO::Traits::Jacobian M
Definition: ovlp_amg_dg_backend.hh:460
const istl::ParallelHelper< DGGO::Traits::TrialGridFunctionSpace > & parallelHelper() const
Definition: istl/ovlpistlsolverbackend.hh:353
M::block_type B
Definition: ovlp_amg_dg_backend.hh:136
bool fixedsize(int dim, int codim) const
returns true if size per entity of given dim and codim is a constant
Definition: ovlp_amg_dg_backend.hh:54
Backend using (possibly nested) ISTL BCRSMatrices.
Definition: bcrsmatrixbackend.hh:187
static const int dim
Definition: adaptivity.hh:83
Dune::TransposedMatMultMatResult< P, Matrix >::type PTADG
Definition: ovlp_amg_dg_backend.hh:479
Dune::PDELab::Backend::Vector< CGGFS, field_type > CGV
Definition: ovlp_amg_dg_backend.hh:467
Definition: istl/ovlpistlsolverbackend.hh:41
M::BaseT Matrix
Definition: ovlp_amg_dg_backend.hh:462
const std::string s
Definition: function.hh:1102
Definition: ovlp_amg_dg_backend.hh:28
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:207
Definition: ovlp_amg_dg_backend.hh:298
V::BaseT Vector
Definition: ovlp_amg_dg_backend.hh:463
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:52
PMatrix & prolongation_matrix()
Definition: ovlp_amg_dg_backend.hh:510
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::domain_type::field_type > CGV
Definition: ovlp_amg_dg_backend.hh:319
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::domain_type::field_type > V
Definition: ovlp_amg_dg_backend.hh:315
std::pair< IdType, B > DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:139
Dune::PDELab::istl::BCRSMatrixBackend MBE
Definition: ovlp_amg_dg_backend.hh:471
Definition: ovlp_amg_dg_backend.hh:121
Dune::PDELab::EmptyTransformation CC
Definition: ovlp_amg_dg_backend.hh:472
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:142
DGGO::Traits::Domain V
Definition: ovlp_amg_dg_backend.hh:461
virtual void apply(V &x, const W &b)
Apply the precondioner.
Definition: ovlp_amg_dg_backend.hh:363
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::range_type::field_type > W
Definition: ovlp_amg_dg_backend.hh:316
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:44
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:35
Definition: ovlp_amg_dg_backend.hh:451
Default class for additional methods in instationary local operators.
Definition: idefault.hh:89
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:143
void apply(M &A, V &z, V &r, typename V::ElementType reduction)
solve the given linear system
Definition: ovlp_amg_dg_backend.hh:596