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) PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_IJVectorGetObject,jac->hmnull[nvec]->ij,(void**)&jac->phmnull[nvec]); 284 nvec++; 285 } 286 PetscCallExternal(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 PetscCallExternal(HYPRE_AMSSetDimension,jac->hsolver,jac->dim); 300 } 301 if (jac->constants[0]) { 302 HYPRE_ParVector ozz,zoz,zzo = NULL; 303 PetscCallExternal(HYPRE_IJVectorGetObject,jac->constants[0]->ij,(void**)(&ozz)); 304 PetscCallExternal(HYPRE_IJVectorGetObject,jac->constants[1]->ij,(void**)(&zoz)); 305 if (jac->constants[2]) { 306 PetscCallExternal(HYPRE_IJVectorGetObject,jac->constants[2]->ij,(void**)(&zzo)); 307 } 308 PetscCallExternal(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]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[0]->ij,(void**)(&coords[0])); 316 if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[1]->ij,(void**)(&coords[1])); 317 if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[2]->ij,(void**)(&coords[2])); 318 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr)); 323 PetscCallExternal(HYPRE_AMSSetDiscreteGradient,jac->hsolver,parcsr); 324 if (jac->alpha_Poisson) { 325 hm = (Mat_HYPRE*)(jac->alpha_Poisson->data); 326 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr)); 327 PetscCallExternal(HYPRE_AMSSetAlphaPoissonMatrix,jac->hsolver,parcsr); 328 } 329 if (jac->ams_beta_is_zero) { 330 PetscCallExternal(HYPRE_AMSSetBetaPoissonMatrix,jac->hsolver,NULL); 331 } else if (jac->beta_Poisson) { 332 hm = (Mat_HYPRE*)(jac->beta_Poisson->data); 333 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr)); 334 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsr[i])); 349 } else { 350 nd_parcsr[i] = NULL; 351 } 352 } 353 PetscCallExternal(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]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[0]->ij,(void**)(&coords[0])); 372 if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[1]->ij,(void**)(&coords[1])); 373 if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject,jac->coords[2]->ij,(void**)(&coords[2])); 374 PetscCallExternal(HYPRE_ADSSetCoordinateVectors,jac->hsolver,coords[0],coords[1],coords[2]); 375 } 376 hm = (Mat_HYPRE*)(jac->G->data); 377 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr)); 378 PetscCallExternal(HYPRE_ADSSetDiscreteGradient,jac->hsolver,parcsr); 379 hm = (Mat_HYPRE*)(jac->C->data); 380 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&parcsr)); 381 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hm->ij,(void**)(&nd_parcsr[i])); 410 } else { 411 nd_parcsr[i] = NULL; 412 } 413 } 414 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat); 418 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&bv); 419 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&xv); 420 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat); 438 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&jbv); 439 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&jxv); 440 PetscStackCallExternalVoid("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 PetscCallExternal(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) PetscCallExternal(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,"PCHYPRESetDiscreteGradient_C",NULL)); 504 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL)); 505 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL)); 506 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL)); 507 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL)); 508 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_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 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",NULL)); 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) PetscCallExternal(HYPRE_ParCSRPilutSetMaxIter,jac->hsolver,jac->maxiter); 527 PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag)); 528 if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetDropTolerance,jac->hsolver,jac->tol); 529 PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag)); 530 if (flag) PetscCallExternal(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) PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_IJMatrixGetObject,hjac->ij,(void**)&hmat); 628 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->b->ij,(void**)&jbv); 629 PetscCallExternal(HYPRE_IJVectorGetObject,hjac->x->ij,(void**)&jxv); 630 631 PetscStackCallExternalVoid("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 PetscCallExternal(HYPRE_SetSpGemmUseCusparse,1); 662 PetscFunctionReturn(0); 663 } 664 PetscCall(PetscStrcmp("hypre",jac->spgemm_type,&flag)); 665 if (flag) { 666 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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) PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_BoomerAMGSetSmoothType,jac->hsolver,indx+6); 827 jac->smoothnumlevels = 25; 828 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_BoomerAMGSetRelaxType,jac->hsolver, indx); 865 /* by default, coarse type set to 9 */ 866 jac->relaxtype[2] = 9; 867 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_BoomerAMGSetSmoothType,jac->hsolver,6); 1008 PetscCallExternal(HYPRE_BoomerAMGSetDomainType,jac->hsolver,1); 1009 PetscCallExternal(HYPRE_BoomerAMGSetOverlap,jac->hsolver,0); 1010 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,its*jac->maxiter); 1035 PetscCallExternal(HYPRE_BoomerAMGSetTol,jac->hsolver,rtol); 1036 jac->applyrichardson = PETSC_TRUE; 1037 PetscCall(PCApply_HYPRE(pc,b,y)); 1038 jac->applyrichardson = PETSC_FALSE; 1039 PetscCallExternal(HYPRE_BoomerAMGGetNumIterations,jac->hsolver,&oits); 1040 *outits = oits; 1041 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1042 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1043 PetscCallExternal(HYPRE_BoomerAMGSetTol,jac->hsolver,jac->tol); 1044 PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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 PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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) PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(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_USE_64BIT_INDICES) 1813 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Hypre Euclid does not support 64 bit indices"); 1814 #endif 1815 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc),&jac->comm_hypre)); 1816 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_ParaSailsSetParams,jac->hsolver,jac->threshold,jac->nlevels); 1846 PetscCallExternal(HYPRE_ParaSailsSetFilter,jac->hsolver,jac->filter); 1847 PetscCallExternal(HYPRE_ParaSailsSetLoadbal,jac->hsolver,jac->loadbal); 1848 PetscCallExternal(HYPRE_ParaSailsSetLogging,jac->hsolver,jac->logging); 1849 PetscCallExternal(HYPRE_ParaSailsSetReuse,jac->hsolver,jac->ruse); 1850 PetscCallExternal(HYPRE_ParaSailsSetSym,jac->hsolver,jac->symt); 1851 PetscFunctionReturn(0); 1852 } 1853 PetscCall(PetscStrcmp("boomeramg",jac->hypre_type,&flag)); 1854 if (flag) { 1855 PetscCallExternal(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 PetscCallExternal(HYPRE_BoomerAMGSetCycleType,jac->hsolver,jac->cycletype); 1927 PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels,jac->hsolver,jac->maxlevels); 1928 PetscCallExternal(HYPRE_BoomerAMGSetMaxIter,jac->hsolver,jac->maxiter); 1929 PetscCallExternal(HYPRE_BoomerAMGSetTol,jac->hsolver,jac->tol); 1930 PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor,jac->hsolver,jac->truncfactor); 1931 PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold,jac->hsolver,jac->strongthreshold); 1932 PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum,jac->hsolver,jac->maxrowsum); 1933 PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType,jac->hsolver,jac->coarsentype); 1934 PetscCallExternal(HYPRE_BoomerAMGSetMeasureType,jac->hsolver,jac->measuretype); 1935 PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder,jac->hsolver, jac->relaxorder); 1936 PetscCallExternal(HYPRE_BoomerAMGSetInterpType,jac->hsolver,jac->interptype); 1937 PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels,jac->hsolver,jac->agg_nl); 1938 PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType,jac->hsolver,jac->agg_interptype); 1939 PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts,jac->hsolver,jac->pmax); 1940 PetscCallExternal(HYPRE_BoomerAMGSetNumPaths,jac->hsolver,jac->agg_num_paths); 1941 PetscCallExternal(HYPRE_BoomerAMGSetRelaxType,jac->hsolver, jac->relaxtype[0]); /* defaults coarse to 9 */ 1942 PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps,jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */ 1943 PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize,jac->hsolver, jac->maxc); 1944 PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize,jac->hsolver, jac->minc); 1945 /* GPU */ 1946 #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0) 1947 PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose,jac->hsolver,jac->keeptranspose ? 1 : 0); 1948 PetscCallExternal(HYPRE_BoomerAMGSetRAP2,jac->hsolver, jac->rap2); 1949 PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2,jac->hsolver, jac->mod_rap2); 1950 #endif 1951 1952 /* AIR */ 1953 #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0) 1954 PetscCallExternal(HYPRE_BoomerAMGSetRestriction,jac->hsolver,jac->Rtype); 1955 PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR,jac->hsolver,jac->Rstrongthreshold); 1956 PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR,jac->hsolver,jac->Rfilterthreshold); 1957 PetscCallExternal(HYPRE_BoomerAMGSetADropTol,jac->hsolver,jac->Adroptol); 1958 PetscCallExternal(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 PetscCallExternal(HYPRE_AMSSetPrintLevel,jac->hsolver,jac->as_print); 1998 PetscCallExternal(HYPRE_AMSSetMaxIter,jac->hsolver,jac->as_max_iter); 1999 PetscCallExternal(HYPRE_AMSSetCycleType,jac->hsolver,jac->ams_cycle_type); 2000 PetscCallExternal(HYPRE_AMSSetTol,jac->hsolver,jac->as_tol); 2001 PetscCallExternal(HYPRE_AMSSetSmoothingOptions,jac->hsolver,jac->as_relax_type, 2002 jac->as_relax_times, 2003 jac->as_relax_weight, 2004 jac->as_omega); 2005 PetscCallExternal(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 PetscCallExternal(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 PetscCallExternal(HYPRE_ADSSetPrintLevel,jac->hsolver,jac->as_print); 2059 PetscCallExternal(HYPRE_ADSSetMaxIter,jac->hsolver,jac->as_max_iter); 2060 PetscCallExternal(HYPRE_ADSSetCycleType,jac->hsolver,jac->ams_cycle_type); 2061 PetscCallExternal(HYPRE_ADSSetTol,jac->hsolver,jac->as_tol); 2062 PetscCallExternal(HYPRE_ADSSetSmoothingOptions,jac->hsolver,jac->as_relax_type, 2063 jac->as_relax_times, 2064 jac->as_relax_weight, 2065 jac->as_omega); 2066 PetscCallExternal(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 PetscCallExternal(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) PetscCall(pc->ops->setfromoptions(PetscOptionsObject,pc)); 2106 PetscOptionsHeadEnd(); 2107 PetscFunctionReturn(0); 2108 } 2109 2110 /*@C 2111 PCHYPRESetType - Sets which hypre preconditioner you wish to use 2112 2113 Input Parameters: 2114 + pc - the preconditioner context 2115 - name - either euclid, pilut, parasails, boomeramg, ams, ads 2116 2117 Options Database Keys: 2118 -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads 2119 2120 Level: intermediate 2121 2122 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 2123 `PCHYPRE` 2124 2125 @*/ 2126 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 2127 { 2128 PetscFunctionBegin; 2129 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2130 PetscValidCharPointer(name,2); 2131 PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name)); 2132 PetscFunctionReturn(0); 2133 } 2134 2135 /*@C 2136 PCHYPREGetType - Gets which hypre preconditioner you are using 2137 2138 Input Parameter: 2139 . pc - the preconditioner context 2140 2141 Output Parameter: 2142 . name - either euclid, pilut, parasails, boomeramg, ams, ads 2143 2144 Level: intermediate 2145 2146 .seealso: `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, 2147 `PCHYPRE` 2148 2149 @*/ 2150 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 2151 { 2152 PetscFunctionBegin; 2153 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2154 PetscValidPointer(name,2); 2155 PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name)); 2156 PetscFunctionReturn(0); 2157 } 2158 2159 /*@C 2160 PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use 2161 2162 Logically Collective on PC 2163 2164 Input Parameters: 2165 + pc - the hypre context 2166 - type - one of 'cusparse', 'hypre' 2167 2168 Options Database Key: 2169 . -pc_mg_galerkin_mat_product_algorithm <cusparse,hypre> - Type of SpGEMM to use in hypre 2170 2171 Level: intermediate 2172 2173 .seealso: `PCMGGalerkinGetMatProductAlgorithm()` 2174 2175 @*/ 2176 PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc,const char name[]) 2177 { 2178 PetscFunctionBegin; 2179 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2180 PetscTryMethod(pc,"PCMGGalerkinSetMatProductAlgorithm_C",(PC,const char[]),(pc,name)); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 /*@C 2185 PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre 2186 2187 Not Collective 2188 2189 Input Parameter: 2190 . pc - the multigrid context 2191 2192 Output Parameter: 2193 . name - one of 'cusparse', 'hypre' 2194 2195 Level: intermediate 2196 2197 .seealso: `PCMGGalerkinSetMatProductAlgorithm()` 2198 2199 @*/ 2200 PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc,const char *name[]) 2201 { 2202 PetscFunctionBegin; 2203 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2204 PetscTryMethod(pc,"PCMGGalerkinGetMatProductAlgorithm_C",(PC,const char*[]),(pc,name)); 2205 PetscFunctionReturn(0); 2206 } 2207 2208 /*MC 2209 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 2210 2211 Options Database Keys: 2212 + -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads 2213 . -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal()) 2214 . -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant()) 2215 - Many others, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX preconditioner 2216 2217 Level: intermediate 2218 2219 Notes: 2220 Apart from pc_hypre_type (for which there is PCHYPRESetType()), 2221 the many hypre options can ONLY be set via the options database (e.g. the command line 2222 or with PetscOptionsSetValue(), there are no functions to set them) 2223 2224 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations 2225 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if 2226 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner 2227 (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of 2228 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations 2229 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 2230 then AT MOST twenty V-cycles of boomeramg will be called. 2231 2232 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation 2233 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 2234 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi. 2235 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly 2236 and use -ksp_max_it to control the number of V-cycles. 2237 (see the PETSc FAQ.html at the PETSc website under the Documentation tab). 2238 2239 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option 2240 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L. 2241 2242 MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use 2243 the following two options: 2244 2245 See PCPFMG for access to the hypre Struct PFMG solver 2246 2247 GPU Notes: 2248 To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda 2249 Then pass VECCUDA vectors and MATAIJCUSPARSE matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers. 2250 2251 To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip 2252 Then pass VECHIP vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers. 2253 2254 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, 2255 `PCHYPRESetType()`, `PCPFMG` 2256 2257 M*/ 2258 2259 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 2260 { 2261 PC_HYPRE *jac; 2262 2263 PetscFunctionBegin; 2264 PetscCall(PetscNewLog(pc,&jac)); 2265 2266 pc->data = jac; 2267 pc->ops->reset = PCReset_HYPRE; 2268 pc->ops->destroy = PCDestroy_HYPRE; 2269 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 2270 pc->ops->setup = PCSetUp_HYPRE; 2271 pc->ops->apply = PCApply_HYPRE; 2272 jac->comm_hypre = MPI_COMM_NULL; 2273 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE)); 2274 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE)); 2275 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE)); 2276 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE)); 2277 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE)); 2278 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE)); 2279 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE)); 2280 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE)); 2281 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinSetMatProductAlgorithm_C",PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG)); 2282 PetscCall(PetscObjectComposeFunction((PetscObject)pc,"PCMGGalerkinGetMatProductAlgorithm_C",PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG)); 2283 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2284 #if defined(HYPRE_USING_HIP) 2285 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); 2286 #endif 2287 #if defined(HYPRE_USING_CUDA) 2288 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); 2289 #endif 2290 #endif 2291 PetscFunctionReturn(0); 2292 } 2293 2294 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 2295 2296 typedef struct { 2297 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 2298 HYPRE_StructSolver hsolver; 2299 2300 /* keep copy of PFMG options used so may view them */ 2301 PetscInt its; 2302 double tol; 2303 PetscInt relax_type; 2304 PetscInt rap_type; 2305 PetscInt num_pre_relax,num_post_relax; 2306 PetscInt max_levels; 2307 PetscInt skip_relax; 2308 PetscBool print_statistics; 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) PetscCallExternal(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 PetscCall(PetscViewerASCIIPrintf(viewer," skip relax %" PetscInt_FMT "\n",ex->skip_relax)); 2341 } 2342 PetscFunctionReturn(0); 2343 } 2344 2345 PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2346 { 2347 PC_PFMG *ex = (PC_PFMG*) pc->data; 2348 2349 PetscFunctionBegin; 2350 PetscOptionsHeadBegin(PetscOptionsObject,"PFMG options"); 2351 PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",ex->print_statistics,&ex->print_statistics,NULL)); 2352 PetscCall(PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL)); 2353 PetscCallExternal(HYPRE_StructPFMGSetMaxIter,ex->hsolver,ex->its); 2354 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)); 2355 PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax,ex->hsolver,ex->num_pre_relax); 2356 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)); 2357 PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax,ex->hsolver,ex->num_post_relax); 2358 2359 PetscCall(PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL)); 2360 PetscCallExternal(HYPRE_StructPFMGSetMaxLevels,ex->hsolver,ex->max_levels); 2361 2362 PetscCall(PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL)); 2363 PetscCallExternal(HYPRE_StructPFMGSetTol,ex->hsolver,ex->tol); 2364 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)); 2365 PetscCallExternal(HYPRE_StructPFMGSetRelaxType,ex->hsolver, ex->relax_type); 2366 PetscCall(PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL)); 2367 PetscCallExternal(HYPRE_StructPFMGSetRAPType,ex->hsolver, ex->rap_type); 2368 PetscCall(PetscOptionsInt("-pc_pfmg_skip_relax","Skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations when the underlying problem is isotropic","HYPRE_StructPFMGSetSkipRelax",ex->skip_relax,&ex->skip_relax,NULL)); 2369 PetscCallExternal(HYPRE_StructPFMGSetSkipRelax,ex->hsolver, ex->skip_relax); 2370 PetscOptionsHeadEnd(); 2371 PetscFunctionReturn(0); 2372 } 2373 2374 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 2375 { 2376 PC_PFMG *ex = (PC_PFMG*) pc->data; 2377 PetscScalar *yy; 2378 const PetscScalar *xx; 2379 PetscInt ilower[3],iupper[3]; 2380 HYPRE_Int hlower[3],hupper[3]; 2381 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2382 2383 PetscFunctionBegin; 2384 PetscCall(PetscCitationsRegister(hypreCitation,&cite)); 2385 PetscCall(DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2])); 2386 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2387 iupper[0] += ilower[0] - 1; 2388 iupper[1] += ilower[1] - 1; 2389 iupper[2] += ilower[2] - 1; 2390 hlower[0] = (HYPRE_Int)ilower[0]; 2391 hlower[1] = (HYPRE_Int)ilower[1]; 2392 hlower[2] = (HYPRE_Int)ilower[2]; 2393 hupper[0] = (HYPRE_Int)iupper[0]; 2394 hupper[1] = (HYPRE_Int)iupper[1]; 2395 hupper[2] = (HYPRE_Int)iupper[2]; 2396 2397 /* copy x values over to hypre */ 2398 PetscCallExternal(HYPRE_StructVectorSetConstantValues,mx->hb,0.0); 2399 PetscCall(VecGetArrayRead(x,&xx)); 2400 PetscCallExternal(HYPRE_StructVectorSetBoxValues,mx->hb,hlower,hupper,(HYPRE_Complex*)xx); 2401 PetscCall(VecRestoreArrayRead(x,&xx)); 2402 PetscCallExternal(HYPRE_StructVectorAssemble,mx->hb); 2403 PetscCallExternal(HYPRE_StructPFMGSolve,ex->hsolver,mx->hmat,mx->hb,mx->hx); 2404 2405 /* copy solution values back to PETSc */ 2406 PetscCall(VecGetArray(y,&yy)); 2407 PetscCallExternal(HYPRE_StructVectorGetBoxValues,mx->hx,hlower,hupper,(HYPRE_Complex*)yy); 2408 PetscCall(VecRestoreArray(y,&yy)); 2409 PetscFunctionReturn(0); 2410 } 2411 2412 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) 2413 { 2414 PC_PFMG *jac = (PC_PFMG*)pc->data; 2415 HYPRE_Int oits; 2416 2417 PetscFunctionBegin; 2418 PetscCall(PetscCitationsRegister(hypreCitation,&cite)); 2419 PetscCallExternal(HYPRE_StructPFMGSetMaxIter,jac->hsolver,its*jac->its); 2420 PetscCallExternal(HYPRE_StructPFMGSetTol,jac->hsolver,rtol); 2421 2422 PetscCall(PCApply_PFMG(pc,b,y)); 2423 PetscCallExternal(HYPRE_StructPFMGGetNumIterations,jac->hsolver,&oits); 2424 *outits = oits; 2425 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2426 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2427 PetscCallExternal(HYPRE_StructPFMGSetTol,jac->hsolver,jac->tol); 2428 PetscCallExternal(HYPRE_StructPFMGSetMaxIter,jac->hsolver,jac->its); 2429 PetscFunctionReturn(0); 2430 } 2431 2432 PetscErrorCode PCSetUp_PFMG(PC pc) 2433 { 2434 PC_PFMG *ex = (PC_PFMG*) pc->data; 2435 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2436 PetscBool flg; 2437 2438 PetscFunctionBegin; 2439 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg)); 2440 PetscCheck(flg,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 2441 2442 /* create the hypre solver object and set its information */ 2443 if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy,ex->hsolver); 2444 PetscCallExternal(HYPRE_StructPFMGCreate,ex->hcomm,&ex->hsolver); 2445 2446 // Print Hypre statistics about the solve process 2447 if (ex->print_statistics) PetscCallExternal(HYPRE_StructPFMGSetPrintLevel,ex->hsolver,3); 2448 2449 // The hypre options must be repeated here because the StructPFMG was destroyed and recreated 2450 PetscCallExternal(HYPRE_StructPFMGSetMaxIter ,ex->hsolver, ex->its); 2451 PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax ,ex->hsolver, ex->num_pre_relax); 2452 PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax,ex->hsolver, ex->num_post_relax); 2453 PetscCallExternal(HYPRE_StructPFMGSetMaxLevels ,ex->hsolver, ex->max_levels); 2454 PetscCallExternal(HYPRE_StructPFMGSetTol ,ex->hsolver, ex->tol); 2455 PetscCallExternal(HYPRE_StructPFMGSetRelaxType ,ex->hsolver, ex->relax_type); 2456 PetscCallExternal(HYPRE_StructPFMGSetRAPType ,ex->hsolver, ex->rap_type); 2457 2458 PetscCallExternal(HYPRE_StructPFMGSetup,ex->hsolver,mx->hmat,mx->hb,mx->hx); 2459 PetscCallExternal(HYPRE_StructPFMGSetZeroGuess,ex->hsolver); 2460 PetscFunctionReturn(0); 2461 } 2462 2463 /*MC 2464 PCPFMG - the hypre PFMG multigrid solver 2465 2466 Level: advanced 2467 2468 Options Database: 2469 + -pc_pfmg_its <its> - number of iterations of PFMG to use as preconditioner 2470 . -pc_pfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid solve 2471 . -pc_pfmg_num_post_relax <steps> - number of smoothing steps after coarse grid solve 2472 . -pc_pfmg_tol <tol> - tolerance of PFMG 2473 . -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 2474 . -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 2475 - -pc_pfmg_skip_relax - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations when the underlying problem is isotropic, one of 0,1 2476 2477 Notes: 2478 This is for CELL-centered descretizations 2479 2480 This must be used with the MATHYPRESTRUCT matrix type. 2481 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 2482 2483 .seealso: `PCMG`, `MATHYPRESTRUCT` 2484 M*/ 2485 2486 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 2487 { 2488 PC_PFMG *ex; 2489 2490 PetscFunctionBegin; 2491 PetscCall(PetscNew(&ex)); \ 2492 pc->data = ex; 2493 2494 ex->its = 1; 2495 ex->tol = 1.e-8; 2496 ex->relax_type = 1; 2497 ex->rap_type = 0; 2498 ex->num_pre_relax = 1; 2499 ex->num_post_relax = 1; 2500 ex->max_levels = 0; 2501 ex->skip_relax = 0; 2502 ex->print_statistics = PETSC_FALSE; 2503 2504 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 2505 pc->ops->view = PCView_PFMG; 2506 pc->ops->destroy = PCDestroy_PFMG; 2507 pc->ops->apply = PCApply_PFMG; 2508 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 2509 pc->ops->setup = PCSetUp_PFMG; 2510 2511 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc),&ex->hcomm)); 2512 PetscCallExternal(HYPRE_StructPFMGCreate,ex->hcomm,&ex->hsolver); 2513 PetscFunctionReturn(0); 2514 } 2515 2516 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 2517 2518 /* we know we are working with a HYPRE_SStructMatrix */ 2519 typedef struct { 2520 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 2521 HYPRE_SStructSolver ss_solver; 2522 2523 /* keep copy of SYSPFMG options used so may view them */ 2524 PetscInt its; 2525 double tol; 2526 PetscInt relax_type; 2527 PetscInt num_pre_relax,num_post_relax; 2528 } PC_SysPFMG; 2529 2530 PetscErrorCode PCDestroy_SysPFMG(PC pc) 2531 { 2532 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2533 2534 PetscFunctionBegin; 2535 if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy,ex->ss_solver); 2536 PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc),&ex->hcomm)); 2537 PetscCall(PetscFree(pc->data)); 2538 PetscFunctionReturn(0); 2539 } 2540 2541 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 2542 2543 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 2544 { 2545 PetscBool iascii; 2546 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2547 2548 PetscFunctionBegin; 2549 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 2550 if (iascii) { 2551 PetscCall(PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n")); 2552 PetscCall(PetscViewerASCIIPrintf(viewer," max iterations %" PetscInt_FMT "\n",ex->its)); 2553 PetscCall(PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol)); 2554 PetscCall(PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type])); 2555 PetscCall(PetscViewerASCIIPrintf(viewer," number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n",ex->num_pre_relax,ex->num_post_relax)); 2556 } 2557 PetscFunctionReturn(0); 2558 } 2559 2560 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2561 { 2562 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2563 PetscBool flg = PETSC_FALSE; 2564 2565 PetscFunctionBegin; 2566 PetscOptionsHeadBegin(PetscOptionsObject,"SysPFMG options"); 2567 PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL)); 2568 if (flg) { 2569 PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel,ex->ss_solver,3); 2570 } 2571 PetscCall(PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL)); 2572 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter,ex->ss_solver,ex->its); 2573 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)); 2574 PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax,ex->ss_solver,ex->num_pre_relax); 2575 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)); 2576 PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax,ex->ss_solver,ex->num_post_relax); 2577 2578 PetscCall(PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL)); 2579 PetscCallExternal(HYPRE_SStructSysPFMGSetTol,ex->ss_solver,ex->tol); 2580 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)); 2581 PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType,ex->ss_solver, ex->relax_type); 2582 PetscOptionsHeadEnd(); 2583 PetscFunctionReturn(0); 2584 } 2585 2586 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 2587 { 2588 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2589 PetscScalar *yy; 2590 const PetscScalar *xx; 2591 PetscInt ilower[3],iupper[3]; 2592 HYPRE_Int hlower[3],hupper[3]; 2593 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2594 PetscInt ordering= mx->dofs_order; 2595 PetscInt nvars = mx->nvars; 2596 PetscInt part = 0; 2597 PetscInt size; 2598 PetscInt i; 2599 2600 PetscFunctionBegin; 2601 PetscCall(PetscCitationsRegister(hypreCitation,&cite)); 2602 PetscCall(DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2])); 2603 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2604 iupper[0] += ilower[0] - 1; 2605 iupper[1] += ilower[1] - 1; 2606 iupper[2] += ilower[2] - 1; 2607 hlower[0] = (HYPRE_Int)ilower[0]; 2608 hlower[1] = (HYPRE_Int)ilower[1]; 2609 hlower[2] = (HYPRE_Int)ilower[2]; 2610 hupper[0] = (HYPRE_Int)iupper[0]; 2611 hupper[1] = (HYPRE_Int)iupper[1]; 2612 hupper[2] = (HYPRE_Int)iupper[2]; 2613 2614 size = 1; 2615 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 2616 2617 /* copy x values over to hypre for variable ordering */ 2618 if (ordering) { 2619 PetscCallExternal(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0); 2620 PetscCall(VecGetArrayRead(x,&xx)); 2621 for (i= 0; i< nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))); 2622 PetscCall(VecRestoreArrayRead(x,&xx)); 2623 PetscCallExternal(HYPRE_SStructVectorAssemble,mx->ss_b); 2624 PetscCallExternal(HYPRE_SStructMatrixMatvec,1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x); 2625 PetscCallExternal(HYPRE_SStructSysPFMGSolve,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x); 2626 2627 /* copy solution values back to PETSc */ 2628 PetscCall(VecGetArray(y,&yy)); 2629 for (i= 0; i< nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))); 2630 PetscCall(VecRestoreArray(y,&yy)); 2631 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2632 PetscScalar *z; 2633 PetscInt j, k; 2634 2635 PetscCall(PetscMalloc1(nvars*size,&z)); 2636 PetscCallExternal(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0); 2637 PetscCall(VecGetArrayRead(x,&xx)); 2638 2639 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2640 for (i= 0; i< size; i++) { 2641 k= i*nvars; 2642 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 2643 } 2644 for (i= 0; i< nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))); 2645 PetscCall(VecRestoreArrayRead(x,&xx)); 2646 PetscCallExternal(HYPRE_SStructVectorAssemble,mx->ss_b); 2647 PetscCallExternal(HYPRE_SStructSysPFMGSolve,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x); 2648 2649 /* copy solution values back to PETSc */ 2650 PetscCall(VecGetArray(y,&yy)); 2651 for (i= 0; i< nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))); 2652 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2653 for (i= 0; i< size; i++) { 2654 k= i*nvars; 2655 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 2656 } 2657 PetscCall(VecRestoreArray(y,&yy)); 2658 PetscCall(PetscFree(z)); 2659 } 2660 PetscFunctionReturn(0); 2661 } 2662 2663 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) 2664 { 2665 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 2666 HYPRE_Int oits; 2667 2668 PetscFunctionBegin; 2669 PetscCall(PetscCitationsRegister(hypreCitation,&cite)); 2670 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter,jac->ss_solver,its*jac->its); 2671 PetscCallExternal(HYPRE_SStructSysPFMGSetTol,jac->ss_solver,rtol); 2672 PetscCall(PCApply_SysPFMG(pc,b,y)); 2673 PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations,jac->ss_solver,&oits); 2674 *outits = oits; 2675 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2676 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2677 PetscCallExternal(HYPRE_SStructSysPFMGSetTol,jac->ss_solver,jac->tol); 2678 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter,jac->ss_solver,jac->its); 2679 PetscFunctionReturn(0); 2680 } 2681 2682 PetscErrorCode PCSetUp_SysPFMG(PC pc) 2683 { 2684 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2685 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2686 PetscBool flg; 2687 2688 PetscFunctionBegin; 2689 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg)); 2690 PetscCheck(flg,PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 2691 2692 /* create the hypre sstruct solver object and set its information */ 2693 if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy,ex->ss_solver); 2694 PetscCallExternal(HYPRE_SStructSysPFMGCreate,ex->hcomm,&ex->ss_solver); 2695 PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess,ex->ss_solver); 2696 PetscCallExternal(HYPRE_SStructSysPFMGSetup,ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x); 2697 PetscFunctionReturn(0); 2698 } 2699 2700 /*MC 2701 PCSysPFMG - the hypre SysPFMG multigrid solver 2702 2703 Level: advanced 2704 2705 Options Database: 2706 + -pc_syspfmg_its <its> - number of iterations of SysPFMG to use as preconditioner 2707 . -pc_syspfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid 2708 . -pc_syspfmg_num_post_relax <steps> - number of smoothing steps after coarse grid 2709 . -pc_syspfmg_tol <tol> - tolerance of SysPFMG 2710 - -pc_syspfmg_relax_type <Weighted-Jacobi,Red/Black-Gauss-Seidel> - relaxation type for the up and down cycles 2711 2712 Notes: 2713 This is for CELL-centered descretizations 2714 2715 This must be used with the MATHYPRESSTRUCT matrix type. 2716 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 2717 Also, only cell-centered variables. 2718 2719 .seealso: `PCMG`, `MATHYPRESSTRUCT` 2720 M*/ 2721 2722 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2723 { 2724 PC_SysPFMG *ex; 2725 2726 PetscFunctionBegin; 2727 PetscCall(PetscNew(&ex)); \ 2728 pc->data = ex; 2729 2730 ex->its = 1; 2731 ex->tol = 1.e-8; 2732 ex->relax_type = 1; 2733 ex->num_pre_relax = 1; 2734 ex->num_post_relax = 1; 2735 2736 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2737 pc->ops->view = PCView_SysPFMG; 2738 pc->ops->destroy = PCDestroy_SysPFMG; 2739 pc->ops->apply = PCApply_SysPFMG; 2740 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2741 pc->ops->setup = PCSetUp_SysPFMG; 2742 2743 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc),&ex->hcomm)); 2744 PetscCallExternal(HYPRE_SStructSysPFMGCreate,ex->hcomm,&ex->ss_solver); 2745 PetscFunctionReturn(0); 2746 } 2747