OpenFPM_pdata  1.1.0
Project that contain the implementation of distributed structures
 All Data Structures Namespaces Functions Variables Typedefs Enumerations Friends Pages
petsc_solver_ASM_AMG_to_do.hpp
1 /*
2  * petsc_solver_ASM_AMG_to_do.hpp
3  *
4  * Created on: Jun 3, 2016
5  * Author: i-bird
6  */
7 
8 #ifndef OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_ASM_AMG_TO_DO_HPP_
9 #define OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_ASM_AMG_TO_DO_HPP_
10 
11 
17  void try_solve_complex_bj(Mat & A_, const Vec & b_, Vec & x_)
18  {
19  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,5));
20  PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
21 
22  solve_complex(A_,b_,x_);
23  }
24 
25 
31  void try_solve_ASM(Mat & A_, const Vec & b_, Vec & x_)
32  {
33  PETSC_SAFE_CALL(KSPSetTolerances(ksp,rtol,abstol,dtol,5));
34  for (size_t i = 0 ; i < solvs.size() ; i++)
35  {
36  setSolver(solvs.get(i).c_str());
37 
38  // Disable convergence check
39  PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
40 
41 // solve_simple(A_,b_,x_);
42  solve_ASM(A_,b_,x_);
43  }
44  }
45 
46 
53  void solve_complex(Mat & A_, const Vec & b_, Vec & x_)
54  {
55  // We get the size of the Matrix A
56  PetscInt row;
57  PetscInt col;
58  PetscInt row_loc;
59  PetscInt col_loc;
60  PetscInt * blks;
61  PetscInt nlocal;
62  PetscInt first;
63  KSP *subksp;
64  PC subpc;
65 
66  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
67  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
68 
69  // Set the preconditioner to Block-jacobi
70  PC pc;
71  KSPGetPC(ksp,&pc);
72  PCSetType(pc,PCBJACOBI);
73  PETSC_SAFE_CALL(KSPSetType(ksp,KSPGMRES));
74 
75  PetscInt m = 1;
76 
77  PetscMalloc1(m,&blks);
78  for (size_t i = 0 ; i < m ; i++) blks[i] = row_loc;
79 
80  PCBJacobiSetLocalBlocks(pc,m,blks);
81  PetscFree(blks);
82 
83  KSPSetUp(ksp);
84  PCSetUp(pc);
85 
86  PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);
87 
88  for (size_t i = 0; i < nlocal; i++)
89  {
90  KSPGetPC(subksp[i],&subpc);
91 
92  PCSetType(subpc,PCLU);
93 
94 // PCSetType(subpc,PCJACOBI);
95  KSPSetType(subksp[i],KSPPREONLY);
96 // PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedDefault,NULL,NULL));
97 // KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
98 
99 /* if (!rank)
100  {
101  if (i%2)
102  {
103  PCSetType(subpc,PCILU);
104  }
105  else
106  {
107  PCSetType(subpc,PCNONE);
108  KSPSetType(subksp[i],KSPBCGS);
109  KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
110  }
111  }
112  else
113  {
114  PCSetType(subpc,PCJACOBI);
115  KSPSetType(subksp[i],KSPGMRES);
116  KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
117  }*/
118  }
119 
120 
121  KSPSolve(ksp,b_,x_);
122 
123  auto & v_cl = create_vcluster();
124  if (try_solve == true)
125  {
126  // calculate error statistic about the solution
127  solError err = statSolutionError(A_,b_,x_);
128 
129  if (v_cl.getProcessUnitID() == 0)
130  {
131  std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
132  std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
133  }
134  }
135  }
136 
137 
138  void solve_ASM(Mat & A_, const Vec & b_, Vec & x_)
139  {
140  PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
141 
142  // We set the size of x according to the Matrix A
143  PetscInt row;
144  PetscInt col;
145  PetscInt row_loc;
146  PetscInt col_loc;
147 
148  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
149  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
150 
151  PC pc;
152 
153  // We set the Matrix operators
154  PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
155 
156  // We set the pre-conditioner
157  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
158  PETSC_SAFE_CALL(PCSetType(pc,PCASM));
159 
161 
163 
164 
165 // PCASMSetLocalSubdomains(pc);
166 
167  PCASMSetOverlap(pc,5);
168 
169  KSP *subksp; /* array of KSP contexts for local subblocks */
170  PetscInt nlocal,first; /* number of local subblocks, first local subblock */
171  PC subpc; /* PC context for subblock */
172  PetscBool isasm;
173 
174  /*
175  Set runtime options
176  */
177 // KSPSetFromOptions(ksp);
178 
179  /*
180  Flag an error if PCTYPE is changed from the runtime options
181  */
182 // PetscObjectTypeCompare((PetscObject)pc,PCASM,&isasm);
183 // if (!isasm) SETERRQ(PETSC_COMM_WORLD,1,"Cannot Change the PCTYPE when manually changing the subdomain solver settings");
184 
185  /*
186  Call KSPSetUp() to set the block Jacobi data structures (including
187  creation of an internal KSP context for each block).
188 
189  Note: KSPSetUp() MUST be called before PCASMGetSubKSP().
190  */
191  KSPSetUp(ksp);
192 
193  /*
194  Extract the array of KSP contexts for the local blocks
195  */
196  PCASMGetSubKSP(pc,&nlocal,&first,&subksp);
197 
198  /*
199  Loop over the local blocks, setting various KSP options
200  for each block.
201  */
202  for (size_t i = 0 ; i < nlocal ; i++)
203  {
204  KSPGetPC(subksp[i],&subpc);
205 // PCFactorSetFill(subpc,30);
206  PCFactorSetLevels(subpc,5);
207  PCSetType(subpc,PCILU);
208  KSPSetType(subksp[i],KSPRICHARDSON);
209  KSPSetTolerances(subksp[i],1.e-3,0.1,PETSC_DEFAULT,PETSC_DEFAULT);
210  }
211 
212  // if we are on on best solve set-up a monitor function
213 
214  if (try_solve == true)
215  {
216  // for bench-mark we are interested in non-preconditioned norm
217  PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
218 
219  // Disable convergence check
220  PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
221  }
222 
223 // KSPGMRESSetRestart(ksp,100);
224 
225  // Solve the system
226  PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
227 
228  KSPConvergedReason reason;
229  KSPGetConvergedReason(ksp,&reason);
230 
231  std::cout << "Reason: " << reason << std::endl;
232 
233  auto & v_cl = create_vcluster();
234 // if (try_solve == true)
235 // {
236  // calculate error statistic about the solution
237  solError err = statSolutionError(A_,b_,x_);
238 
239  if (v_cl.getProcessUnitID() == 0)
240  {
241  std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
242  std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
243  }
244 // }
245 
246  }
247 
248  void solve_AMG(Mat & A_, const Vec & b_, Vec & x_)
249  {
250  // PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
251  PETSC_SAFE_CALL(KSPSetType(ksp,KSPRICHARDSON));
252 
253  // -pc_hypre_boomeramg_tol <tol>
254 
255  // We set the size of x according to the Matrix A
256  PetscInt row;
257  PetscInt col;
258  PetscInt row_loc;
259  PetscInt col_loc;
260 
261  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
262  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
263 
264  PC pc;
265 
266  // We set the Matrix operators
267  PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
268 
269  // We set the pre-conditioner
270  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
271 
272  // PETSC_SAFE_CALL(PCSetType(pc,PCJACOBI));
273 
274  PETSC_SAFE_CALL(PCSetType(pc,PCHYPRE));
275  // PCGAMGSetNSmooths(pc,0);
276  PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
277  PCFactorSetShiftAmount(pc, PETSC_DECIDE);
278  PCHYPRESetType(pc, "boomeramg");
279  MatSetBlockSize(A_,4);
280  PetscOptionsSetValue("-pc_hypre_boomeramg_print_statistics","2");
281  PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter","1000");
282  PetscOptionsSetValue("-pc_hypre_boomeramg_nodal_coarsen","true");
283  PetscOptionsSetValue("-pc_hypre_boomeramg_relax_type_all","SOR/Jacobi");
284  PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type","Falgout");
285  PetscOptionsSetValue("-pc_hypre_boomeramg_cycle_type","W");
286  PetscOptionsSetValue("-pc_hypre_boomeramg_max_levels","20");
287  PetscOptionsSetValue("-pc_hypre_boomeramg_vec_interp_variant","0");
288  KSPSetFromOptions(ksp);
289 
290  // if we are on on best solve set-up a monitor function
291 
292  if (try_solve == true)
293  {
294  // for bench-mark we are interested in non-preconditioned norm
295  PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
296 
297  // Disable convergence check
298  PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
299  }
300 
301  // Solve the system
302  PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
303 
304  auto & v_cl = create_vcluster();
305  // if (try_solve == true)
306  // {
307  // calculate error statistic about the solution
308  solError err = statSolutionError(A_,b_,x_);
309 
310  if (v_cl.getProcessUnitID() == 0)
311  {
312  std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
313  std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
314  }
315  // }
316  }
317 
318 
322  void solve_Krylov_simple(Mat & A_, const Vec & b_, Vec & x_)
323  {
324  PETSC_SAFE_CALL(KSPSetType(ksp,s_type));
325 
326  // We set the size of x according to the Matrix A
327  PetscInt row;
328  PetscInt col;
329  PetscInt row_loc;
330  PetscInt col_loc;
331 
332  PETSC_SAFE_CALL(MatGetSize(A_,&row,&col));
333  PETSC_SAFE_CALL(MatGetLocalSize(A_,&row_loc,&col_loc));
334 
335  PC pc;
336 
337  // We set the Matrix operators
338  PETSC_SAFE_CALL(KSPSetOperators(ksp,A_,A_));
339 
340  // We set the pre-conditioner
341  PETSC_SAFE_CALL(KSPGetPC(ksp,&pc));
342 
343  // PETSC_SAFE_CALL(PCSetType(pc,PCJACOBI));
344 
345  PETSC_SAFE_CALL(PCSetType(pc,PCHYPRE));
346 // PCGAMGSetNSmooths(pc,0);
347  PCFactorSetShiftType(pc, MAT_SHIFT_NONZERO);
348  PCFactorSetShiftAmount(pc, PETSC_DECIDE);
349  PCHYPRESetType(pc, "boomeramg");
350  MatSetBlockSize(A_,4);
351  PetscOptionsSetValue("-pc_hypre_boomeramg_print_statistics","2");
352  PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter","1000");
353  PetscOptionsSetValue("-pc_hypre_boomeramg_nodal_coarsen","true");
354  PetscOptionsSetValue("-pc_hypre_boomeramg_relax_type_all","SOR/Jacobi");
355  PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type","Falgout");
356  PetscOptionsSetValue("-pc_hypre_boomeramg_cycle_type","W");
357  PetscOptionsSetValue("-pc_hypre_boomeramg_max_levels","10");
358  KSPSetFromOptions(ksp);
359  // if we are on on best solve set-up a monitor function
360 
361  if (try_solve == true)
362  {
363  // for bench-mark we are interested in non-preconditioned norm
364  PETSC_SAFE_CALL(KSPMonitorSet(ksp,monitor,&vres,NULL));
365 
366  // Disable convergence check
367  PETSC_SAFE_CALL(KSPSetConvergenceTest(ksp,KSPConvergedSkip,NULL,NULL));
368  }
369 
370  // Solve the system
371  PETSC_SAFE_CALL(KSPSolve(ksp,b_,x_));
372 
373  auto & v_cl = create_vcluster();
374 // if (try_solve == true)
375 // {
376  // calculate error statistic about the solution
377  solError err = statSolutionError(A_,b_,x_);
378 
379  if (v_cl.getProcessUnitID() == 0)
380  {
381  std::cout << "Method: " << s_type << " " << " pre-conditoner: " << PCJACOBI << " iterations: " << err.its << std::endl;
382  std::cout << "Norm of error: " << err.err_norm << " Norm infinity: " << err.err_inf << std::endl;
383  }
384 // }
385  }
386 
388 
389  IS *is,*is_local;
390  PetscInt Nsub;
391 
392  PCASMCreateSubdomains2D(128,128,4,4,1,1,&Nsub,&is,&is_local);
393 
394 
395  if (create_vcluster().getProcessUnitID() == 1)
396  ISView(is_local[1],PETSC_VIEWER_STDOUT_SELF);
397 
399 
400 #endif /* OPENFPM_NUMERICS_SRC_SOLVERS_PETSC_SOLVER_ASM_AMG_TO_DO_HPP_ */
solError
It contain statistic of the error of the calculated solution.
Definition: petsc_solver.hpp:80
solError::its
PetscInt its
Number of iterations.
Definition: petsc_solver.hpp:89
solError::err_norm
PetscReal err_norm
L1 norm of the error.
Definition: petsc_solver.hpp:86
solError::err_inf
PetscReal err_inf
infinity norm of the error
Definition: petsc_solver.hpp:83