dune-pdelab 2.7-git
Loading...
Searching...
No Matches
ovlp_amg_dg_backend.hh
Go to the documentation of this file.
1#ifndef DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
2#define DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
3
4#include <dune/common/parametertree.hh>
5#include <dune/common/power.hh>
6
7#include <dune/istl/matrixmatrix.hh>
8
9#include <dune/grid/common/datahandleif.hh>
19
20namespace Dune {
21 namespace PDELab {
22 //***********************************************************
23 // A data handle / function to communicate matrix entries
24 // in the overlap. It turned out that it is actually not
25 // required to do that, but anyway it is there now.
26 //***********************************************************
27
30 template<class GFS>
32 : public Dune::CommDataHandleIF<LocalGlobalMapDataHandle<GFS>,int>
33 {
34 public:
35 // some types from the grid view
36 typedef typename GFS::Traits::GridView GV;
37 typedef typename GV::IndexSet IndexSet;
38 typedef typename IndexSet::IndexType IndexType;
39 typedef typename GV::Grid Grid;
40 typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
41 typedef typename GlobalIdSet::IdType IdType;
42
44 typedef int DataType;
45
46 // map types
47 typedef std::map<IndexType,IdType> LocalToGlobalMap;
48 typedef std::map<IdType,IndexType> GlobalToLocalMap;
49
51 bool contains (int dim, int codim) const
52 {
53 return (codim==0);
54 }
55
57 bool fixedSize (int dim, int codim) const
58 {
59 return true;
60 }
61
66 template<class EntityType>
67 size_t size (EntityType& e) const
68 {
69 return 1;
70 }
71
73 template<class MessageBuffer, class EntityType>
74 void gather (MessageBuffer& buff, const EntityType& e) const
75 {
76 // fill map
77 IndexType myindex = indexSet.index(e);
78 IdType myid = globalIdSet.id(e);
79 l2g[myindex] = myid;
80 g2l[myid] = myindex;
81 //std::cout << gv.comm().rank() << ": gather local=" << myindex << " global=" << myid << std::endl;
82
83 buff.write(0); // this is a dummy, we are not interested in the data
84 }
85
90 template<class MessageBuffer, class EntityType>
91 void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
92 {
93 DataType x;
94 buff.read(x); // this is a dummy, we are not interested in the data
95
96 IndexType myindex = indexSet.index(e);
97 IdType myid = globalIdSet.id(e);
98 l2g[myindex] = myid;
99 g2l[myid] = myindex;
100 //std::cout << gv.comm().rank() << ": scatter local=" << myindex << " global=" << myid << std::endl;
101 }
102
105 : gfs(gfs_), gv(gfs.gridView()), indexSet(gv.indexSet()),
106 grid(gv.grid()), globalIdSet(grid.globalIdSet()),
107 l2g(l2g_), g2l(g2l_)
108 {
109 }
110
111 private:
112 const GFS& gfs;
113 GV gv;
114 const IndexSet& indexSet;
115 const Grid& grid;
116 const GlobalIdSet& globalIdSet;
117 LocalToGlobalMap& l2g;
118 GlobalToLocalMap& g2l;
119 };
120
121
122 // A DataHandle class to exchange rows of a sparse matrix
123 template<class GFS, class M> // mapper type and vector type
125 : public Dune::CommDataHandleIF<MatrixExchangeDataHandle<GFS,M>,
126 std::pair<typename GFS::Traits::GridView::Grid::Traits::GlobalIdSet::IdType,
127 typename M::block_type> >
128 {
129 public:
130 // some types from the grid view
131 typedef typename GFS::Traits::GridView GV;
132 typedef typename GV::IndexSet IndexSet;
133 typedef typename IndexSet::IndexType IndexType;
134 typedef typename GV::Grid Grid;
135 typedef typename Grid::Traits::GlobalIdSet GlobalIdSet;
136 typedef typename GlobalIdSet::IdType IdType;
137
138 // some types from the matrix
139 typedef typename M::block_type B;
140
142 typedef std::pair<IdType,B> DataType;
143
144 // map types
145 typedef std::map<IndexType,IdType> LocalToGlobalMap;
146 typedef std::map<IdType,IndexType> GlobalToLocalMap;
147
149 bool contains (int dim, int codim) const
150 {
151 return (codim==0);
152 }
153
155 bool fixedSize (int dim, int codim) const
156 {
157 return false;
158 }
159
164 template<class EntityType>
165 size_t size (EntityType& e) const
166 {
167 IndexType myindex = indexSet.index(e);
168 typename M::row_type myrow = m[myindex];
169 typename M::row_type::iterator endit=myrow.end();
170 size_t count=0;
171 for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
172 {
173 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
174 if (find!=l2g.end())
175 {
176 count++;
177 }
178 }
179 //std::cout << gv.comm().rank() << ": row=" << myindex << " size=" << count << std::endl;
180 return count;
181 }
182
184 template<class MessageBuffer, class EntityType>
185 void gather (MessageBuffer& buff, const EntityType& e) const
186 {
187 IndexType myindex = indexSet.index(e);
188 //std::cout << gv.comm().rank() << ": begin gather local=" << myindex << std::endl;
189 typename M::row_type myrow = m[myindex];
190 typename M::row_type::iterator endit=myrow.end();
191 size_t count = 0;
192 for (typename M::row_type::iterator it=myrow.begin(); it!=endit; ++it)
193 {
194 typename LocalToGlobalMap::const_iterator find=l2g.find(it.index());
195 if (find!=l2g.end())
196 {
197 buff.write(std::make_pair(find->second,*it));
198 //std::cout << gv.comm().rank() << ": ==> local=" << find->first << " global=" << find->second << std::endl;
199 count++;
200 }
201 }
202 //std::cout << gv.comm().rank() << ": end gather row=" << myindex << " size=" << count << std::endl;
203 }
204
209 template<class MessageBuffer, class EntityType>
210 void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
211 {
212 IndexType myindex = indexSet.index(e);
213 std::cout << gv.comm().rank() << ": begin scatter local=" << myindex << " size=" << n << std::endl;
214 DataType x;
215 size_t count = 0;
216 for (size_t i=0; i<n; ++i)
217 {
218 buff.read(x);
219 std::cout << gv.comm().rank() << ": --> received global=" << x.first << std::endl;
220 typename GlobalToLocalMap::const_iterator find=g2l.find(x.first);
221 if (find!=g2l.end())
222 {
223 IndexType nbindex=find->second;
224 if (m.exists(myindex,nbindex))
225 {
226 m[myindex][nbindex] = x.second;
227 B block(x.second);
228 block -= m2[myindex][nbindex];
229 std::cout << gv.comm().rank() << ": compare i=" << myindex << " j=" << nbindex
230 << " norm=" << block.infinity_norm() << std::endl;
231
232 count++;
233 //std::cout << gv.comm().rank() << ": --> local=" << find->first << " global=" << find->second << std::endl;
234 }
235 }
236 }
237 //std::cout << gv.comm().rank() << ": end scatter row=" << myindex << " size=" << count << std::endl;
238 }
239
241 MatrixExchangeDataHandle (const GFS& gfs_, M& m_, const LocalToGlobalMap& l2g_, const GlobalToLocalMap& g2l_,M& m2_)
242 : gfs(gfs_), m(m_), gv(gfs.gridView()), indexSet(gv.indexSet()),
243 grid(gv.grid()), globalIdSet(grid.globalIdSet()),
244 l2g(l2g_), g2l(g2l_), m2(m2_)
245 {
246 }
247
248 private:
249 const GFS& gfs;
250 M& m;
251 GV gv;
252 const IndexSet& indexSet;
253 const Grid& grid;
254 const GlobalIdSet& globalIdSet;
255 const LocalToGlobalMap& l2g;
256 const GlobalToLocalMap& g2l;
257 M& m2;
258 };
259
260
263 template<class GFS, class T, class A, int n, int m>
264 void restore_overlap_entries (const GFS& gfs, Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix,
265 Dune::BCRSMatrix<Dune::FieldMatrix<T,n,m>,A>& matrix2)
266 {
267 typedef Dune::FieldMatrix<T,n,m> B;
268 typedef Dune::BCRSMatrix<B,A> M; // m is of type M
269 typedef typename LocalGlobalMapDataHandle<GFS>::LocalToGlobalMap LocalToGlobalMap;
270 typedef typename LocalGlobalMapDataHandle<GFS>::GlobalToLocalMap GlobalToLocalMap;
271
272 // build up two maps to associate local indices and global ids
273 LocalToGlobalMap l2g;
274 GlobalToLocalMap g2l;
275 LocalGlobalMapDataHandle<GFS> lgdh(gfs,l2g,g2l);
276 if (gfs.gridView().comm().size()>1)
277 gfs.gridView().communicate(lgdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
278
279 // exchange matrix entries
280 MatrixExchangeDataHandle<GFS,M> mexdh(gfs,matrix,l2g,g2l,matrix2);
281 if (gfs.gridView().comm().size()>1)
282 gfs.gridView().communicate(mexdh,Dune::InteriorBorder_All_Interface,Dune::ForwardCommunication);
283 }
284
285 //***********************************************************
286 // The DG/AMG preconditioner in the overlapping case
287 //***********************************************************
288
298 template<class DGGFS, class DGMatrix, class DGPrec, class DGCC,
299 class CGGFS, class CGPrec, class CGCC,
300 class P, class DGHelper, class Comm>
302 : public Dune::Preconditioner<Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::domain_type::field_type>,
303 Dune::PDELab::Backend::Vector<DGGFS,typename DGPrec::range_type::field_type>>
304 {
305 const DGGFS& dggfs;
306 DGMatrix& dgmatrix;
307 DGPrec& dgprec;
308 const DGCC& dgcc;
309 const CGGFS& cggfs;
310 CGPrec& cgprec;
311 const CGCC& cgcc;
312 P& p;
313 const DGHelper& dghelper;
314 const Comm& comm;
315 int n1,n2;
316
317 public:
324
325 // define the category
326 SolverCategory::Category category() const override
327 {
328 return SolverCategory::overlapping;
329 }
330
338 OvlpDGAMGPrec (const DGGFS& dggfs_, DGMatrix& dgmatrix_, DGPrec& dgprec_, const DGCC& dgcc_,
339 const CGGFS& cggfs_, CGPrec& cgprec_, const CGCC& cgcc_, P& p_,
340 const DGHelper& dghelper_, const Comm& comm_, int n1_, int n2_)
341 : dggfs(dggfs_), dgmatrix(dgmatrix_), dgprec(dgprec_), dgcc(dgcc_),
342 cggfs(cggfs_), cgprec(cgprec_), cgcc(cgcc_), p(p_), dghelper(dghelper_),
343 comm(comm_), n1(n1_), n2(n2_)
344 {
345 }
346
352 virtual void pre (V& x, W& b) override
353 {
354 using Backend::native;
355 dgprec.pre(native(x),native(b));
356 CGW cgd(cggfs,0.0);
357 CGV cgv(cggfs,0.0);
358 cgprec.pre(native(cgv),native(cgd));
359 }
360
366 virtual void apply (V& x, const W& b) override
367 {
368 using Backend::native;
369 // need local copies to store defect and solution
370 W d(b);
372 V v(x); // only to get correct size
373
374 // pre-smoothing on DG matrix
375 for (int i=0; i<n1; i++)
376 {
377 using Backend::native;
378 v = 0.0;
379 dgprec.apply(native(v),native(d));
381 if (dggfs.gridView().comm().size()>1)
382 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
383 dgmatrix.mmv(native(v),native(d));
385 x += v;
386 }
387
388 // restrict defect to CG subspace
389 dghelper.maskForeignDOFs(d); // DG defect is additive for overlap 1, but in case we use more
390 CGW cgd(cggfs,0.0);
391 p.mtv(native(d),native(cgd));
393 if (cggfs.gridView().comm().size()>1)
394 cggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication); // now we have consistent defect on coarse grid
396 comm.project(native(cgd));
397 CGV cgv(cggfs,0.0);
398
399
400 // call preconditioner
401 cgprec.apply(native(cgv),native(cgd));
402
403 // prolongate correction
404 p.mv(native(cgv),native(v));
405 dgmatrix.mmv(native(v),native(d));
407 x += v;
408
409 // post-smoothing on DG matrix
410 for (int i=0; i<n2; i++)
411 {
412 v = 0.0;
413 dgprec.apply(native(v),native(d));
415 if (dggfs.gridView().comm().size()>1)
416 dggfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
417 dgmatrix.mmv(native(v),native(d));
419 x += v;
420 }
421 }
422
428 virtual void post (V& x) override
429 {
430 dgprec.post(Backend::native(x));
431 CGV cgv(cggfs,0.0);
432 cgprec.post(Backend::native(cgv));
433 }
434 };
435
436//***********************************************************
437// The DG/AMG linear solver backend in the overlapping case
438//***********************************************************
439
463template<class DGGO, class DGCC, class CGGFS, class CGCC, class TransferLOP,
464 template<class,class,class,int> class DGPrec, template<class> class Solver, int s=96>
466 public Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>,
468{
469public:
470 // DG grid function space
471 typedef typename DGGO::Traits::TrialGridFunctionSpace GFS;
472
473 // vectors and matrices on DG level
474 typedef typename DGGO::Traits::Jacobian M; // wrapped istl DG matrix
475 typedef typename DGGO::Traits::Domain V; // wrapped istl DG vector
476 typedef Backend::Native<M> Matrix; // istl DG matrix
477 typedef Backend::Native<V> Vector; // istl DG vector
478 typedef typename Vector::field_type field_type;
479
480 // vectors and matrices on CG level
481 using CGV = Dune::PDELab::Backend::Vector<CGGFS,field_type>; // wrapped istl CG vector
482 typedef Backend::Native<CGV> CGVector; // istl CG vector
483
484 // prolongation matrix
487 typedef TransferLOP CGTODGLOP; // local operator
489 typedef typename PGO::Jacobian PMatrix; // wrapped ISTL prolongation matrix
490 typedef Backend::Native<PMatrix> P; // ISTL prolongation matrix
491
492 // CG subspace matrix
493 typedef typename Dune::TransposedMatMultMatResult<P,Matrix>::type PTADG;
494 typedef typename Dune::MatMultMatResult<PTADG,P>::type CGMatrix; // istl coarse space matrix
495
496 // AMG in CG-subspace
498 typedef Dune::OverlappingSchwarzOperator<CGMatrix,CGVector,CGVector,Comm> ParCGOperator;
499 typedef Dune::SeqSSOR<CGMatrix,CGVector,CGVector,1> Smoother;
500 typedef Dune::BlockPreconditioner<CGVector,CGVector,Comm,Smoother> ParSmoother;
501 typedef Dune::Amg::AMG<ParCGOperator,CGVector,ParSmoother,Comm> AMG;
502 typedef Dune::Amg::Parameters Parameters;
503
504private:
505
506 const GFS& gfs;
507 DGGO& dggo;
508 const DGCC& dgcc;
509 CGGFS& cggfs;
510 const CGCC& cgcc;
511 std::shared_ptr<AMG> amg;
512 Parameters amg_parameters;
513 unsigned maxiter;
514 int verbose;
515 bool reuse;
516 bool firstapply;
517 bool usesuperlu;
518 std::size_t low_order_space_entries_per_row;
519
520 // Parameters for DG smoother
521 double smoother_relaxation = 0.92; // Relaxation parameter of DG smoother
522 int n1=3; // Number of DG pre-smoothing steps
523 int n2=3; // Number of DG post-smoothing steps
524
525 CGTODGLOP cgtodglop; // local operator to assemble prolongation matrix
526 PGO pgo; // grid operator to assemble prolongation matrix
527 PMatrix pmatrix; // wrapped prolongation matrix
528
531 class EmptyLop : public Dune::PDELab::NumericalJacobianApplyVolume<EmptyLop>,
535 {
536 };
537
538 // grid operator with empty local operator for constraints assembly in parallel
540 // CG-subspace matrix
541 typedef typename CGGO::Jacobian CGM;
542 CGM acg;
543
544public:
545
546 // access to prolongation matrix
548 {
549 return pmatrix;
550 }
551
556 void setParameters(const Parameters& amg_parameters_)
557 {
558 amg_parameters = amg_parameters_;
559 }
560
568 const Parameters& parameters() const
569 {
570 return amg_parameters;
571 }
572
574 void setReuse(bool reuse_)
575 {
576 reuse = reuse_;
577 }
578
580 bool getReuse() const
581 {
582 return reuse;
583 }
584
587 ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
588 unsigned maxiter_=5000, int verbose_=1, bool reuse_=false,
589 bool usesuperlu_=true)
590 : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
591 , gfs(dggo_.trialGridFunctionSpace())
592 , dggo(dggo_)
593 , dgcc(dgcc_)
594 , cggfs(cggfs_)
595 , cgcc(cgcc_)
596 , amg_parameters(15,2000)
597 , maxiter(maxiter_)
598 , verbose(verbose_)
599 , reuse(reuse_)
600 , firstapply(true)
601 , usesuperlu(usesuperlu_)
602 , low_order_space_entries_per_row(StaticPower<3,GFS::Traits::GridView::dimension>::power)
603 , cgtodglop()
604 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
605 , pmatrix(pgo)
606 , acg(Backend::attached_container())
607 {
608 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
609 amg_parameters.setDebugLevel(verbose_);
610#if !HAVE_SUPERLU
611 if (usesuperlu == true)
612 {
613 if (gfs.gridView().comm().rank()==0)
614 std::cout << "WARNING: You are using AMG without SuperLU!"
615 << " Please consider installing SuperLU,"
616 << " or set the usesuperlu flag to false"
617 << " to suppress this warning." << std::endl;
618 }
619#endif
620
621 // assemble prolongation matrix; this will not change from one apply to the next
622 pmatrix = 0.0;
623 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
624 CGV cgx(cggfs,0.0); // need vector to call jacobian
625 pgo.jacobian(cgx,pmatrix);
626 }
627
628
631 ISTLBackend_OVLP_AMG_4_DG(DGGO& dggo_, const DGCC& dgcc_, CGGFS& cggfs_, const CGCC& cgcc_,
632 const ParameterTree& params)
633 : Dune::PDELab::OVLPScalarProductImplementation<typename DGGO::Traits::TrialGridFunctionSpace>(dggo_.trialGridFunctionSpace())
634 , gfs(dggo_.trialGridFunctionSpace())
635 , dggo(dggo_)
636 , dgcc(dgcc_)
637 , cggfs(cggfs_)
638 , cgcc(cgcc_)
639 , maxiter(params.get<int>("max_iterations",5000))
640 , amg_parameters(15,2000)
641 , verbose(params.get<int>("verbose",1))
642 , reuse(params.get<bool>("reuse",false))
643 , firstapply(true)
644 , usesuperlu(params.get<bool>("use_superlu",true))
645 , low_order_space_entries_per_row(params.get<std::size_t>("low_order_space.entries_per_row",StaticPower<3,GFS::Traits::GridView::dimension>::power))
646 , cgtodglop()
647 , pgo(cggfs,dggo.trialGridFunctionSpace(),cgtodglop,MBE(low_order_space_entries_per_row))
648 , pmatrix(pgo)
649 , acg(Backend::attached_container())
650 {
651 amg_parameters.setDefaultValuesIsotropic(GFS::Traits::GridViewType::Traits::Grid::dimension);
652 amg_parameters.setDebugLevel(params.get<int>("verbose",1));
653#if !HAVE_SUPERLU
654 if (usesuperlu == true)
655 {
656 if (gfs.gridView().comm().rank()==0)
657 std::cout << "WARNING: You are using AMG without SuperLU!"
658 << " Please consider installing SuperLU,"
659 << " or set the usesuperlu flag to false"
660 << " to suppress this warning." << std::endl;
661 }
662#endif
663
664 // assemble prolongation matrix; this will not change from one apply to the next
665 pmatrix = 0.0;
666 if (verbose>0 && gfs.gridView().comm().rank()==0) std::cout << "allocated prolongation matrix of size " << pmatrix.N() << " x " << pmatrix.M() << std::endl;
667 CGV cgx(cggfs,0.0); // need vector to call jacobian
668 pgo.jacobian(cgx,pmatrix);
669 }
670
671
673 void setDGSmootherRelaxation(double relaxation_)
674 {
675 smoother_relaxation = relaxation_;
676 }
677
680 {
681 n1 = n1_;
682 }
683
686 {
687 n2 = n2_;
688 }
689
690
698 void apply (M& A, V& z, V& r, typename Dune::template FieldTraits<typename V::ElementType >::real_type reduction)
699 {
700 using Backend::native;
701 // make operator and scalar product for overlapping solver
703 POP pop(dgcc,A);
705 PSP psp(*this);
706
707 // compute CG matrix
708 // make grid operator with empty local operator => matrix data type and constraints assembly
709 EmptyLop emptylop;
710 CGGO cggo(cggfs,cgcc,cggfs,cgcc,emptylop,MBE(low_order_space_entries_per_row));
711
712 // do triple matrix product ACG = P^T ADG P; this is purely local
713 Dune::Timer watch;
714 watch.reset();
715 // only do triple matrix product if the matrix changes
716 double triple_product_time = 0.0;
717 // no need to set acg here back to zero, this is done in matMultmat
718 if(reuse == false || firstapply == true) {
719 PTADG ptadg;
720 Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A)); // 1a
721 //Dune::transposeMatMultMat(ptadg,native(pmatrix),native(A2)); // 1b
722 Dune::matMultMat(native(acg),ptadg,native(pmatrix));
723 triple_product_time = watch.elapsed();
724 if (verbose>0 && gfs.gridView().comm().rank()==0)
725 std::cout << "=== triple matrix product " << triple_product_time << " s" << std::endl;
726 //Dune::printmatrix(std::cout,native(acg),"triple product matrix","row",10,2);
727 CGV cgx(cggfs,0.0); // need vector to call jacobian
728 cggo.jacobian(cgx,acg); // insert trivial rows at processor boundaries
729 //std::cout << "CG constraints: " << cgcc.size() << " out of " << cggfs.globalSize() << std::endl;
730 }
731 else if(verbose>0 && gfs.gridView().comm().rank()==0)
732 std::cout << "=== reuse CG matrix, SKIPPING triple matrix product " << std::endl;
733
734 // NOW we need to insert the processor boundary conditions in DG matrix
736 DGGOEmpty dggoempty(gfs,dgcc,gfs,dgcc,emptylop,MBE(1 << GFS::Traits::GridView::dimension));
737 dggoempty.jacobian(z,A);
738
739 // and in the residual
741
742 // now set up parallel AMG solver for the CG subspace
743 Comm oocc(gfs.gridView().comm());
745 CGHELPER cghelper(cggfs,2);
746 cghelper.createIndexSetAndProjectForAMG(acg,oocc);
747 ParCGOperator paroop(native(acg),oocc);
748 Dune::OverlappingSchwarzScalarProduct<CGVector,Comm> sp(oocc);
749
750 typedef typename Dune::Amg::SmootherTraits<ParSmoother>::Arguments SmootherArgs;
751 SmootherArgs smootherArgs;
752 smootherArgs.iterations = 1;
753 smootherArgs.relaxationFactor = 1.0;
754 typedef Dune::Amg::CoarsenCriterion<Dune::Amg::SymmetricCriterion<CGMatrix,Dune::Amg::FirstDiagonal> > Criterion;
755 Criterion criterion(amg_parameters);
756 watch.reset();
757
758 // only construct a new AMG for the CG-subspace if the matrix changes
759 double amg_setup_time = 0.0;
760 if(reuse == false || firstapply == true) {
761 amg.reset(new AMG(paroop,criterion,smootherArgs,oocc));
762 firstapply = false;
763 amg_setup_time = watch.elapsed();
764 if (verbose>0 && gfs.gridView().comm().rank()==0)
765 std::cout << "=== AMG setup " <<amg_setup_time << " s" << std::endl;
766 }
767 else if (verbose>0 && gfs.gridView().comm().rank()==0)
768 std::cout << "=== reuse CG matrix, SKIPPING AMG setup " << std::endl;
769
770 // set up hybrid DG/CG preconditioner
771 typedef DGPrec<Matrix,Vector,Vector,1> DGPrecType;
772 DGPrecType dgprec(native(A),1,smoother_relaxation);
775 HybridPrec hybridprec(gfs,native(A),dgprec,dgcc,cggfs,*amg,cgcc,native(pmatrix),
776 this->parallelHelper(),oocc,n1,n2);
777
778 // /********/
779 // /* Test */
780 // /********/
781 // Solver<CGVector> testsolver(paroop,sp,amg,1e-8,100,2);
782 // CGV cgxx(cggfs,0.0);
783 // CGV cgdd(cggfs,1.0);
784 // Dune::InverseOperatorResult statstat;
785 // testsolver.apply(native(cgxx),native(cgdd),statstat);
786 // /********/
787
788 // set up solver
789 int verb=verbose;
790 if (gfs.gridView().comm().rank()>0) verb=0;
791 Solver<V> solver(pop,psp,hybridprec,reduction,maxiter,verb);
792
793 // solve
794 Dune::InverseOperatorResult stat;
795 watch.reset();
796 solver.apply(z,r,stat);
797 double amg_solve_time = watch.elapsed();
798 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;
799 res.converged = stat.converged;
800 res.iterations = stat.iterations;
801 res.elapsed = amg_solve_time+amg_setup_time+triple_product_time;
802 res.reduction = stat.reduction;
803 res.conv_rate = stat.conv_rate;
804 }
805
806};
807}
808}
809#endif // DUNE_PDELAB_BACKEND_ISTL_OVLP_AMG_DG_BACKEND_HH
const std::string s
Definition: function.hh:843
static const int dim
Definition: adaptivity.hh:84
const P & p
Definition: constraints.hh:148
void set_constrained_dofs(const CG &cg, typename XG::ElementType x, XG &xg)
construct constraints from given boundary condition function
Definition: constraints.hh:796
STL namespace.
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
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:264
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: ovlp_amg_dg_backend.hh:33
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:37
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:47
LocalGlobalMapDataHandle(const GFS &gfs_, LocalToGlobalMap &l2g_, GlobalToLocalMap &g2l_)
constructor
Definition: ovlp_amg_dg_backend.hh:104
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:57
GlobalIdSet::IdType IdType
Definition: ovlp_amg_dg_backend.hh:41
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:36
int DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:44
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:48
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:74
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:67
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:51
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:40
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:39
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:38
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:91
Definition: ovlp_amg_dg_backend.hh:128
bool contains(int dim, int codim) const
returns true if data for this codim should be communicated
Definition: ovlp_amg_dg_backend.hh:149
MatrixExchangeDataHandle(const GFS &gfs_, M &m_, const LocalToGlobalMap &l2g_, const GlobalToLocalMap &g2l_, M &m2_)
constructor
Definition: ovlp_amg_dg_backend.hh:241
std::map< IndexType, IdType > LocalToGlobalMap
Definition: ovlp_amg_dg_backend.hh:145
M::block_type B
Definition: ovlp_amg_dg_backend.hh:139
GV::Grid Grid
Definition: ovlp_amg_dg_backend.hh:134
void scatter(MessageBuffer &buff, const EntityType &e, size_t n)
Definition: ovlp_amg_dg_backend.hh:210
void gather(MessageBuffer &buff, const EntityType &e) const
pack data from user to message buffer
Definition: ovlp_amg_dg_backend.hh:185
Grid::Traits::GlobalIdSet GlobalIdSet
Definition: ovlp_amg_dg_backend.hh:135
std::map< IdType, IndexType > GlobalToLocalMap
Definition: ovlp_amg_dg_backend.hh:146
std::pair< IdType, B > DataType
export type of data for message buffer
Definition: ovlp_amg_dg_backend.hh:142
GV::IndexSet IndexSet
Definition: ovlp_amg_dg_backend.hh:132
IndexSet::IndexType IndexType
Definition: ovlp_amg_dg_backend.hh:133
GlobalIdSet::IdType IdType
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:155
GFS::Traits::GridView GV
Definition: ovlp_amg_dg_backend.hh:131
size_t size(EntityType &e) const
Definition: ovlp_amg_dg_backend.hh:165
Definition: ovlp_amg_dg_backend.hh:304
virtual void pre(V &x, W &b) override
Prepare the preconditioner.
Definition: ovlp_amg_dg_backend.hh:352
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::domain_type::field_type > CGV
Definition: ovlp_amg_dg_backend.hh:322
virtual void apply(V &x, const W &b) override
Apply the precondioner.
Definition: ovlp_amg_dg_backend.hh:366
SolverCategory::Category category() const override
Definition: ovlp_amg_dg_backend.hh:326
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::range_type::field_type > W
Definition: ovlp_amg_dg_backend.hh:319
Backend::Native< V > X
Definition: ovlp_amg_dg_backend.hh:320
Backend::Native< W > Y
Definition: ovlp_amg_dg_backend.hh:321
Dune::PDELab::Backend::Vector< CGGFS, typename CGPrec::range_type::field_type > CGW
Definition: ovlp_amg_dg_backend.hh:323
virtual void post(V &x) override
Clean up.
Definition: ovlp_amg_dg_backend.hh:428
Dune::PDELab::Backend::Vector< DGGFS, typename DGPrec::domain_type::field_type > V
Definition: ovlp_amg_dg_backend.hh:318
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:338
Definition: ovlp_amg_dg_backend.hh:468
Dune::SeqSSOR< CGMatrix, CGVector, CGVector, 1 > Smoother
Definition: ovlp_amg_dg_backend.hh:499
Dune::TransposedMatMultMatResult< P, Matrix >::type PTADG
Definition: ovlp_amg_dg_backend.hh:493
Dune::BlockPreconditioner< CGVector, CGVector, Comm, Smoother > ParSmoother
Definition: ovlp_amg_dg_backend.hh:500
Dune::Amg::AMG< ParCGOperator, CGVector, ParSmoother, Comm > AMG
Definition: ovlp_amg_dg_backend.hh:501
DGGO::Traits::TrialGridFunctionSpace GFS
Definition: ovlp_amg_dg_backend.hh:471
DGGO::Traits::Jacobian M
Definition: ovlp_amg_dg_backend.hh:474
Backend::Native< CGV > CGVector
Definition: ovlp_amg_dg_backend.hh:482
void setDGSmootherRelaxation(double relaxation_)
set number of presmoothing steps on the DG level
Definition: ovlp_amg_dg_backend.hh:673
Dune::PDELab::ISTL::BCRSMatrixBackend MBE
Definition: ovlp_amg_dg_backend.hh:485
Backend::Native< V > Vector
Definition: ovlp_amg_dg_backend.hh:477
Dune::PDELab::Backend::Vector< CGGFS, field_type > CGV
Definition: ovlp_amg_dg_backend.hh:481
Backend::Native< M > Matrix
Definition: ovlp_amg_dg_backend.hh:476
void setNoDGPostSmoothSteps(int n2_)
set number of postsmoothing steps on the DG level
Definition: ovlp_amg_dg_backend.hh:685
bool getReuse() const
Return whether the AMG is reused during call to apply()
Definition: ovlp_amg_dg_backend.hh:580
Vector::field_type field_type
Definition: ovlp_amg_dg_backend.hh:478
Backend::Native< PMatrix > P
Definition: ovlp_amg_dg_backend.hh:490
DGGO::Traits::Domain V
Definition: ovlp_amg_dg_backend.hh:475
Dune::OverlappingSchwarzOperator< CGMatrix, CGVector, CGVector, Comm > ParCGOperator
Definition: ovlp_amg_dg_backend.hh:498
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, unsigned maxiter_=5000, int verbose_=1, bool reuse_=false, bool usesuperlu_=true)
Definition: ovlp_amg_dg_backend.hh:587
Dune::Amg::Parameters Parameters
Definition: ovlp_amg_dg_backend.hh:502
void setNoDGPreSmoothSteps(int n1_)
set number of presmoothing steps on the DG level
Definition: ovlp_amg_dg_backend.hh:679
Dune::PDELab::GridOperator< CGGFS, GFS, CGTODGLOP, MBE, field_type, field_type, field_type, CC, CC > PGO
Definition: ovlp_amg_dg_backend.hh:488
void apply(M &A, V &z, V &r, typename Dune::template FieldTraits< typename V::ElementType >::real_type reduction)
solve the given linear system
Definition: ovlp_amg_dg_backend.hh:698
Dune::MatMultMatResult< PTADG, P >::type CGMatrix
Definition: ovlp_amg_dg_backend.hh:494
PMatrix & prolongation_matrix()
Definition: ovlp_amg_dg_backend.hh:547
void setReuse(bool reuse_)
Set whether the AMG should be reused again during call to apply().
Definition: ovlp_amg_dg_backend.hh:574
TransferLOP CGTODGLOP
Definition: ovlp_amg_dg_backend.hh:487
void setParameters(const Parameters &amg_parameters_)
set AMG parameters
Definition: ovlp_amg_dg_backend.hh:556
Dune::PDELab::ISTL::CommSelector< s, Dune::MPIHelper::isFake >::type Comm
Definition: ovlp_amg_dg_backend.hh:497
Dune::PDELab::EmptyTransformation CC
Definition: ovlp_amg_dg_backend.hh:486
const Parameters & parameters() const
Get the parameters describing the behaviuour of AMG.
Definition: ovlp_amg_dg_backend.hh:568
PGO::Jacobian PMatrix
Definition: ovlp_amg_dg_backend.hh:489
ISTLBackend_OVLP_AMG_4_DG(DGGO &dggo_, const DGCC &dgcc_, CGGFS &cggfs_, const CGCC &cgcc_, const ParameterTree &params)
Definition: ovlp_amg_dg_backend.hh:631
Definition: ovlpistlsolverbackend.hh:43
Definition: ovlpistlsolverbackend.hh:384
const ISTL::ParallelHelper< DGGO::Traits::TrialGridFunctionSpace > & parallelHelper() const
Definition: ovlpistlsolverbackend.hh:414
Definition: ovlpistlsolverbackend.hh:434
Definition: parallelhelper.hh:51
Dune::Amg::SequentialInformation type
Definition: parallelhelper.hh:429
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
Definition: solver.hh:54
Dune::PDELab::LinearSolverResult< double > res
Definition: solver.hh:63
Definition: constraintstransformation.hh:112
Definition: genericdatahandle.hh:667
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
Default flags for all local operators.
Definition: flags.hh:19
Default class for additional methods in instationary local operators.
Definition: idefault.hh:90
Implements linear and nonlinear versions of jacobian_apply_volume() based on alpha_volume()
Definition: numericaljacobianapply.hh:34
sparsity pattern generator
Definition: pattern.hh:14
Definition: recipe-operator-splitting.cc:108