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