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