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