OpenFPM_pdata  4.1.0
Project that contain the implementation of distributed structures
 
Loading...
Searching...
No Matches
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:81
solError::err_inf
PetscReal err_inf
infinity norm of the error
Definition petsc_solver.hpp:83
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