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