dune-composites (2.5.1)

two_level_schwarz.hh
1#ifndef TWO_LEVEL_SCHWARZ_HH
2#define TWO_LEVEL_SCHWARZ_HH
3
4#if HAVE_ARPACKPP
5
6#include <dune/pdelab/boilerplate/pdelab.hh>
7
8#include <dune/common/timer.hh>
9
10#include "coarsespace.hh"
11#include "coarse_prec.hh"
12#include "../hypre/hypreinterface_sequential.hh"
13
14#include <dune/common/deprecated.hh>
15#include <dune/common/parallel/mpihelper.hh>
16
17#include <dune/istl/owneroverlapcopy.hh>
18#include <dune/istl/solvercategory.hh>
19#include <dune/istl/operators.hh>
20#include <dune/istl/solvers.hh>
21#include <dune/istl/preconditioners.hh>
22#include <dune/istl/scalarproducts.hh>
23#include <dune/istl/paamg/amg.hh>
24#include <dune/istl/paamg/pinfo.hh>
25#include <dune/istl/io.hh>
26#include <dune/istl/superlu.hh>
27#include <dune/istl/umfpack.hh>
28
29#include <dune/pdelab/constraints/common/constraints.hh>
30#include <dune/pdelab/gridfunctionspace/genericdatahandle.hh>
31#include <dune/pdelab/backend/solver.hh>
32#include <dune/pdelab/backend/istl/vector.hh>
33#include <dune/pdelab/backend/istl/bcrsmatrix.hh>
34
35// myModel.template solve<GO,V,GFS,C,Constraints,MBE,LOP>(go,u,gfs,cg,constraints,mbe,lop);
36
37namespace Dune{
38 namespace Geneo{
39
42 template<class GFS, class M, class X, class Y>
43 class TwoLevelOverlappingAdditiveSchwarzDC
44 : public Dune::Preconditioner<X,Y>
45 {
46 public:
47 typedef Dune::PDELab::Backend::Native<M> ISTLM;
48
49 typedef Dune::BlockVector<Dune::FieldVector<double,1> > COARSE_V;
50 typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > COARSE_M;
51
52 // define the category
53 virtual Dune::SolverCategory::Category category() const
54 {
55 return Dune::SolverCategory::overlapping;
56 }
57
65 TwoLevelOverlappingAdditiveSchwarzDC (const GFS& gfs_, const M& AF, std::shared_ptr<CoarseSpace<M,X> > coarse_space, bool coarseSpaceActive_ = true)
66 : gfs(gfs_),
67 solverf(Dune::PDELab::Backend::native(AF),false),
68 my_rank(gfs.gridView().comm().rank()),
69 coarse_space_(coarse_space) ,
70 coarse_solver_(*coarse_space_->get_coarse_system()),
71 coarseSpaceActive(coarseSpaceActive_)
72 {}
73
79 virtual void pre (X& x, Y& b)
80 { }
81
87 double coarse_time = 0.0;
88 int apply_calls = 0;
89
90 virtual void apply (X& v, const Y& d)
91 {
92 // first the subdomain solves
93 Y b(d); // need copy, since solver overwrites right hand side
94 Dune::InverseOperatorResult result;
95 solverf.apply(v,b,result);
96 if (!coarseSpaceActive) {
97
98 Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,v);
99 // Just add local results and return in 1-level Schwarz case
100 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
101
102 } else {
103
104 MPI_Barrier(gfs.gridView().comm());
105 Dune::Timer timer_setup;
106
107 Dune::InverseOperatorResult result;
108
109 auto coarse_defect = coarse_space_->restrict_defect (d);
110
111 // Solve coarse system
112 COARSE_V v0(coarse_space_->basis_size(),coarse_space_->basis_size());
113
114 //coarse_solver_2->apply(v0, *coarse_defect, result);
115 coarse_solver_.apply(v0,*coarse_defect,result);
116
117 // Prolongate coarse solution on local domain
118 auto coarse_correction = coarse_space_->prolongate_defect (v0);
119 v += *coarse_correction;
120
121 coarse_time += timer_setup.elapsed();
122 apply_calls++;
123
124 Dune::PDELab::AddDataHandle<GFS,X> result_addh(gfs,v);
125 gfs.gridView().communicate(result_addh,Dune::All_All_Interface,Dune::ForwardCommunication);
126 }
127 }
128
134 virtual void post (X& x) {
135 if (my_rank == 0) std::cout << "Coarse time CT=" << coarse_time << std::endl;
136 if (my_rank == 0) std::cout << "Coarse time per apply CTA=" << coarse_time / (double)apply_calls << std::endl;
137 }
138
139 private:
140
141 const GFS& gfs;
142 Dune::UMFPack<ISTLM> solverf;
143 int my_rank;
144 std::shared_ptr<CoarseSpace<M,X> > coarse_space_;
145 Dune::UMFPack<COARSE_M> coarse_solver_;
146 bool coarseSpaceActive;
147 };
148
149
152 template<class GFS, class M, class X, class Y>
153 class TwoLevelOverlappingAdditiveSchwarzAMG
154 : public Dune::Preconditioner<X,Y>
155 {
156 public:
157 typedef Dune::PDELab::Backend::Native<M> ISTLM;
158
159 typedef Dune::BlockVector<Dune::FieldVector<double,1> > COARSE_V;
160 typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > COARSE_M;
161
162 // define the category
163 virtual Dune::SolverCategory::Category category() const
164 {
165 return Dune::SolverCategory::overlapping;
166 }
167
175 TwoLevelOverlappingAdditiveSchwarzAMG (const GFS& gfs_, const M& AF, std::shared_ptr<CoarseSpace<M,X> > coarse_space, bool coarseSpaceActive_ = true)
176 : gfs(gfs_),
177 solverf(Dune::PDELab::Backend::native(AF),false),
178 my_rank(gfs.gridView().comm().rank()),
179 coarse_space_(coarse_space) ,
180 coarseSpaceActive(coarseSpaceActive_)
181 {}
182
188 virtual void pre (X& x, Y& b)
189 { }
190
196 double coarse_time = 0.0;
197 int apply_calls = 0;
198
199 virtual void apply (X& v, const Y& d)
200 {
201 // first the subdomain solves
202 Y b(d); // need copy, since solver overwrites right hand side
203 Dune::InverseOperatorResult result;
204 solverf.apply(v,b,result);
205 if (!coarseSpaceActive) {
206
207 Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,v);
208 // Just add local results and return in 1-level Schwarz case
209 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
210
211 } else {
212
213 MPI_Barrier(gfs.gridView().comm());
214 Dune::Timer timer_setup;
215
216 Dune::InverseOperatorResult result;
217
218 auto coarse_defect = coarse_space_->restrict_defect (d);
219
220 // Solve coarse system
221 COARSE_V v0(coarse_space_->basis_size(),coarse_space_->basis_size());
222
223 COARSE_M a = *coarse_space_->get_coarse_system();
224 typedef typename GFS::Traits::GridView GV;
225 GV gv = gfs.gridView();
226#if HAVE_HYPRE
227 SequentialHypreSolver<GV,COARSE_M, COARSE_V> solver(gv, v0,a,*coarse_defect, 1e-3);
228 solver.solve(v0);
229#endif
230 // Prolongate coarse solution on local domain
231 auto coarse_correction = coarse_space_->prolongate_defect (v0);
232 v += *coarse_correction;
233
234 coarse_time += timer_setup.elapsed();
235 apply_calls++;
236
237 Dune::PDELab::AddDataHandle<GFS,X> result_addh(gfs,v);
238 gfs.gridView().communicate(result_addh,Dune::All_All_Interface,Dune::ForwardCommunication);
239 }
240 }
241
247 virtual void post (X& x) {
248 if (my_rank == 0) std::cout << "Coarse time CT=" << coarse_time << std::endl;
249 if (my_rank == 0) std::cout << "Coarse time per apply CTA=" << coarse_time / (double)apply_calls << std::endl;
250 }
251
252 private:
253
254 const GFS& gfs;
255 Dune::UMFPack<ISTLM> solverf;
256 int my_rank;
257 std::shared_ptr<CoarseSpace<M,X> > coarse_space_;
258 bool coarseSpaceActive;
259 };
260
261
264 template<class GFS, class M, class X, class Y>
265 class TwoLevelOverlappingAdditiveSchwarzPCGC
266 : public Dune::Preconditioner<X,Y>
267 {
268 public:
269 typedef Dune::PDELab::Backend::Native<M> ISTLM;
270
271 typedef Dune::BlockVector<Dune::FieldVector<double,1> > COARSE_V;
272 typedef Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> > COARSE_M;
273
274
275 //Set up coarse solver
276 typedef typename Dune::CGSolver<COARSE_V> CSOLVER;
277 typedef CoarsePrec<GFS,COARSE_M,COARSE_V, COARSE_V,M,X,Y> PREC;
278
279 // define the category
280 virtual Dune::SolverCategory::Category category() const
281 {
282 return Dune::SolverCategory::overlapping;
283 }
284
292 TwoLevelOverlappingAdditiveSchwarzPCGC (const GFS& gfs_, const M& AF, std::shared_ptr<CoarseSpace<M,X> > coarse_space, std::shared_ptr<CoarseSpace<M,X>> coarse_space_small, bool coarseSpaceActive_ = true)
293 : gfs(gfs_),
294 solverf(Dune::PDELab::Backend::native(AF),false),
295 my_rank(gfs.gridView().comm().rank()),
296 coarse_space_(coarse_space),
297 coarse_space_small_(coarse_space_small),
298 prec(gfs,*coarse_space_small_->get_coarse_system(), coarse_space_,coarse_space_small_),
299 opa(*coarse_space_->get_coarse_system()),
300 coarseSpaceActive(coarseSpaceActive_)
301 {
302 using Dune::PDELab::Backend::native;
303 int verb = 0;
304 if(my_rank == 0) verb = 1;
305 coarse_solver_2 = std::make_shared<CSOLVER>(opa,prec,1e-3,5000,verb);
306 }
307
313 virtual void pre (X& x, Y& b)
314 { }
315
321 double coarse_time = 0.0;
322 int apply_calls = 0;
323
324 virtual void apply (X& v, const Y& d)
325 {
326 // first the subdomain solves
327 Y b(d); // need copy, since solver overwrites right hand side
328 Dune::InverseOperatorResult result;
329 solverf.apply(v,b,result);
330
331 if (!coarseSpaceActive) {
332
333 Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,v);
334 // Just add local results and return in 1-level Schwarz case
335 gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
336
337 } else {
338
339 MPI_Barrier(gfs.gridView().comm());
340 Dune::Timer timer_setup;
341
342 Dune::InverseOperatorResult result;
343
344 auto coarse_defect = coarse_space_->restrict_defect (d);
345
346 // Solve coarse system
347 COARSE_V v0(coarse_space_->basis_size(),coarse_space_->basis_size());
348
349 coarse_solver_2->apply(v0, *coarse_defect, result);
350
351 // Prolongate coarse solution on local domain
352 auto coarse_correction = coarse_space_->prolongate_defect (v0);
353 v += *coarse_correction;
354
355 coarse_time += timer_setup.elapsed();
356 apply_calls++;
357
358 Dune::PDELab::AddDataHandle<GFS,X> result_addh(gfs,v);
359 gfs.gridView().communicate(result_addh,Dune::All_All_Interface,Dune::ForwardCommunication);
360 }
361 }
362
368 virtual void post (X& x) {
369 if (my_rank == 0) std::cout << "Coarse time CT=" << coarse_time << std::endl;
370 if (my_rank == 0) std::cout << "Coarse time per apply CTA=" << coarse_time / (double)apply_calls << std::endl;
371 }
372
373 private:
374
375 const GFS& gfs;
376 Dune::UMFPack<ISTLM> solverf;
377 int my_rank;
378 std::shared_ptr<CoarseSpace<M,X> > coarse_space_;
379 std::shared_ptr<CoarseSpace<M,X> > coarse_space_small_;
380 Dune::MatrixAdapter<COARSE_M,COARSE_V,COARSE_V> opa;
381 PREC prec;
382 std::shared_ptr<CSOLVER> coarse_solver_2;
383 bool coarseSpaceActive;
384 };
385 }
386}
387
388
389#endif // have arpack
390
391#endif
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 5, 22:35, 2025)