Actual source code: hypre.c
1: /*
2: Provides an interface to the LLNL package hypre
3: */
5: #include <petscpkg_version.h>
6: #include <petsc/private/pcimpl.h>
7: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
8: #include <petsc/private/matimpl.h>
9: #include <petsc/private/vecimpl.h>
10: #include <../src/vec/vec/impls/hypre/vhyp.h>
11: #include <../src/mat/impls/hypre/mhypre.h>
12: #include <../src/dm/impls/da/hypre/mhyp.h>
13: #include <_hypre_parcsr_ls.h>
14: #include <petscmathypre.h>
16: #if defined(PETSC_HAVE_HYPRE_DEVICE)
17: #include <petsc/private/deviceimpl.h>
18: #endif
20: static PetscBool cite = PETSC_FALSE;
21: static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = {\\url{https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods}}\n}\n";
23: /*
24: Private context (data structure) for the preconditioner.
25: */
26: typedef struct {
27: HYPRE_Solver hsolver;
28: Mat hpmat; /* MatHYPRE */
30: HYPRE_Int (*destroy)(HYPRE_Solver);
31: HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
32: HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
34: MPI_Comm comm_hypre;
35: char *hypre_type;
37: /* options for Pilut and BoomerAMG*/
38: PetscInt maxiter;
39: PetscReal tol;
41: /* options for Pilut */
42: PetscInt factorrowsize;
44: /* options for ParaSails */
45: PetscInt nlevels;
46: PetscReal threshold;
47: PetscReal filter;
48: PetscReal loadbal;
49: PetscInt logging;
50: PetscInt ruse;
51: PetscInt symt;
53: /* options for BoomerAMG */
54: PetscBool printstatistics;
56: /* options for BoomerAMG */
57: PetscInt cycletype;
58: PetscInt maxlevels;
59: PetscReal strongthreshold;
60: PetscReal maxrowsum;
61: PetscInt gridsweeps[3];
62: PetscInt coarsentype;
63: PetscInt measuretype;
64: PetscInt smoothtype;
65: PetscInt smoothnumlevels;
66: PetscInt eu_level; /* Number of levels for ILU(k) in Euclid */
67: PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
68: PetscInt eu_bj; /* Defines use of Block Jacobi ILU in Euclid */
69: PetscInt relaxtype[3];
70: PetscReal relaxweight;
71: PetscReal outerrelaxweight;
72: PetscInt relaxorder;
73: PetscReal truncfactor;
74: PetscBool applyrichardson;
75: PetscInt pmax;
76: PetscInt interptype;
77: PetscInt maxc;
78: PetscInt minc;
79: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
80: char *spgemm_type; // this is a global hypre parameter but is closely associated with BoomerAMG
81: #endif
82: /* GPU */
83: PetscBool keeptranspose;
84: PetscInt rap2;
85: PetscInt mod_rap2;
87: /* AIR */
88: PetscInt Rtype;
89: PetscReal Rstrongthreshold;
90: PetscReal Rfilterthreshold;
91: PetscInt Adroptype;
92: PetscReal Adroptol;
94: PetscInt agg_nl;
95: PetscInt agg_interptype;
96: PetscInt agg_num_paths;
97: PetscBool nodal_relax;
98: PetscInt nodal_relax_levels;
100: PetscInt nodal_coarsening;
101: PetscInt nodal_coarsening_diag;
102: PetscInt vec_interp_variant;
103: PetscInt vec_interp_qmax;
104: PetscBool vec_interp_smooth;
105: PetscInt interp_refine;
107: /* NearNullSpace support */
108: VecHYPRE_IJVector *hmnull;
109: HYPRE_ParVector *phmnull;
110: PetscInt n_hmnull;
111: Vec hmnull_constant;
113: /* options for AS (Auxiliary Space preconditioners) */
114: PetscInt as_print;
115: PetscInt as_max_iter;
116: PetscReal as_tol;
117: PetscInt as_relax_type;
118: PetscInt as_relax_times;
119: PetscReal as_relax_weight;
120: PetscReal as_omega;
121: PetscInt as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
122: PetscReal as_amg_alpha_theta; /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
123: PetscInt as_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
124: PetscReal as_amg_beta_theta; /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS) */
125: PetscInt ams_cycle_type;
126: PetscInt ads_cycle_type;
128: /* additional data */
129: Mat G; /* MatHYPRE */
130: Mat C; /* MatHYPRE */
131: Mat alpha_Poisson; /* MatHYPRE */
132: Mat beta_Poisson; /* MatHYPRE */
134: /* extra information for AMS */
135: PetscInt dim; /* geometrical dimension */
136: VecHYPRE_IJVector coords[3];
137: VecHYPRE_IJVector constants[3];
138: Mat RT_PiFull, RT_Pi[3];
139: Mat ND_PiFull, ND_Pi[3];
140: PetscBool ams_beta_is_zero;
141: PetscBool ams_beta_is_zero_part;
142: PetscInt ams_proj_freq;
143: } PC_HYPRE;
145: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
146: {
147: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
149: *hsolver = jac->hsolver;
150: return 0;
151: }
153: /*
154: Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
155: is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
156: It is used in PCHMG. Other users should avoid using this function.
157: */
158: static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc,PetscInt *nlevels,Mat *operators[])
159: {
160: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
161: PetscBool same = PETSC_FALSE;
162: PetscInt num_levels,l;
163: Mat *mattmp;
164: hypre_ParCSRMatrix **A_array;
166: PetscStrcmp(jac->hypre_type,"boomeramg",&same);
168: num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
169: PetscMalloc1(num_levels,&mattmp);
170: A_array = hypre_ParAMGDataAArray((hypre_ParAMGData*) (jac->hsolver));
171: for (l=1; l<num_levels; l++) {
172: MatCreateFromParCSR(A_array[l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[num_levels-1-l]));
173: /* We want to own the data, and HYPRE can not touch this matrix any more */
174: A_array[l] = NULL;
175: }
176: *nlevels = num_levels;
177: *operators = mattmp;
178: return 0;
179: }
181: /*
182: Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
183: is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
184: It is used in PCHMG. Other users should avoid using this function.
185: */
186: static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc,PetscInt *nlevels,Mat *interpolations[])
187: {
188: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
189: PetscBool same = PETSC_FALSE;
190: PetscInt num_levels,l;
191: Mat *mattmp;
192: hypre_ParCSRMatrix **P_array;
194: PetscStrcmp(jac->hypre_type,"boomeramg",&same);
196: num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
197: PetscMalloc1(num_levels,&mattmp);
198: P_array = hypre_ParAMGDataPArray((hypre_ParAMGData*) (jac->hsolver));
199: for (l=1; l<num_levels; l++) {
200: MatCreateFromParCSR(P_array[num_levels-1-l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[l-1]));
201: /* We want to own the data, and HYPRE can not touch this matrix any more */
202: P_array[num_levels-1-l] = NULL;
203: }
204: *nlevels = num_levels;
205: *interpolations = mattmp;
206: return 0;
207: }
209: /* Resets (frees) Hypre's representation of the near null space */
210: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
211: {
212: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
213: PetscInt i;
215: for (i=0; i<jac->n_hmnull; i++) {
216: VecHYPRE_IJVectorDestroy(&jac->hmnull[i]);
217: }
218: PetscFree(jac->hmnull);
219: PetscFree(jac->phmnull);
220: VecDestroy(&jac->hmnull_constant);
221: jac->n_hmnull = 0;
222: return 0;
223: }
225: static PetscErrorCode PCSetUp_HYPRE(PC pc)
226: {
227: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
228: Mat_HYPRE *hjac;
229: HYPRE_ParCSRMatrix hmat;
230: HYPRE_ParVector bv,xv;
231: PetscBool ishypre;
233: if (!jac->hypre_type) {
234: PCHYPRESetType(pc,"boomeramg");
235: }
237: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
238: if (!ishypre) {
239: MatDestroy(&jac->hpmat);
240: MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);
241: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hpmat);
242: } else {
243: PetscObjectReference((PetscObject)pc->pmat);
244: MatDestroy(&jac->hpmat);
245: jac->hpmat = pc->pmat;
246: }
247: /* allow debug */
248: MatViewFromOptions(jac->hpmat,NULL,"-pc_hypre_mat_view");
249: hjac = (Mat_HYPRE*)(jac->hpmat->data);
251: /* special case for BoomerAMG */
252: if (jac->setup == HYPRE_BoomerAMGSetup) {
253: MatNullSpace mnull;
254: PetscBool has_const;
255: PetscInt bs,nvec,i;
256: const Vec *vecs;
258: MatGetBlockSize(pc->pmat,&bs);
259: if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,jac->hsolver,bs);
260: MatGetNearNullSpace(pc->mat, &mnull);
261: if (mnull) {
262: PCHYPREResetNearNullSpace_Private(pc);
263: MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
264: PetscMalloc1(nvec+1,&jac->hmnull);
265: PetscMalloc1(nvec+1,&jac->phmnull);
266: for (i=0; i<nvec; i++) {
267: VecHYPRE_IJVectorCreate(vecs[i]->map,&jac->hmnull[i]);
268: VecHYPRE_IJVectorCopy(vecs[i],jac->hmnull[i]);
269: PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->hmnull[i]->ij,(void**)&jac->phmnull[i]);
270: }
271: if (has_const) {
272: MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
273: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hmnull_constant);
274: VecSet(jac->hmnull_constant,1);
275: VecNormalize(jac->hmnull_constant,NULL);
276: VecHYPRE_IJVectorCreate(jac->hmnull_constant->map,&jac->hmnull[nvec]);
277: VecHYPRE_IJVectorCopy(jac->hmnull_constant,jac->hmnull[nvec]);
278: PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->hmnull[nvec]->ij,(void**)&jac->phmnull[nvec]);
279: nvec++;
280: }
281: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,jac->hsolver,nvec,jac->phmnull);
282: jac->n_hmnull = nvec;
283: }
284: }
286: /* special case for AMS */
287: if (jac->setup == HYPRE_AMSSetup) {
288: Mat_HYPRE *hm;
289: HYPRE_ParCSRMatrix parcsr;
290: if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
291: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the edge constant vectors via PCHYPRESetEdgeConstantVectors() or the interpolation matrix via PCHYPRESetInterpolations");
292: }
293: if (jac->dim) {
294: PetscStackCallStandard(HYPRE_AMSSetDimension,jac->hsolver,jac->dim);
295: }
296: if (jac->constants[0]) {
297: HYPRE_ParVector ozz,zoz,zzo = NULL;
298: PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->constants[0]->ij,(void**)(&ozz));
299: PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->constants[1]->ij,(void**)(&zoz));
300: if (jac->constants[2]) {
301: PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->constants[2]->ij,(void**)(&zzo));
302: }
303: PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,jac->hsolver,ozz,zoz,zzo);
304: }
305: if (jac->coords[0]) {
306: HYPRE_ParVector coords[3];
307: coords[0] = NULL;
308: coords[1] = NULL;
309: coords[2] = NULL;
310: if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[0]->ij,(void**)(&coords[0]));
311: if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[1]->ij,(void**)(&coords[1]));
312: if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[2]->ij,(void**)(&coords[2]));
313: PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,jac->hsolver,coords[0],coords[1],coords[2]);
314: }
316: hm = (Mat_HYPRE*)(jac->G->data);
317: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr));
318: PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,jac->hsolver,parcsr);
319: if (jac->alpha_Poisson) {
320: hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
321: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr));
322: PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,jac->hsolver,parcsr);
323: }
324: if (jac->ams_beta_is_zero) {
325: PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,jac->hsolver,NULL);
326: } else if (jac->beta_Poisson) {
327: hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
328: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr));
329: PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,jac->hsolver,parcsr);
330: }
331: if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
332: PetscInt i;
333: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
334: if (jac->ND_PiFull) {
335: hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
336: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsrfull));
337: } else {
338: nd_parcsrfull = NULL;
339: }
340: for (i=0;i<3;++i) {
341: if (jac->ND_Pi[i]) {
342: hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
343: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsr[i]));
344: } else {
345: nd_parcsr[i] = NULL;
346: }
347: }
348: PetscStackCallStandard(HYPRE_AMSSetInterpolations,jac->hsolver,nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]);
349: }
350: }
351: /* special case for ADS */
352: if (jac->setup == HYPRE_ADSSetup) {
353: Mat_HYPRE *hm;
354: HYPRE_ParCSRMatrix parcsr;
355: if (!jac->coords[0] && !((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])))) {
356: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
357: }
361: if (jac->coords[0]) {
362: HYPRE_ParVector coords[3];
363: coords[0] = NULL;
364: coords[1] = NULL;
365: coords[2] = NULL;
366: if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[0]->ij,(void**)(&coords[0]));
367: if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[1]->ij,(void**)(&coords[1]));
368: if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,jac->coords[2]->ij,(void**)(&coords[2]));
369: PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,jac->hsolver,coords[0],coords[1],coords[2]);
370: }
371: hm = (Mat_HYPRE*)(jac->G->data);
372: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr));
373: PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,jac->hsolver,parcsr);
374: hm = (Mat_HYPRE*)(jac->C->data);
375: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr));
376: PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,jac->hsolver,parcsr);
377: if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
378: PetscInt i;
379: HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
380: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
381: if (jac->RT_PiFull) {
382: hm = (Mat_HYPRE*)(jac->RT_PiFull->data);
383: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&rt_parcsrfull));
384: } else {
385: rt_parcsrfull = NULL;
386: }
387: for (i=0;i<3;++i) {
388: if (jac->RT_Pi[i]) {
389: hm = (Mat_HYPRE*)(jac->RT_Pi[i]->data);
390: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&rt_parcsr[i]));
391: } else {
392: rt_parcsr[i] = NULL;
393: }
394: }
395: if (jac->ND_PiFull) {
396: hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
397: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsrfull));
398: } else {
399: nd_parcsrfull = NULL;
400: }
401: for (i=0;i<3;++i) {
402: if (jac->ND_Pi[i]) {
403: hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
404: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsr[i]));
405: } else {
406: nd_parcsr[i] = NULL;
407: }
408: }
409: PetscStackCallStandard(HYPRE_ADSSetInterpolations,jac->hsolver,rt_parcsrfull,rt_parcsr[0],rt_parcsr[1],rt_parcsr[2],nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]);
410: }
411: }
412: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat);
413: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&bv);
414: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&xv);
415: PetscStackCallStandard(jac->setup,jac->hsolver,hmat,bv,xv);
416: return 0;
417: }
419: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
420: {
421: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
422: Mat_HYPRE *hjac = (Mat_HYPRE*)(jac->hpmat->data);
423: HYPRE_ParCSRMatrix hmat;
424: HYPRE_ParVector jbv,jxv;
426: PetscCitationsRegister(hypreCitation,&cite);
427: if (!jac->applyrichardson) VecSet(x,0.0);
428: VecHYPRE_IJVectorPushVecRead(hjac->b,b);
429: if (jac->applyrichardson) VecHYPRE_IJVectorPushVec(hjac->x,x);
430: else VecHYPRE_IJVectorPushVecWrite(hjac->x,x);
431: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat);
432: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&jbv);
433: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&jxv);
434: PetscStackCall("Hypre solve",do {
435: HYPRE_Int h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
436: if (hierr) {
438: hypre__global_error = 0;
439: }
440: } while (0));
442: if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
443: PetscStackCallStandard(HYPRE_AMSProjectOutGradients,jac->hsolver,jxv);
444: }
445: VecHYPRE_IJVectorPopVec(hjac->x);
446: VecHYPRE_IJVectorPopVec(hjac->b);
447: return 0;
448: }
450: static PetscErrorCode PCReset_HYPRE(PC pc)
451: {
452: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
454: MatDestroy(&jac->hpmat);
455: MatDestroy(&jac->G);
456: MatDestroy(&jac->C);
457: MatDestroy(&jac->alpha_Poisson);
458: MatDestroy(&jac->beta_Poisson);
459: MatDestroy(&jac->RT_PiFull);
460: MatDestroy(&jac->RT_Pi[0]);
461: MatDestroy(&jac->RT_Pi[1]);
462: MatDestroy(&jac->RT_Pi[2]);
463: MatDestroy(&jac->ND_PiFull);
464: MatDestroy(&jac->ND_Pi[0]);
465: MatDestroy(&jac->ND_Pi[1]);
466: MatDestroy(&jac->ND_Pi[2]);
467: VecHYPRE_IJVectorDestroy(&jac->coords[0]);
468: VecHYPRE_IJVectorDestroy(&jac->coords[1]);
469: VecHYPRE_IJVectorDestroy(&jac->coords[2]);
470: VecHYPRE_IJVectorDestroy(&jac->constants[0]);
471: VecHYPRE_IJVectorDestroy(&jac->constants[1]);
472: VecHYPRE_IJVectorDestroy(&jac->constants[2]);
473: PCHYPREResetNearNullSpace_Private(pc);
474: jac->ams_beta_is_zero = PETSC_FALSE;
475: jac->dim = 0;
476: return 0;
477: }
479: static PetscErrorCode PCDestroy_HYPRE(PC pc)
480: {
481: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
483: PCReset_HYPRE(pc);
484: if (jac->destroy) PetscStackCallStandard(jac->destroy,jac->hsolver);
485: PetscFree(jac->hypre_type);
486: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
487: PetscFree(jac->spgemm_type);
488: #endif
489: if (jac->comm_hypre != MPI_COMM_NULL) PetscCommRestoreComm(PetscObjectComm((PetscObject)pc),&jac->comm_hypre);
490: PetscFree(pc->data);
492: PetscObjectChangeTypeName((PetscObject)pc,0);
493: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
494: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
495: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
496: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
497: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
498: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);
499: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
500: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
501: PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);
502: PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);
503: PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinSetMatProductAlgorithm_C",NULL);
504: PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinGetMatProductAlgorithm_C",NULL);
506: return 0;
507: }
509: /* --------------------------------------------------------------------------------------------*/
510: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
511: {
512: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
513: PetscBool flag;
515: PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
516: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
517: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,jac->hsolver,jac->maxiter);
518: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
519: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,jac->hsolver,jac->tol);
520: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
521: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,jac->hsolver,jac->factorrowsize);
522: PetscOptionsTail();
523: return 0;
524: }
526: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
527: {
528: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
529: PetscBool iascii;
531: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
532: if (iascii) {
533: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
534: if (jac->maxiter != PETSC_DEFAULT) {
535: PetscViewerASCIIPrintf(viewer," maximum number of iterations %d\n",jac->maxiter);
536: } else {
537: PetscViewerASCIIPrintf(viewer," default maximum number of iterations \n");
538: }
539: if (jac->tol != PETSC_DEFAULT) {
540: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)jac->tol);
541: } else {
542: PetscViewerASCIIPrintf(viewer," default drop tolerance \n");
543: }
544: if (jac->factorrowsize != PETSC_DEFAULT) {
545: PetscViewerASCIIPrintf(viewer," factor row size %d\n",jac->factorrowsize);
546: } else {
547: PetscViewerASCIIPrintf(viewer," default factor row size \n");
548: }
549: }
550: return 0;
551: }
553: /* --------------------------------------------------------------------------------------------*/
554: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
555: {
556: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
557: PetscBool flag,eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE;
559: PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");
560: PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);
561: if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,jac->hsolver,jac->eu_level);
563: PetscOptionsReal("-pc_hypre_euclid_droptolerance","Drop tolerance for ILU(k) in Euclid","None",jac->eu_droptolerance,&jac->eu_droptolerance,&flag);
564: if (flag) {
565: PetscMPIInt size;
567: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
569: PetscStackCallStandard(HYPRE_EuclidSetILUT,jac->hsolver,jac->eu_droptolerance);
570: }
572: PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj,&eu_bj,&flag);
573: if (flag) {
574: jac->eu_bj = eu_bj ? 1 : 0;
575: PetscStackCallStandard(HYPRE_EuclidSetBJ,jac->hsolver,jac->eu_bj);
576: }
577: PetscOptionsTail();
578: return 0;
579: }
581: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
582: {
583: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
584: PetscBool iascii;
586: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
587: if (iascii) {
588: PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");
589: if (jac->eu_level != PETSC_DEFAULT) {
590: PetscViewerASCIIPrintf(viewer," factorization levels %d\n",jac->eu_level);
591: } else {
592: PetscViewerASCIIPrintf(viewer," default factorization levels \n");
593: }
594: PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)jac->eu_droptolerance);
595: PetscViewerASCIIPrintf(viewer," use Block-Jacobi? %D\n",jac->eu_bj);
596: }
597: return 0;
598: }
600: /* --------------------------------------------------------------------------------------------*/
602: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
603: {
604: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
605: Mat_HYPRE *hjac = (Mat_HYPRE*)(jac->hpmat->data);
606: HYPRE_ParCSRMatrix hmat;
607: HYPRE_ParVector jbv,jxv;
609: PetscCitationsRegister(hypreCitation,&cite);
610: VecSet(x,0.0);
611: VecHYPRE_IJVectorPushVecRead(hjac->x,b);
612: VecHYPRE_IJVectorPushVecWrite(hjac->b,x);
614: PetscStackCallStandard(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat);
615: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&jbv);
616: PetscStackCallStandard(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&jxv);
618: PetscStackCall("Hypre Transpose solve",do {
619: HYPRE_Int hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
620: if (hierr) {
621: /* error code of 1 in BoomerAMG merely means convergence not achieved */
623: hypre__global_error = 0;
624: }
625: } while (0));
627: VecHYPRE_IJVectorPopVec(hjac->x);
628: VecHYPRE_IJVectorPopVec(hjac->b);
629: return 0;
630: }
632: /* static array length */
633: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
635: static PetscErrorCode PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc,const char name[])
636: {
637: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
638: PetscBool flag;
640: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
641: if (jac->spgemm_type) {
642: PetscStrcmp(jac->spgemm_type,name,&flag);
644: return 0;
645: } else {
646: PetscStrallocpy(name, &jac->spgemm_type);
647: }
648: PetscStrcmp("cusparse",jac->spgemm_type,&flag);
649: if (flag) {
650: PetscStackCallStandard(HYPRE_SetSpGemmUseCusparse,1);
651: return 0;
652: }
653: PetscStrcmp("hypre",jac->spgemm_type,&flag);
654: if (flag) {
655: PetscStackCallStandard(HYPRE_SetSpGemmUseCusparse,0);
656: return 0;
657: }
658: jac->spgemm_type = NULL;
659: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE SpGEM type %s; Choices are cusparse, hypre",name);
660: #endif
661: }
663: static PetscErrorCode PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char *spgemm[])
664: {
665: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
668: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
669: *spgemm = jac->spgemm_type;
670: #endif
671: return 0;
672: }
674: static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
675: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
676: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
677: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
678: static const char *HYPREBoomerAMGSmoothType[] = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
679: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
680: "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
681: "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
682: "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
683: "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
684: static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
685: "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1",
686: "ext", "ad-wts", "ext-mm", "ext+i-mm", "ext+e-mm"};
687: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
688: {
689: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
690: PetscInt bs,n,indx,level;
691: PetscBool flg, tmp_truth;
692: double tmpdbl, twodbl[2];
693: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
694: const char *PCHYPRESpgemmTypes[] = {"cusparse","hypre"};
696: PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
697: PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
698: if (flg) {
699: jac->cycletype = indx+1;
700: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,jac->hsolver,jac->cycletype);
701: }
702: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
703: if (flg) {
705: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,jac->hsolver,jac->maxlevels);
706: }
707: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
708: if (flg) {
710: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,jac->maxiter);
711: }
712: PetscOptionsReal("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
713: if (flg) {
715: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,jac->hsolver,jac->tol);
716: }
717: bs = 1;
718: if (pc->pmat) {
719: MatGetBlockSize(pc->pmat,&bs);
720: }
721: PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);
722: if (flg) {
723: PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,jac->hsolver,bs);
724: }
726: PetscOptionsReal("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
727: if (flg) {
729: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,jac->hsolver,jac->truncfactor);
730: }
732: PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
733: if (flg) {
735: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,jac->hsolver,jac->pmax);
736: }
738: PetscOptionsRangeInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg,0,jac->maxlevels);
739: if (flg) PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,jac->hsolver,jac->agg_nl);
741: PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
742: if (flg) {
744: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,jac->hsolver,jac->agg_num_paths);
745: }
747: PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
748: if (flg) {
750: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,jac->hsolver,jac->strongthreshold);
751: }
752: PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
753: if (flg) {
756: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,jac->hsolver,jac->maxrowsum);
757: }
759: /* Grid sweeps */
760: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
761: if (flg) {
762: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,jac->hsolver,indx);
763: /* modify the jac structure so we can view the updated options with PC_View */
764: jac->gridsweeps[0] = indx;
765: jac->gridsweeps[1] = indx;
766: /*defaults coarse to 1 */
767: jac->gridsweeps[2] = 1;
768: }
769: PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);
770: if (flg) {
771: PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,jac->hsolver,jac->nodal_coarsening);
772: }
773: PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen_diag","Diagonal in strength matrix for nodal based coarsening 0-2","HYPRE_BoomerAMGSetNodalDiag",jac->nodal_coarsening_diag,&jac->nodal_coarsening_diag,&flg);
774: if (flg) {
775: PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,jac->hsolver,jac->nodal_coarsening_diag);
776: }
777: PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
778: if (flg) {
779: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,jac->hsolver,jac->vec_interp_variant);
780: }
781: PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_qmax","Max elements per row for each Q","HYPRE_BoomerAMGSetInterpVecQMax",jac->vec_interp_qmax, &jac->vec_interp_qmax,&flg);
782: if (flg) {
783: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,jac->hsolver,jac->vec_interp_qmax);
784: }
785: PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);
786: if (flg) {
787: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,jac->hsolver,jac->vec_interp_smooth);
788: }
789: PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);
790: if (flg) {
791: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,jac->hsolver,jac->interp_refine);
792: }
793: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
794: if (flg) {
795: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,jac->hsolver,indx, 1);
796: jac->gridsweeps[0] = indx;
797: }
798: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
799: if (flg) {
800: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,jac->hsolver,indx, 2);
801: jac->gridsweeps[1] = indx;
802: }
803: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
804: if (flg) {
805: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,jac->hsolver,indx, 3);
806: jac->gridsweeps[2] = indx;
807: }
809: /* Smooth type */
810: PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
811: if (flg) {
812: jac->smoothtype = indx;
813: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,jac->hsolver,indx+6);
814: jac->smoothnumlevels = 25;
815: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,jac->hsolver,25);
816: }
818: /* Number of smoothing levels */
819: PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
820: if (flg && (jac->smoothtype != -1)) {
821: jac->smoothnumlevels = indx;
822: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,jac->hsolver,indx);
823: }
825: /* Number of levels for ILU(k) for Euclid */
826: PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
827: if (flg && (jac->smoothtype == 3)) {
828: jac->eu_level = indx;
829: PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,jac->hsolver,indx);
830: }
832: /* Filter for ILU(k) for Euclid */
833: double droptolerance;
834: PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
835: if (flg && (jac->smoothtype == 3)) {
836: jac->eu_droptolerance = droptolerance;
837: PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,jac->hsolver,droptolerance);
838: }
840: /* Use Block Jacobi ILUT for Euclid */
841: PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
842: if (flg && (jac->smoothtype == 3)) {
843: jac->eu_bj = tmp_truth;
844: PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,jac->hsolver,jac->eu_bj);
845: }
847: /* Relax type */
848: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
849: if (flg) {
850: jac->relaxtype[0] = jac->relaxtype[1] = indx;
851: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,jac->hsolver, indx);
852: /* by default, coarse type set to 9 */
853: jac->relaxtype[2] = 9;
854: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,jac->hsolver, 9, 3);
855: }
856: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
857: if (flg) {
858: jac->relaxtype[0] = indx;
859: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,jac->hsolver, indx, 1);
860: }
861: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
862: if (flg) {
863: jac->relaxtype[1] = indx;
864: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,jac->hsolver, indx, 2);
865: }
866: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
867: if (flg) {
868: jac->relaxtype[2] = indx;
869: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,jac->hsolver, indx, 3);
870: }
872: /* Relaxation Weight */
873: PetscOptionsReal("-pc_hypre_boomeramg_relax_weight_all","Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)","None",jac->relaxweight, &tmpdbl,&flg);
874: if (flg) {
875: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,jac->hsolver,tmpdbl);
876: jac->relaxweight = tmpdbl;
877: }
879: n = 2;
880: twodbl[0] = twodbl[1] = 1.0;
881: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
882: if (flg) {
883: if (n == 2) {
884: indx = (int)PetscAbsReal(twodbl[1]);
885: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,jac->hsolver,twodbl[0],indx);
886: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
887: }
889: /* Outer relaxation Weight */
890: PetscOptionsReal("-pc_hypre_boomeramg_outer_relax_weight_all","Outer relaxation weight for all levels (-k = determined with k CG steps)","None",jac->outerrelaxweight, &tmpdbl,&flg);
891: if (flg) {
892: PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,jac->hsolver, tmpdbl);
893: jac->outerrelaxweight = tmpdbl;
894: }
896: n = 2;
897: twodbl[0] = twodbl[1] = 1.0;
898: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
899: if (flg) {
900: if (n == 2) {
901: indx = (int)PetscAbsReal(twodbl[1]);
902: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,jac->hsolver, twodbl[0], indx);
903: } else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
904: }
906: /* the Relax Order */
907: PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
909: if (flg && tmp_truth) {
910: jac->relaxorder = 0;
911: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,jac->hsolver, jac->relaxorder);
912: }
913: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
914: if (flg) {
915: jac->measuretype = indx;
916: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,jac->hsolver,jac->measuretype);
917: }
918: /* update list length 3/07 */
919: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
920: if (flg) {
921: jac->coarsentype = indx;
922: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,jac->hsolver,jac->coarsentype);
923: }
925: PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg);
926: if (flg) {
927: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,jac->hsolver, jac->maxc);
928: }
929: PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg);
930: if (flg) {
931: PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,jac->hsolver, jac->minc);
932: }
933: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
934: // global parameter but is closely associated with BoomerAMG
935: PetscOptionsEList("-pc_mg_galerkin_mat_product_algorithm","Type of SpGEMM to use in hypre (only for now)","PCMGGalerkinSetMatProductAlgorithm",PCHYPRESpgemmTypes,ALEN(PCHYPRESpgemmTypes),PCHYPRESpgemmTypes[0],&indx,&flg);
936: if (!flg) indx = 0;
937: PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc,PCHYPRESpgemmTypes[indx]);
938: #endif
939: /* AIR */
940: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
941: PetscOptionsInt("-pc_hypre_boomeramg_restriction_type", "Type of AIR method (distance 1 or 2, 0 means no AIR)", "None", jac->Rtype, &jac->Rtype, NULL);
942: PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,jac->hsolver,jac->Rtype);
943: if (jac->Rtype) {
944: jac->interptype = 100; /* no way we can pass this with strings... Set it as default as in MFEM, then users can still customize it back to a different one */
946: PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR","Threshold for R","None",jac->Rstrongthreshold,&jac->Rstrongthreshold,NULL);
947: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,jac->hsolver,jac->Rstrongthreshold);
949: PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR","Filter threshold for R","None",jac->Rfilterthreshold,&jac->Rfilterthreshold,NULL);
950: PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,jac->hsolver,jac->Rfilterthreshold);
952: PetscOptionsReal("-pc_hypre_boomeramg_Adroptol","Defines the drop tolerance for the A-matrices from the 2nd level of AMG","None",jac->Adroptol,&jac->Adroptol,NULL);
953: PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,jac->hsolver,jac->Adroptol);
955: PetscOptionsInt("-pc_hypre_boomeramg_Adroptype","Drops the entries that are not on the diagonal and smaller than its row norm: type 1: 1-norm, 2: 2-norm, -1: infinity norm","None",jac->Adroptype,&jac->Adroptype,NULL);
956: PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,jac->hsolver,jac->Adroptype);
957: }
958: #endif
960: #if PETSC_PKG_HYPRE_VERSION_LE(9,9,9)
962: #endif
964: /* new 3/07 */
965: PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
966: if (flg || jac->Rtype) {
967: if (flg) jac->interptype = indx;
968: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,jac->hsolver,jac->interptype);
969: }
971: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
972: if (flg) {
973: level = 3;
974: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);
976: jac->printstatistics = PETSC_TRUE;
977: PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,jac->hsolver,level);
978: }
980: PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
981: if (flg) {
982: level = 3;
983: PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);
985: jac->printstatistics = PETSC_TRUE;
986: PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,jac->hsolver,level);
987: }
989: PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
990: if (flg && tmp_truth) {
991: PetscInt tmp_int;
992: PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
993: if (flg) jac->nodal_relax_levels = tmp_int;
994: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,jac->hsolver,6);
995: PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,jac->hsolver,1);
996: PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,jac->hsolver,0);
997: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,jac->hsolver,jac->nodal_relax_levels);
998: }
1000: PetscOptionsBool("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL);
1001: PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,jac->hsolver,jac->keeptranspose ? 1 : 0);
1003: /* options for ParaSails solvers */
1004: PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flg);
1005: if (flg) {
1006: jac->symt = indx;
1007: PetscStackCallStandard(HYPRE_BoomerAMGSetSym,jac->hsolver,jac->symt);
1008: }
1010: PetscOptionsTail();
1011: return 0;
1012: }
1014: static PetscErrorCode PCApplyRichardson_HYPRE_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
1015: {
1016: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1017: HYPRE_Int oits;
1019: PetscCitationsRegister(hypreCitation,&cite);
1020: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,its*jac->maxiter);
1021: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,jac->hsolver,rtol);
1022: jac->applyrichardson = PETSC_TRUE;
1023: PCApply_HYPRE(pc,b,y);
1024: jac->applyrichardson = PETSC_FALSE;
1025: PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,jac->hsolver,&oits);
1026: *outits = oits;
1027: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1028: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1029: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,jac->hsolver,jac->tol);
1030: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,jac->maxiter);
1031: return 0;
1032: }
1034: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
1035: {
1036: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1037: PetscBool iascii;
1039: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1040: if (iascii) {
1041: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
1042: PetscViewerASCIIPrintf(viewer," Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
1043: PetscViewerASCIIPrintf(viewer," Maximum number of levels %D\n",jac->maxlevels);
1044: PetscViewerASCIIPrintf(viewer," Maximum number of iterations PER hypre call %D\n",jac->maxiter);
1045: PetscViewerASCIIPrintf(viewer," Convergence tolerance PER hypre call %g\n",(double)jac->tol);
1046: PetscViewerASCIIPrintf(viewer," Threshold for strong coupling %g\n",(double)jac->strongthreshold);
1047: PetscViewerASCIIPrintf(viewer," Interpolation truncation factor %g\n",(double)jac->truncfactor);
1048: PetscViewerASCIIPrintf(viewer," Interpolation: max elements per row %D\n",jac->pmax);
1049: if (jac->interp_refine) {
1050: PetscViewerASCIIPrintf(viewer," Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);
1051: }
1052: PetscViewerASCIIPrintf(viewer," Number of levels of aggressive coarsening %D\n",jac->agg_nl);
1053: PetscViewerASCIIPrintf(viewer," Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);
1055: PetscViewerASCIIPrintf(viewer," Maximum row sums %g\n",(double)jac->maxrowsum);
1057: PetscViewerASCIIPrintf(viewer," Sweeps down %D\n",jac->gridsweeps[0]);
1058: PetscViewerASCIIPrintf(viewer," Sweeps up %D\n",jac->gridsweeps[1]);
1059: PetscViewerASCIIPrintf(viewer," Sweeps on coarse %D\n",jac->gridsweeps[2]);
1061: PetscViewerASCIIPrintf(viewer," Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
1062: PetscViewerASCIIPrintf(viewer," Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
1063: PetscViewerASCIIPrintf(viewer," Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
1065: PetscViewerASCIIPrintf(viewer," Relax weight (all) %g\n",(double)jac->relaxweight);
1066: PetscViewerASCIIPrintf(viewer," Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);
1068: if (jac->relaxorder) {
1069: PetscViewerASCIIPrintf(viewer," Using CF-relaxation\n");
1070: } else {
1071: PetscViewerASCIIPrintf(viewer," Not using CF-relaxation\n");
1072: }
1073: if (jac->smoothtype!=-1) {
1074: PetscViewerASCIIPrintf(viewer," Smooth type %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
1075: PetscViewerASCIIPrintf(viewer," Smooth num levels %D\n",jac->smoothnumlevels);
1076: } else {
1077: PetscViewerASCIIPrintf(viewer," Not using more complex smoothers.\n");
1078: }
1079: if (jac->smoothtype==3) {
1080: PetscViewerASCIIPrintf(viewer," Euclid ILU(k) levels %D\n",jac->eu_level);
1081: PetscViewerASCIIPrintf(viewer," Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);
1082: PetscViewerASCIIPrintf(viewer," Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);
1083: }
1084: PetscViewerASCIIPrintf(viewer," Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
1085: PetscViewerASCIIPrintf(viewer," Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
1086: PetscViewerASCIIPrintf(viewer," Interpolation type %s\n",jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt");
1087: if (jac->nodal_coarsening) {
1088: PetscViewerASCIIPrintf(viewer," Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
1089: }
1090: if (jac->vec_interp_variant) {
1091: PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
1092: PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);
1093: PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);
1094: }
1095: if (jac->nodal_relax) {
1096: PetscViewerASCIIPrintf(viewer," Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);
1097: }
1098: #if PETSC_PKG_HYPRE_VERSION_GE(2,23,0)
1099: PetscViewerASCIIPrintf(viewer," SpGEMM type %s\n",jac->spgemm_type);
1100: #endif
1101: /* AIR */
1102: if (jac->Rtype) {
1103: PetscViewerASCIIPrintf(viewer," Using approximate ideal restriction type %D\n",jac->Rtype);
1104: PetscViewerASCIIPrintf(viewer," Threshold for R %g\n",(double)jac->Rstrongthreshold);
1105: PetscViewerASCIIPrintf(viewer," Filter for R %g\n",(double)jac->Rfilterthreshold);
1106: PetscViewerASCIIPrintf(viewer," A drop tolerance %g\n",(double)jac->Adroptol);
1107: PetscViewerASCIIPrintf(viewer," A drop type %D\n",jac->Adroptype);
1108: }
1109: }
1110: return 0;
1111: }
1113: /* --------------------------------------------------------------------------------------------*/
1114: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
1115: {
1116: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1117: PetscInt indx;
1118: PetscBool flag;
1119: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
1121: PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
1122: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
1123: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshold,&jac->threshold,&flag);
1124: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,jac->hsolver,jac->threshold,jac->nlevels);
1126: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
1127: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,jac->hsolver,jac->filter);
1129: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
1130: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,jac->hsolver,jac->loadbal);
1132: PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
1133: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,jac->hsolver,jac->logging);
1135: PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
1136: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,jac->hsolver,jac->ruse);
1138: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
1139: if (flag) {
1140: jac->symt = indx;
1141: PetscStackCallStandard(HYPRE_ParaSailsSetSym,jac->hsolver,jac->symt);
1142: }
1144: PetscOptionsTail();
1145: return 0;
1146: }
1148: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
1149: {
1150: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1151: PetscBool iascii;
1152: const char *symt = 0;
1154: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1155: if (iascii) {
1156: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
1157: PetscViewerASCIIPrintf(viewer," nlevels %d\n",jac->nlevels);
1158: PetscViewerASCIIPrintf(viewer," threshold %g\n",(double)jac->threshold);
1159: PetscViewerASCIIPrintf(viewer," filter %g\n",(double)jac->filter);
1160: PetscViewerASCIIPrintf(viewer," load balance %g\n",(double)jac->loadbal);
1161: PetscViewerASCIIPrintf(viewer," reuse nonzero structure %s\n",PetscBools[jac->ruse]);
1162: PetscViewerASCIIPrintf(viewer," print info to screen %s\n",PetscBools[jac->logging]);
1163: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1164: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1165: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1166: else SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1167: PetscViewerASCIIPrintf(viewer," %s\n",symt);
1168: }
1169: return 0;
1170: }
1171: /* --------------------------------------------------------------------------------------------*/
1172: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
1173: {
1174: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1175: PetscInt n;
1176: PetscBool flag,flag2,flag3,flag4;
1178: PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
1179: PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
1180: if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,jac->hsolver,jac->as_print);
1181: PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1182: if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,jac->hsolver,jac->as_max_iter);
1183: PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
1184: if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,jac->hsolver,jac->ams_cycle_type);
1185: PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1186: if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,jac->hsolver,jac->as_tol);
1187: PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1188: PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1189: PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1190: PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1191: if (flag || flag2 || flag3 || flag4) {
1192: PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,jac->hsolver,jac->as_relax_type,
1193: jac->as_relax_times,
1194: jac->as_relax_weight,
1195: jac->as_omega);
1196: }
1197: PetscOptionsReal("-pc_hypre_ams_amg_alpha_theta","Threshold for strong coupling of vector Poisson AMG solver","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);
1198: n = 5;
1199: PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
1200: if (flag || flag2) {
1201: PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1202: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1203: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1204: jac->as_amg_alpha_theta,
1205: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1206: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
1207: }
1208: PetscOptionsReal("-pc_hypre_ams_amg_beta_theta","Threshold for strong coupling of scalar Poisson AMG solver","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);
1209: n = 5;
1210: PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
1211: if (flag || flag2) {
1212: PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1213: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1214: jac->as_amg_beta_opts[2], /* AMG relax_type */
1215: jac->as_amg_beta_theta,
1216: jac->as_amg_beta_opts[3], /* AMG interp_type */
1217: jac->as_amg_beta_opts[4]); /* AMG Pmax */
1218: }
1219: PetscOptionsInt("-pc_hypre_ams_projection_frequency","Frequency at which a projection onto the compatible subspace for problems with zero conductivity regions is performed","None",jac->ams_proj_freq,&jac->ams_proj_freq,&flag);
1220: if (flag) { /* override HYPRE's default only if the options is used */
1221: PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,jac->hsolver,jac->ams_proj_freq);
1222: }
1223: PetscOptionsTail();
1224: return 0;
1225: }
1227: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
1228: {
1229: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1230: PetscBool iascii;
1232: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1233: if (iascii) {
1234: PetscViewerASCIIPrintf(viewer," HYPRE AMS preconditioning\n");
1235: PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);
1236: PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);
1237: PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);
1238: PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);
1239: PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);
1240: PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);
1241: PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);
1242: if (jac->alpha_Poisson) {
1243: PetscViewerASCIIPrintf(viewer," vector Poisson solver (passed in by user)\n");
1244: } else {
1245: PetscViewerASCIIPrintf(viewer," vector Poisson solver (computed) \n");
1246: }
1247: PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1248: PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1249: PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1250: PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1251: PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1252: PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1253: if (!jac->ams_beta_is_zero) {
1254: if (jac->beta_Poisson) {
1255: PetscViewerASCIIPrintf(viewer," scalar Poisson solver (passed in by user)\n");
1256: } else {
1257: PetscViewerASCIIPrintf(viewer," scalar Poisson solver (computed) \n");
1258: }
1259: PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1260: PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1261: PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1262: PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1263: PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1264: PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1265: if (jac->ams_beta_is_zero_part) {
1266: PetscViewerASCIIPrintf(viewer," compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1267: }
1268: } else {
1269: PetscViewerASCIIPrintf(viewer," scalar Poisson solver not used (zero-conductivity everywhere) \n");
1270: }
1271: }
1272: return 0;
1273: }
1275: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1276: {
1277: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1278: PetscInt n;
1279: PetscBool flag,flag2,flag3,flag4;
1281: PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1282: PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1283: if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,jac->hsolver,jac->as_print);
1284: PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1285: if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,jac->hsolver,jac->as_max_iter);
1286: PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1287: if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,jac->hsolver,jac->ads_cycle_type);
1288: PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1289: if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,jac->hsolver,jac->as_tol);
1290: PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1291: PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1292: PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1293: PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1294: if (flag || flag2 || flag3 || flag4) {
1295: PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,jac->hsolver,jac->as_relax_type,
1296: jac->as_relax_times,
1297: jac->as_relax_weight,
1298: jac->as_omega);
1299: }
1300: PetscOptionsReal("-pc_hypre_ads_ams_theta","Threshold for strong coupling of AMS solver inside ADS","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);
1301: n = 5;
1302: PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1303: PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1304: if (flag || flag2 || flag3) {
1305: PetscStackCallStandard(HYPRE_ADSSetAMSOptions,jac->hsolver,jac->ams_cycle_type, /* AMS cycle type */
1306: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1307: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1308: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1309: jac->as_amg_alpha_theta,
1310: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1311: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
1312: }
1313: PetscOptionsReal("-pc_hypre_ads_amg_theta","Threshold for strong coupling of vector AMG solver inside ADS","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);
1314: n = 5;
1315: PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1316: if (flag || flag2) {
1317: PetscStackCallStandard(HYPRE_ADSSetAMGOptions,jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1318: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1319: jac->as_amg_beta_opts[2], /* AMG relax_type */
1320: jac->as_amg_beta_theta,
1321: jac->as_amg_beta_opts[3], /* AMG interp_type */
1322: jac->as_amg_beta_opts[4]); /* AMG Pmax */
1323: }
1324: PetscOptionsTail();
1325: return 0;
1326: }
1328: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1329: {
1330: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1331: PetscBool iascii;
1333: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1334: if (iascii) {
1335: PetscViewerASCIIPrintf(viewer," HYPRE ADS preconditioning\n");
1336: PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);
1337: PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ads_cycle_type);
1338: PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);
1339: PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);
1340: PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);
1341: PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);
1342: PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);
1343: PetscViewerASCIIPrintf(viewer," AMS solver using boomerAMG\n");
1344: PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);
1345: PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1346: PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1347: PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1348: PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1349: PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1350: PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_alpha_theta);
1351: PetscViewerASCIIPrintf(viewer," vector Poisson solver using boomerAMG\n");
1352: PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_beta_opts[0]);
1353: PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1354: PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_beta_opts[2]);
1355: PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_beta_opts[3]);
1356: PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1357: PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_beta_theta);
1358: }
1359: return 0;
1360: }
1362: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1363: {
1364: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1365: PetscBool ishypre;
1367: PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1368: if (ishypre) {
1369: PetscObjectReference((PetscObject)G);
1370: MatDestroy(&jac->G);
1371: jac->G = G;
1372: } else {
1373: MatDestroy(&jac->G);
1374: MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);
1375: }
1376: return 0;
1377: }
1379: /*@
1380: PCHYPRESetDiscreteGradient - Set discrete gradient matrix
1382: Collective on PC
1384: Input Parameters:
1385: + pc - the preconditioning context
1386: - G - the discrete gradient
1388: Level: intermediate
1390: Notes:
1391: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1393: Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix entries are +1 and -1 depending on edge orientation
1395: .seealso: PCHYPRESetDiscreteCurl()
1396: @*/
1397: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1398: {
1402: PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1403: return 0;
1404: }
1406: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1407: {
1408: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1409: PetscBool ishypre;
1411: PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1412: if (ishypre) {
1413: PetscObjectReference((PetscObject)C);
1414: MatDestroy(&jac->C);
1415: jac->C = C;
1416: } else {
1417: MatDestroy(&jac->C);
1418: MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);
1419: }
1420: return 0;
1421: }
1423: /*@
1424: PCHYPRESetDiscreteCurl - Set discrete curl matrix
1426: Collective on PC
1428: Input Parameters:
1429: + pc - the preconditioning context
1430: - C - the discrete curl
1432: Level: intermediate
1434: Notes:
1435: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1437: Each row of G has as many nonzeros as the number of edges of a face, with column indexes being the global indexes of the corresponding edge: matrix entries are +1 and -1 depending on edge orientation with respect to the face orientation
1439: .seealso: PCHYPRESetDiscreteGradient()
1440: @*/
1441: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1442: {
1446: PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1447: return 0;
1448: }
1450: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1451: {
1452: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1453: PetscBool ishypre;
1454: PetscInt i;
1456: MatDestroy(&jac->RT_PiFull);
1457: MatDestroy(&jac->ND_PiFull);
1458: for (i=0;i<3;++i) {
1459: MatDestroy(&jac->RT_Pi[i]);
1460: MatDestroy(&jac->ND_Pi[i]);
1461: }
1463: jac->dim = dim;
1464: if (RT_PiFull) {
1465: PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);
1466: if (ishypre) {
1467: PetscObjectReference((PetscObject)RT_PiFull);
1468: jac->RT_PiFull = RT_PiFull;
1469: } else {
1470: MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);
1471: }
1472: }
1473: if (RT_Pi) {
1474: for (i=0;i<dim;++i) {
1475: if (RT_Pi[i]) {
1476: PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);
1477: if (ishypre) {
1478: PetscObjectReference((PetscObject)RT_Pi[i]);
1479: jac->RT_Pi[i] = RT_Pi[i];
1480: } else {
1481: MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);
1482: }
1483: }
1484: }
1485: }
1486: if (ND_PiFull) {
1487: PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);
1488: if (ishypre) {
1489: PetscObjectReference((PetscObject)ND_PiFull);
1490: jac->ND_PiFull = ND_PiFull;
1491: } else {
1492: MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);
1493: }
1494: }
1495: if (ND_Pi) {
1496: for (i=0;i<dim;++i) {
1497: if (ND_Pi[i]) {
1498: PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);
1499: if (ishypre) {
1500: PetscObjectReference((PetscObject)ND_Pi[i]);
1501: jac->ND_Pi[i] = ND_Pi[i];
1502: } else {
1503: MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);
1504: }
1505: }
1506: }
1507: }
1509: return 0;
1510: }
1512: /*@
1513: PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner
1515: Collective on PC
1517: Input Parameters:
1518: + pc - the preconditioning context
1519: - dim - the dimension of the problem, only used in AMS
1520: - RT_PiFull - Raviart-Thomas interpolation matrix
1521: - RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1522: - ND_PiFull - Nedelec interpolation matrix
1523: - ND_Pi - x/y/z component of Nedelec interpolation matrix
1525: Notes:
1526: For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1528: For ADS, both type of interpolation matrices are needed.
1530: Level: intermediate
1532: @*/
1533: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1534: {
1535: PetscInt i;
1538: if (RT_PiFull) {
1541: }
1542: if (RT_Pi) {
1544: for (i=0;i<dim;++i) {
1545: if (RT_Pi[i]) {
1548: }
1549: }
1550: }
1551: if (ND_PiFull) {
1554: }
1555: if (ND_Pi) {
1557: for (i=0;i<dim;++i) {
1558: if (ND_Pi[i]) {
1561: }
1562: }
1563: }
1564: PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));
1565: return 0;
1566: }
1568: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1569: {
1570: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1571: PetscBool ishypre;
1573: PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1574: if (ishypre) {
1575: if (isalpha) {
1576: PetscObjectReference((PetscObject)A);
1577: MatDestroy(&jac->alpha_Poisson);
1578: jac->alpha_Poisson = A;
1579: } else {
1580: if (A) {
1581: PetscObjectReference((PetscObject)A);
1582: } else {
1583: jac->ams_beta_is_zero = PETSC_TRUE;
1584: }
1585: MatDestroy(&jac->beta_Poisson);
1586: jac->beta_Poisson = A;
1587: }
1588: } else {
1589: if (isalpha) {
1590: MatDestroy(&jac->alpha_Poisson);
1591: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);
1592: } else {
1593: if (A) {
1594: MatDestroy(&jac->beta_Poisson);
1595: MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);
1596: } else {
1597: MatDestroy(&jac->beta_Poisson);
1598: jac->ams_beta_is_zero = PETSC_TRUE;
1599: }
1600: }
1601: }
1602: return 0;
1603: }
1605: /*@
1606: PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix
1608: Collective on PC
1610: Input Parameters:
1611: + pc - the preconditioning context
1612: - A - the matrix
1614: Level: intermediate
1616: Notes:
1617: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements
1619: .seealso: PCHYPRESetDiscreteGradient(), PCHYPRESetDiscreteCurl(), PCHYPRESetBetaPoissonMatrix()
1620: @*/
1621: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1622: {
1626: PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1627: return 0;
1628: }
1630: /*@
1631: PCHYPRESetBetaPoissonMatrix - Set Poisson matrix
1633: Collective on PC
1635: Input Parameters:
1636: + pc - the preconditioning context
1637: - A - the matrix
1639: Level: intermediate
1641: Notes:
1642: A should be obtained by discretizing the Poisson problem with linear finite elements.
1643: Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.
1645: .seealso: PCHYPRESetDiscreteGradient(), PCHYPRESetDiscreteCurl(), PCHYPRESetAlphaPoissonMatrix()
1646: @*/
1647: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1648: {
1650: if (A) {
1653: }
1654: PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1655: return 0;
1656: }
1658: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1659: {
1660: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1662: /* throw away any vector if already set */
1663: VecHYPRE_IJVectorDestroy(&jac->constants[0]);
1664: VecHYPRE_IJVectorDestroy(&jac->constants[1]);
1665: VecHYPRE_IJVectorDestroy(&jac->constants[2]);
1666: VecHYPRE_IJVectorCreate(ozz->map,&jac->constants[0]);
1667: VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1668: VecHYPRE_IJVectorCreate(zoz->map,&jac->constants[1]);
1669: VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1670: jac->dim = 2;
1671: if (zzo) {
1672: VecHYPRE_IJVectorCreate(zzo->map,&jac->constants[2]);
1673: VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1674: jac->dim++;
1675: }
1676: return 0;
1677: }
1679: /*@
1680: PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in the edge element basis
1682: Collective on PC
1684: Input Parameters:
1685: + pc - the preconditioning context
1686: - ozz - vector representing (1,0,0) (or (1,0) in 2D)
1687: - zoz - vector representing (0,1,0) (or (0,1) in 2D)
1688: - zzo - vector representing (0,0,1) (use NULL in 2D)
1690: Level: intermediate
1692: @*/
1693: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1694: {
1702: PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1703: return 0;
1704: }
1706: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1707: {
1708: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1709: Vec tv;
1710: PetscInt i;
1712: /* throw away any coordinate vector if already set */
1713: VecHYPRE_IJVectorDestroy(&jac->coords[0]);
1714: VecHYPRE_IJVectorDestroy(&jac->coords[1]);
1715: VecHYPRE_IJVectorDestroy(&jac->coords[2]);
1716: jac->dim = dim;
1718: /* compute IJ vector for coordinates */
1719: VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1720: VecSetType(tv,VECSTANDARD);
1721: VecSetSizes(tv,nloc,PETSC_DECIDE);
1722: for (i=0;i<dim;i++) {
1723: PetscScalar *array;
1724: PetscInt j;
1726: VecHYPRE_IJVectorCreate(tv->map,&jac->coords[i]);
1727: VecGetArrayWrite(tv,&array);
1728: for (j=0;j<nloc;j++) array[j] = coords[j*dim+i];
1729: VecRestoreArrayWrite(tv,&array);
1730: VecHYPRE_IJVectorCopy(tv,jac->coords[i]);
1731: }
1732: VecDestroy(&tv);
1733: return 0;
1734: }
1736: /* ---------------------------------------------------------------------------------*/
1738: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
1739: {
1740: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1742: *name = jac->hypre_type;
1743: return 0;
1744: }
1746: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
1747: {
1748: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1749: PetscBool flag;
1751: if (jac->hypre_type) {
1752: PetscStrcmp(jac->hypre_type,name,&flag);
1754: return 0;
1755: } else {
1756: PetscStrallocpy(name, &jac->hypre_type);
1757: }
1759: jac->maxiter = PETSC_DEFAULT;
1760: jac->tol = PETSC_DEFAULT;
1761: jac->printstatistics = PetscLogPrintInfo;
1763: PetscStrcmp("pilut",jac->hypre_type,&flag);
1764: if (flag) {
1765: PetscCommGetComm(PetscObjectComm((PetscObject)pc),&jac->comm_hypre);
1766: PetscStackCallStandard(HYPRE_ParCSRPilutCreate,jac->comm_hypre,&jac->hsolver);
1767: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1768: pc->ops->view = PCView_HYPRE_Pilut;
1769: jac->destroy = HYPRE_ParCSRPilutDestroy;
1770: jac->setup = HYPRE_ParCSRPilutSetup;
1771: jac->solve = HYPRE_ParCSRPilutSolve;
1772: jac->factorrowsize = PETSC_DEFAULT;
1773: return 0;
1774: }
1775: PetscStrcmp("euclid",jac->hypre_type,&flag);
1776: if (flag) {
1777: #if defined(PETSC_USE_64BIT_INDICES)
1778: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Hypre Euclid does not support 64 bit indices");
1779: #endif
1780: PetscCommGetComm(PetscObjectComm((PetscObject)pc),&jac->comm_hypre);
1781: PetscStackCallStandard(HYPRE_EuclidCreate,jac->comm_hypre,&jac->hsolver);
1782: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1783: pc->ops->view = PCView_HYPRE_Euclid;
1784: jac->destroy = HYPRE_EuclidDestroy;
1785: jac->setup = HYPRE_EuclidSetup;
1786: jac->solve = HYPRE_EuclidSolve;
1787: jac->factorrowsize = PETSC_DEFAULT;
1788: jac->eu_level = PETSC_DEFAULT; /* default */
1789: return 0;
1790: }
1791: PetscStrcmp("parasails",jac->hypre_type,&flag);
1792: if (flag) {
1793: PetscCommGetComm(PetscObjectComm((PetscObject)pc),&jac->comm_hypre);
1794: PetscStackCallStandard(HYPRE_ParaSailsCreate,jac->comm_hypre,&jac->hsolver);
1795: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1796: pc->ops->view = PCView_HYPRE_ParaSails;
1797: jac->destroy = HYPRE_ParaSailsDestroy;
1798: jac->setup = HYPRE_ParaSailsSetup;
1799: jac->solve = HYPRE_ParaSailsSolve;
1800: /* initialize */
1801: jac->nlevels = 1;
1802: jac->threshold = .1;
1803: jac->filter = .1;
1804: jac->loadbal = 0;
1805: if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1806: else jac->logging = (int) PETSC_FALSE;
1808: jac->ruse = (int) PETSC_FALSE;
1809: jac->symt = 0;
1810: PetscStackCallStandard(HYPRE_ParaSailsSetParams,jac->hsolver,jac->threshold,jac->nlevels);
1811: PetscStackCallStandard(HYPRE_ParaSailsSetFilter,jac->hsolver,jac->filter);
1812: PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,jac->hsolver,jac->loadbal);
1813: PetscStackCallStandard(HYPRE_ParaSailsSetLogging,jac->hsolver,jac->logging);
1814: PetscStackCallStandard(HYPRE_ParaSailsSetReuse,jac->hsolver,jac->ruse);
1815: PetscStackCallStandard(HYPRE_ParaSailsSetSym,jac->hsolver,jac->symt);
1816: return 0;
1817: }
1818: PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1819: if (flag) {
1820: PetscStackCallStandard(HYPRE_BoomerAMGCreate,&jac->hsolver);
1821: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
1822: pc->ops->view = PCView_HYPRE_BoomerAMG;
1823: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
1824: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1825: PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);
1826: PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);
1827: jac->destroy = HYPRE_BoomerAMGDestroy;
1828: jac->setup = HYPRE_BoomerAMGSetup;
1829: jac->solve = HYPRE_BoomerAMGSolve;
1830: jac->applyrichardson = PETSC_FALSE;
1831: /* these defaults match the hypre defaults */
1832: jac->cycletype = 1;
1833: jac->maxlevels = 25;
1834: jac->maxiter = 1;
1835: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1836: jac->truncfactor = 0.0;
1837: jac->strongthreshold = .25;
1838: jac->maxrowsum = .9;
1839: jac->coarsentype = 6;
1840: jac->measuretype = 0;
1841: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1842: jac->smoothtype = -1; /* Not set by default */
1843: jac->smoothnumlevels = 25;
1844: jac->eu_level = 0;
1845: jac->eu_droptolerance = 0;
1846: jac->eu_bj = 0;
1847: jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
1848: jac->relaxtype[2] = 9; /*G.E. */
1849: jac->relaxweight = 1.0;
1850: jac->outerrelaxweight = 1.0;
1851: jac->relaxorder = 1;
1852: jac->interptype = 0;
1853: jac->Rtype = 0;
1854: jac->Rstrongthreshold = 0.25;
1855: jac->Rfilterthreshold = 0.0;
1856: jac->Adroptype = -1;
1857: jac->Adroptol = 0.0;
1858: jac->agg_nl = 0;
1859: jac->agg_interptype = 4;
1860: jac->pmax = 0;
1861: jac->truncfactor = 0.0;
1862: jac->agg_num_paths = 1;
1863: jac->maxc = 9;
1864: jac->minc = 1;
1865: jac->nodal_coarsening = 0;
1866: jac->nodal_coarsening_diag = 0;
1867: jac->vec_interp_variant = 0;
1868: jac->vec_interp_qmax = 0;
1869: jac->vec_interp_smooth = PETSC_FALSE;
1870: jac->interp_refine = 0;
1871: jac->nodal_relax = PETSC_FALSE;
1872: jac->nodal_relax_levels = 1;
1873: jac->rap2 = 0;
1875: /* GPU defaults
1876: from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
1877: and /src/parcsr_ls/par_amg.c */
1878: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1879: jac->keeptranspose = PETSC_TRUE;
1880: jac->mod_rap2 = 1;
1881: jac->coarsentype = 8;
1882: jac->relaxorder = 0;
1883: jac->interptype = 6;
1884: jac->relaxtype[0] = 18;
1885: jac->relaxtype[1] = 18;
1886: jac->agg_interptype = 7;
1887: #else
1888: jac->keeptranspose = PETSC_FALSE;
1889: jac->mod_rap2 = 0;
1890: #endif
1891: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,jac->hsolver,jac->cycletype);
1892: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,jac->hsolver,jac->maxlevels);
1893: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,jac->maxiter);
1894: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,jac->hsolver,jac->tol);
1895: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,jac->hsolver,jac->truncfactor);
1896: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,jac->hsolver,jac->strongthreshold);
1897: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,jac->hsolver,jac->maxrowsum);
1898: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,jac->hsolver,jac->coarsentype);
1899: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,jac->hsolver,jac->measuretype);
1900: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,jac->hsolver, jac->relaxorder);
1901: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,jac->hsolver,jac->interptype);
1902: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,jac->hsolver,jac->agg_nl);
1903: PetscStackCallStandard(HYPRE_BoomerAMGSetAggInterpType,jac->hsolver,jac->agg_interptype);
1904: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,jac->hsolver,jac->pmax);
1905: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,jac->hsolver,jac->agg_num_paths);
1906: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,jac->hsolver, jac->relaxtype[0]); /* defaults coarse to 9 */
1907: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */
1908: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,jac->hsolver, jac->maxc);
1909: PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,jac->hsolver, jac->minc);
1910: /* GPU */
1911: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
1912: PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,jac->hsolver,jac->keeptranspose ? 1 : 0);
1913: PetscStackCallStandard(HYPRE_BoomerAMGSetRAP2,jac->hsolver, jac->rap2);
1914: PetscStackCallStandard(HYPRE_BoomerAMGSetModuleRAP2,jac->hsolver, jac->mod_rap2);
1915: #endif
1917: /* AIR */
1918: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
1919: PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,jac->hsolver,jac->Rtype);
1920: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,jac->hsolver,jac->Rstrongthreshold);
1921: PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,jac->hsolver,jac->Rfilterthreshold);
1922: PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,jac->hsolver,jac->Adroptol);
1923: PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,jac->hsolver,jac->Adroptype);
1924: #endif
1925: return 0;
1926: }
1927: PetscStrcmp("ams",jac->hypre_type,&flag);
1928: if (flag) {
1929: HYPRE_AMSCreate(&jac->hsolver);
1930: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
1931: pc->ops->view = PCView_HYPRE_AMS;
1932: jac->destroy = HYPRE_AMSDestroy;
1933: jac->setup = HYPRE_AMSSetup;
1934: jac->solve = HYPRE_AMSSolve;
1935: jac->coords[0] = NULL;
1936: jac->coords[1] = NULL;
1937: jac->coords[2] = NULL;
1938: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1939: jac->as_print = 0;
1940: jac->as_max_iter = 1; /* used as a preconditioner */
1941: jac->as_tol = 0.; /* used as a preconditioner */
1942: jac->ams_cycle_type = 13;
1943: /* Smoothing options */
1944: jac->as_relax_type = 2;
1945: jac->as_relax_times = 1;
1946: jac->as_relax_weight = 1.0;
1947: jac->as_omega = 1.0;
1948: /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1949: jac->as_amg_alpha_opts[0] = 10;
1950: jac->as_amg_alpha_opts[1] = 1;
1951: jac->as_amg_alpha_opts[2] = 6;
1952: jac->as_amg_alpha_opts[3] = 6;
1953: jac->as_amg_alpha_opts[4] = 4;
1954: jac->as_amg_alpha_theta = 0.25;
1955: /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1956: jac->as_amg_beta_opts[0] = 10;
1957: jac->as_amg_beta_opts[1] = 1;
1958: jac->as_amg_beta_opts[2] = 6;
1959: jac->as_amg_beta_opts[3] = 6;
1960: jac->as_amg_beta_opts[4] = 4;
1961: jac->as_amg_beta_theta = 0.25;
1962: PetscStackCallStandard(HYPRE_AMSSetPrintLevel,jac->hsolver,jac->as_print);
1963: PetscStackCallStandard(HYPRE_AMSSetMaxIter,jac->hsolver,jac->as_max_iter);
1964: PetscStackCallStandard(HYPRE_AMSSetCycleType,jac->hsolver,jac->ams_cycle_type);
1965: PetscStackCallStandard(HYPRE_AMSSetTol,jac->hsolver,jac->as_tol);
1966: PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,jac->hsolver,jac->as_relax_type,
1967: jac->as_relax_times,
1968: jac->as_relax_weight,
1969: jac->as_omega);
1970: PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1971: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1972: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1973: jac->as_amg_alpha_theta,
1974: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1975: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
1976: PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1977: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1978: jac->as_amg_beta_opts[2], /* AMG relax_type */
1979: jac->as_amg_beta_theta,
1980: jac->as_amg_beta_opts[3], /* AMG interp_type */
1981: jac->as_amg_beta_opts[4]); /* AMG Pmax */
1982: /* Zero conductivity */
1983: jac->ams_beta_is_zero = PETSC_FALSE;
1984: jac->ams_beta_is_zero_part = PETSC_FALSE;
1985: return 0;
1986: }
1987: PetscStrcmp("ads",jac->hypre_type,&flag);
1988: if (flag) {
1989: HYPRE_ADSCreate(&jac->hsolver);
1990: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
1991: pc->ops->view = PCView_HYPRE_ADS;
1992: jac->destroy = HYPRE_ADSDestroy;
1993: jac->setup = HYPRE_ADSSetup;
1994: jac->solve = HYPRE_ADSSolve;
1995: jac->coords[0] = NULL;
1996: jac->coords[1] = NULL;
1997: jac->coords[2] = NULL;
1998: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1999: jac->as_print = 0;
2000: jac->as_max_iter = 1; /* used as a preconditioner */
2001: jac->as_tol = 0.; /* used as a preconditioner */
2002: jac->ads_cycle_type = 13;
2003: /* Smoothing options */
2004: jac->as_relax_type = 2;
2005: jac->as_relax_times = 1;
2006: jac->as_relax_weight = 1.0;
2007: jac->as_omega = 1.0;
2008: /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2009: jac->ams_cycle_type = 14;
2010: jac->as_amg_alpha_opts[0] = 10;
2011: jac->as_amg_alpha_opts[1] = 1;
2012: jac->as_amg_alpha_opts[2] = 6;
2013: jac->as_amg_alpha_opts[3] = 6;
2014: jac->as_amg_alpha_opts[4] = 4;
2015: jac->as_amg_alpha_theta = 0.25;
2016: /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2017: jac->as_amg_beta_opts[0] = 10;
2018: jac->as_amg_beta_opts[1] = 1;
2019: jac->as_amg_beta_opts[2] = 6;
2020: jac->as_amg_beta_opts[3] = 6;
2021: jac->as_amg_beta_opts[4] = 4;
2022: jac->as_amg_beta_theta = 0.25;
2023: PetscStackCallStandard(HYPRE_ADSSetPrintLevel,jac->hsolver,jac->as_print);
2024: PetscStackCallStandard(HYPRE_ADSSetMaxIter,jac->hsolver,jac->as_max_iter);
2025: PetscStackCallStandard(HYPRE_ADSSetCycleType,jac->hsolver,jac->ams_cycle_type);
2026: PetscStackCallStandard(HYPRE_ADSSetTol,jac->hsolver,jac->as_tol);
2027: PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,jac->hsolver,jac->as_relax_type,
2028: jac->as_relax_times,
2029: jac->as_relax_weight,
2030: jac->as_omega);
2031: PetscStackCallStandard(HYPRE_ADSSetAMSOptions,jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */
2032: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2033: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
2034: jac->as_amg_alpha_opts[2], /* AMG relax_type */
2035: jac->as_amg_alpha_theta,
2036: jac->as_amg_alpha_opts[3], /* AMG interp_type */
2037: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
2038: PetscStackCallStandard(HYPRE_ADSSetAMGOptions,jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
2039: jac->as_amg_beta_opts[1], /* AMG agg_levels */
2040: jac->as_amg_beta_opts[2], /* AMG relax_type */
2041: jac->as_amg_beta_theta,
2042: jac->as_amg_beta_opts[3], /* AMG interp_type */
2043: jac->as_amg_beta_opts[4]); /* AMG Pmax */
2044: return 0;
2045: }
2046: PetscFree(jac->hypre_type);
2048: jac->hypre_type = NULL;
2049: SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name);
2050: }
2052: /*
2053: It only gets here if the HYPRE type has not been set before the call to
2054: ...SetFromOptions() which actually is most of the time
2055: */
2056: PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
2057: {
2058: PetscInt indx;
2059: const char *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
2060: PetscBool flg;
2062: PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
2063: PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,ALEN(type),"boomeramg",&indx,&flg);
2064: if (flg) {
2065: PCHYPRESetType_HYPRE(pc,type[indx]);
2066: } else {
2067: PCHYPRESetType_HYPRE(pc,"boomeramg");
2068: }
2069: if (pc->ops->setfromoptions) {
2070: pc->ops->setfromoptions(PetscOptionsObject,pc);
2071: }
2072: PetscOptionsTail();
2073: return 0;
2074: }
2076: /*@C
2077: PCHYPRESetType - Sets which hypre preconditioner you wish to use
2079: Input Parameters:
2080: + pc - the preconditioner context
2081: - name - either euclid, pilut, parasails, boomeramg, ams, ads
2083: Options Database Keys:
2084: -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
2086: Level: intermediate
2088: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
2089: PCHYPRE
2091: @*/
2092: PetscErrorCode PCHYPRESetType(PC pc,const char name[])
2093: {
2096: PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
2097: return 0;
2098: }
2100: /*@C
2101: PCHYPREGetType - Gets which hypre preconditioner you are using
2103: Input Parameter:
2104: . pc - the preconditioner context
2106: Output Parameter:
2107: . name - either euclid, pilut, parasails, boomeramg, ams, ads
2109: Level: intermediate
2111: .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
2112: PCHYPRE
2114: @*/
2115: PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
2116: {
2119: PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
2120: return 0;
2121: }
2123: /*@C
2124: PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use
2126: Logically Collective on PC
2128: Input Parameters:
2129: + pc - the hypre context
2130: - type - one of 'cusparse', 'hypre'
2132: Options Database Key:
2133: . -pc_mg_galerkin_mat_product_algorithm <cusparse,hypre> - Type of SpGEMM to use in hypre
2135: Level: intermediate
2137: .seealso: PCMGGalerkinGetMatProductAlgorithm()
2139: @*/
2140: PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc,const char name[])
2141: {
2143: PetscTryMethod(pc,"PCMGGalerkinSetMatProductAlgorithm_C",(PC,const char[]),(pc,name));
2144: return 0;
2145: }
2147: /*@C
2148: PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre
2150: Not Collective
2152: Input Parameter:
2153: . pc - the multigrid context
2155: Output Parameter:
2156: . name - one of 'cusparse', 'hypre'
2158: Level: intermediate
2160: .seealso: PCMGGalerkinSetMatProductAlgorithm()
2162: @*/
2163: PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc,const char *name[])
2164: {
2166: PetscTryMethod(pc,"PCMGGalerkinGetMatProductAlgorithm_C",(PC,const char*[]),(pc,name));
2167: return 0;
2168: }
2170: /*MC
2171: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
2173: Options Database Keys:
2174: + -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
2175: . -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
2176: . -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())
2177: - Many others, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX preconditioner
2179: Level: intermediate
2181: Notes:
2182: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
2183: the many hypre options can ONLY be set via the options database (e.g. the command line
2184: or with PetscOptionsSetValue(), there are no functions to set them)
2186: The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations
2187: (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
2188: -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
2189: (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of
2190: iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
2191: and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
2192: then AT MOST twenty V-cycles of boomeramg will be called.
2194: Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
2195: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
2196: Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
2197: If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
2198: and use -ksp_max_it to control the number of V-cycles.
2199: (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
2201: 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
2202: -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
2204: MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2205: the following two options:
2207: See PCPFMG for access to the hypre Struct PFMG solver
2209: GPU Notes:
2210: To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda
2211: Then pass VECCUDA vectors and MATAIJCUSPARSE matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2213: To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip
2214: Then pass VECHIP vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2216: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
2217: PCHYPRESetType(), PCPFMG
2219: M*/
2221: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2222: {
2223: PC_HYPRE *jac;
2225: PetscNewLog(pc,&jac);
2227: pc->data = jac;
2228: pc->ops->reset = PCReset_HYPRE;
2229: pc->ops->destroy = PCDestroy_HYPRE;
2230: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2231: pc->ops->setup = PCSetUp_HYPRE;
2232: pc->ops->apply = PCApply_HYPRE;
2233: jac->comm_hypre = MPI_COMM_NULL;
2234: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
2235: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
2236: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
2237: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
2238: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
2239: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);
2240: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
2241: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
2242: PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinSetMatProductAlgorithm_C",PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG);
2243: PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinGetMatProductAlgorithm_C",PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG);
2244: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2245: #if defined(HYPRE_USING_HIP)
2246: PetscDeviceInitialize(PETSC_DEVICE_HIP);
2247: #endif
2248: #if defined(HYPRE_USING_CUDA)
2249: PetscDeviceInitialize(PETSC_DEVICE_CUDA);
2250: #endif
2251: #endif
2252: return 0;
2253: }
2255: /* ---------------------------------------------------------------------------------------------------------------------------------*/
2257: typedef struct {
2258: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2259: HYPRE_StructSolver hsolver;
2261: /* keep copy of PFMG options used so may view them */
2262: PetscInt its;
2263: double tol;
2264: PetscInt relax_type;
2265: PetscInt rap_type;
2266: PetscInt num_pre_relax,num_post_relax;
2267: PetscInt max_levels;
2268: } PC_PFMG;
2270: PetscErrorCode PCDestroy_PFMG(PC pc)
2271: {
2272: PC_PFMG *ex = (PC_PFMG*) pc->data;
2274: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,ex->hsolver);
2275: PetscCommRestoreComm(PetscObjectComm((PetscObject)pc),&ex->hcomm);
2276: PetscFree(pc->data);
2277: return 0;
2278: }
2280: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
2281: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
2283: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2284: {
2285: PetscBool iascii;
2286: PC_PFMG *ex = (PC_PFMG*) pc->data;
2288: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2289: if (iascii) {
2290: PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");
2291: PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);
2292: PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);
2293: PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);
2294: PetscViewerASCIIPrintf(viewer," RAP type %s\n",PFMGRAPType[ex->rap_type]);
2295: PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2296: PetscViewerASCIIPrintf(viewer," max levels %d\n",ex->max_levels);
2297: }
2298: return 0;
2299: }
2301: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2302: {
2303: PC_PFMG *ex = (PC_PFMG*) pc->data;
2304: PetscBool flg = PETSC_FALSE;
2306: PetscOptionsHead(PetscOptionsObject,"PFMG options");
2307: PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
2308: if (flg) {
2309: PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,ex->hsolver,3);
2310: }
2311: PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
2312: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,ex->hsolver,ex->its);
2313: PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2314: PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,ex->hsolver,ex->num_pre_relax);
2315: PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2316: PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,ex->hsolver,ex->num_post_relax);
2318: PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
2319: PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,ex->hsolver,ex->max_levels);
2321: PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
2322: PetscStackCallStandard(HYPRE_StructPFMGSetTol,ex->hsolver,ex->tol);
2323: PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,ALEN(PFMGRelaxType),PFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
2324: PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,ex->hsolver, ex->relax_type);
2325: PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
2326: PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,ex->hsolver, ex->rap_type);
2327: PetscOptionsTail();
2328: return 0;
2329: }
2331: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2332: {
2333: PC_PFMG *ex = (PC_PFMG*) pc->data;
2334: PetscScalar *yy;
2335: const PetscScalar *xx;
2336: PetscInt ilower[3],iupper[3];
2337: HYPRE_Int hlower[3],hupper[3];
2338: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2340: PetscCitationsRegister(hypreCitation,&cite);
2341: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2342: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2343: iupper[0] += ilower[0] - 1;
2344: iupper[1] += ilower[1] - 1;
2345: iupper[2] += ilower[2] - 1;
2346: hlower[0] = (HYPRE_Int)ilower[0];
2347: hlower[1] = (HYPRE_Int)ilower[1];
2348: hlower[2] = (HYPRE_Int)ilower[2];
2349: hupper[0] = (HYPRE_Int)iupper[0];
2350: hupper[1] = (HYPRE_Int)iupper[1];
2351: hupper[2] = (HYPRE_Int)iupper[2];
2353: /* copy x values over to hypre */
2354: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,mx->hb,0.0);
2355: VecGetArrayRead(x,&xx);
2356: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,mx->hb,hlower,hupper,(HYPRE_Complex*)xx);
2357: VecRestoreArrayRead(x,&xx);
2358: PetscStackCallStandard(HYPRE_StructVectorAssemble,mx->hb);
2359: PetscStackCallStandard(HYPRE_StructPFMGSolve,ex->hsolver,mx->hmat,mx->hb,mx->hx);
2361: /* copy solution values back to PETSc */
2362: VecGetArray(y,&yy);
2363: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,mx->hx,hlower,hupper,(HYPRE_Complex*)yy);
2364: VecRestoreArray(y,&yy);
2365: return 0;
2366: }
2368: static PetscErrorCode PCApplyRichardson_PFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
2369: {
2370: PC_PFMG *jac = (PC_PFMG*)pc->data;
2371: HYPRE_Int oits;
2373: PetscCitationsRegister(hypreCitation,&cite);
2374: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,jac->hsolver,its*jac->its);
2375: PetscStackCallStandard(HYPRE_StructPFMGSetTol,jac->hsolver,rtol);
2377: PCApply_PFMG(pc,b,y);
2378: PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,jac->hsolver,&oits);
2379: *outits = oits;
2380: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2381: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2382: PetscStackCallStandard(HYPRE_StructPFMGSetTol,jac->hsolver,jac->tol);
2383: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,jac->hsolver,jac->its);
2384: return 0;
2385: }
2387: PetscErrorCode PCSetUp_PFMG(PC pc)
2388: {
2389: PC_PFMG *ex = (PC_PFMG*) pc->data;
2390: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2391: PetscBool flg;
2393: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
2396: /* create the hypre solver object and set its information */
2397: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,ex->hsolver);
2398: PetscStackCallStandard(HYPRE_StructPFMGCreate,ex->hcomm,&ex->hsolver);
2399: PetscStackCallStandard(HYPRE_StructPFMGSetup,ex->hsolver,mx->hmat,mx->hb,mx->hx);
2400: PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,ex->hsolver);
2401: return 0;
2402: }
2404: /*MC
2405: PCPFMG - the hypre PFMG multigrid solver
2407: Level: advanced
2409: Options Database:
2410: + -pc_pfmg_its <its> - number of iterations of PFMG to use as preconditioner
2411: . -pc_pfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid solve
2412: . -pc_pfmg_num_post_relax <steps> - number of smoothing steps after coarse grid solve
2413: . -pc_pfmg_tol <tol> - tolerance of PFMG
2414: . -pc_pfmg_relax_type - relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel
2415: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
2417: Notes:
2418: This is for CELL-centered descretizations
2420: This must be used with the MATHYPRESTRUCT matrix type.
2421: This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
2423: .seealso: PCMG, MATHYPRESTRUCT
2424: M*/
2426: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2427: {
2428: PC_PFMG *ex;
2430: PetscNew(&ex); \
2431: pc->data = ex;
2433: ex->its = 1;
2434: ex->tol = 1.e-8;
2435: ex->relax_type = 1;
2436: ex->rap_type = 0;
2437: ex->num_pre_relax = 1;
2438: ex->num_post_relax = 1;
2439: ex->max_levels = 0;
2441: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
2442: pc->ops->view = PCView_PFMG;
2443: pc->ops->destroy = PCDestroy_PFMG;
2444: pc->ops->apply = PCApply_PFMG;
2445: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2446: pc->ops->setup = PCSetUp_PFMG;
2448: PetscCommGetComm(PetscObjectComm((PetscObject)pc),&ex->hcomm);
2449: PetscStackCallStandard(HYPRE_StructPFMGCreate,ex->hcomm,&ex->hsolver);
2450: return 0;
2451: }
2453: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
2455: /* we know we are working with a HYPRE_SStructMatrix */
2456: typedef struct {
2457: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2458: HYPRE_SStructSolver ss_solver;
2460: /* keep copy of SYSPFMG options used so may view them */
2461: PetscInt its;
2462: double tol;
2463: PetscInt relax_type;
2464: PetscInt num_pre_relax,num_post_relax;
2465: } PC_SysPFMG;
2467: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2468: {
2469: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2471: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,ex->ss_solver);
2472: PetscCommRestoreComm(PetscObjectComm((PetscObject)pc),&ex->hcomm);
2473: PetscFree(pc->data);
2474: return 0;
2475: }
2477: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
2479: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2480: {
2481: PetscBool iascii;
2482: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2484: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2485: if (iascii) {
2486: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");
2487: PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);
2488: PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);
2489: PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);
2490: PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2491: }
2492: return 0;
2493: }
2495: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2496: {
2497: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2498: PetscBool flg = PETSC_FALSE;
2500: PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2501: PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2502: if (flg) {
2503: PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,ex->ss_solver,3);
2504: }
2505: PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2506: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,ex->ss_solver,ex->its);
2507: PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2508: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,ex->ss_solver,ex->num_pre_relax);
2509: PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2510: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,ex->ss_solver,ex->num_post_relax);
2512: PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2513: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,ex->ss_solver,ex->tol);
2514: PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,ALEN(SysPFMGRelaxType),SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
2515: PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,ex->ss_solver, ex->relax_type);
2516: PetscOptionsTail();
2517: return 0;
2518: }
2520: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2521: {
2522: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2523: PetscScalar *yy;
2524: const PetscScalar *xx;
2525: PetscInt ilower[3],iupper[3];
2526: HYPRE_Int hlower[3],hupper[3];
2527: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2528: PetscInt ordering= mx->dofs_order;
2529: PetscInt nvars = mx->nvars;
2530: PetscInt part = 0;
2531: PetscInt size;
2532: PetscInt i;
2534: PetscCitationsRegister(hypreCitation,&cite);
2535: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2536: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2537: iupper[0] += ilower[0] - 1;
2538: iupper[1] += ilower[1] - 1;
2539: iupper[2] += ilower[2] - 1;
2540: hlower[0] = (HYPRE_Int)ilower[0];
2541: hlower[1] = (HYPRE_Int)ilower[1];
2542: hlower[2] = (HYPRE_Int)ilower[2];
2543: hupper[0] = (HYPRE_Int)iupper[0];
2544: hupper[1] = (HYPRE_Int)iupper[1];
2545: hupper[2] = (HYPRE_Int)iupper[2];
2547: size = 1;
2548: for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
2550: /* copy x values over to hypre for variable ordering */
2551: if (ordering) {
2552: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0);
2553: VecGetArrayRead(x,&xx);
2554: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i)));
2555: VecRestoreArrayRead(x,&xx);
2556: PetscStackCallStandard(HYPRE_SStructVectorAssemble,mx->ss_b);
2557: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x);
2558: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x);
2560: /* copy solution values back to PETSc */
2561: VecGetArray(y,&yy);
2562: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i)));
2563: VecRestoreArray(y,&yy);
2564: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2565: PetscScalar *z;
2566: PetscInt j, k;
2568: PetscMalloc1(nvars*size,&z);
2569: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0);
2570: VecGetArrayRead(x,&xx);
2572: /* transform nodal to hypre's variable ordering for sys_pfmg */
2573: for (i= 0; i< size; i++) {
2574: k= i*nvars;
2575: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2576: }
2577: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)));
2578: VecRestoreArrayRead(x,&xx);
2579: PetscStackCallStandard(HYPRE_SStructVectorAssemble,mx->ss_b);
2580: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x);
2582: /* copy solution values back to PETSc */
2583: VecGetArray(y,&yy);
2584: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)));
2585: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2586: for (i= 0; i< size; i++) {
2587: k= i*nvars;
2588: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2589: }
2590: VecRestoreArray(y,&yy);
2591: PetscFree(z);
2592: }
2593: return 0;
2594: }
2596: static PetscErrorCode PCApplyRichardson_SysPFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
2597: {
2598: PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
2599: HYPRE_Int oits;
2601: PetscCitationsRegister(hypreCitation,&cite);
2602: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,jac->ss_solver,its*jac->its);
2603: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,jac->ss_solver,rtol);
2604: PCApply_SysPFMG(pc,b,y);
2605: PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,jac->ss_solver,&oits);
2606: *outits = oits;
2607: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2608: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2609: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,jac->ss_solver,jac->tol);
2610: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,jac->ss_solver,jac->its);
2611: return 0;
2612: }
2614: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2615: {
2616: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2617: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2618: PetscBool flg;
2620: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
2623: /* create the hypre sstruct solver object and set its information */
2624: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,ex->ss_solver);
2625: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,ex->hcomm,&ex->ss_solver);
2626: PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,ex->ss_solver);
2627: PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x);
2628: return 0;
2629: }
2631: /*MC
2632: PCSysPFMG - the hypre SysPFMG multigrid solver
2634: Level: advanced
2636: Options Database:
2637: + -pc_syspfmg_its <its> - number of iterations of SysPFMG to use as preconditioner
2638: . -pc_syspfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid
2639: . -pc_syspfmg_num_post_relax <steps> - number of smoothing steps after coarse grid
2640: . -pc_syspfmg_tol <tol> - tolerance of SysPFMG
2641: - -pc_syspfmg_relax_type <Weighted-Jacobi,Red/Black-Gauss-Seidel> - relaxation type for the up and down cycles
2643: Notes:
2644: This is for CELL-centered descretizations
2646: This must be used with the MATHYPRESSTRUCT matrix type.
2647: This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
2648: Also, only cell-centered variables.
2650: .seealso: PCMG, MATHYPRESSTRUCT
2651: M*/
2653: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2654: {
2655: PC_SysPFMG *ex;
2657: PetscNew(&ex); \
2658: pc->data = ex;
2660: ex->its = 1;
2661: ex->tol = 1.e-8;
2662: ex->relax_type = 1;
2663: ex->num_pre_relax = 1;
2664: ex->num_post_relax = 1;
2666: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
2667: pc->ops->view = PCView_SysPFMG;
2668: pc->ops->destroy = PCDestroy_SysPFMG;
2669: pc->ops->apply = PCApply_SysPFMG;
2670: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2671: pc->ops->setup = PCSetUp_SysPFMG;
2673: PetscCommGetComm(PetscObjectComm((PetscObject)pc),&ex->hcomm);
2674: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,ex->hcomm,&ex->ss_solver);
2675: return 0;
2676: }