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 PetscReal 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 PetscReal 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 1523 PetscFunctionBegin; 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 PetscFunctionReturn(PETSC_SUCCESS); 1577 } 1578 1579 /*@ 1580 PCHYPRESetInterpolations - Set interpolation matrices for `PCHYPRE` type of ams or ads 1581 1582 Collective 1583 1584 Input Parameters: 1585 + pc - the preconditioning context 1586 . dim - the dimension of the problem, only used in AMS 1587 . RT_PiFull - Raviart-Thomas interpolation matrix 1588 . RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix 1589 . ND_PiFull - Nedelec interpolation matrix 1590 - ND_Pi - x/y/z component of Nedelec interpolation matrix 1591 1592 Level: intermediate 1593 1594 Notes: 1595 For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL. 1596 1597 For ADS, both type of interpolation matrices are needed. 1598 1599 Developer Notes: 1600 This automatically converts the matrix to `MATHYPRE` if it is not already of that type 1601 1602 .seealso: [](ch_ksp), `PCHYPRE` 1603 @*/ 1604 PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) 1605 { 1606 PetscInt i; 1607 1608 PetscFunctionBegin; 1609 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1610 if (RT_PiFull) { 1611 PetscValidHeaderSpecific(RT_PiFull, MAT_CLASSID, 3); 1612 PetscCheckSameComm(pc, 1, RT_PiFull, 3); 1613 } 1614 if (RT_Pi) { 1615 PetscAssertPointer(RT_Pi, 4); 1616 for (i = 0; i < dim; ++i) { 1617 if (RT_Pi[i]) { 1618 PetscValidHeaderSpecific(RT_Pi[i], MAT_CLASSID, 4); 1619 PetscCheckSameComm(pc, 1, RT_Pi[i], 4); 1620 } 1621 } 1622 } 1623 if (ND_PiFull) { 1624 PetscValidHeaderSpecific(ND_PiFull, MAT_CLASSID, 5); 1625 PetscCheckSameComm(pc, 1, ND_PiFull, 5); 1626 } 1627 if (ND_Pi) { 1628 PetscAssertPointer(ND_Pi, 6); 1629 for (i = 0; i < dim; ++i) { 1630 if (ND_Pi[i]) { 1631 PetscValidHeaderSpecific(ND_Pi[i], MAT_CLASSID, 6); 1632 PetscCheckSameComm(pc, 1, ND_Pi[i], 6); 1633 } 1634 } 1635 } 1636 PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi)); 1637 PetscFunctionReturn(PETSC_SUCCESS); 1638 } 1639 1640 static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha) 1641 { 1642 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1643 PetscBool ishypre; 1644 1645 PetscFunctionBegin; 1646 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre)); 1647 if (ishypre) { 1648 if (isalpha) { 1649 PetscCall(PetscObjectReference((PetscObject)A)); 1650 PetscCall(MatDestroy(&jac->alpha_Poisson)); 1651 jac->alpha_Poisson = A; 1652 } else { 1653 if (A) { 1654 PetscCall(PetscObjectReference((PetscObject)A)); 1655 } else { 1656 jac->ams_beta_is_zero = PETSC_TRUE; 1657 } 1658 PetscCall(MatDestroy(&jac->beta_Poisson)); 1659 jac->beta_Poisson = A; 1660 } 1661 } else { 1662 if (isalpha) { 1663 PetscCall(MatDestroy(&jac->alpha_Poisson)); 1664 PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson)); 1665 } else { 1666 if (A) { 1667 PetscCall(MatDestroy(&jac->beta_Poisson)); 1668 PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson)); 1669 } else { 1670 PetscCall(MatDestroy(&jac->beta_Poisson)); 1671 jac->ams_beta_is_zero = PETSC_TRUE; 1672 } 1673 } 1674 } 1675 PetscFunctionReturn(PETSC_SUCCESS); 1676 } 1677 1678 /*@ 1679 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix for `PCHYPRE` of type ams 1680 1681 Collective 1682 1683 Input Parameters: 1684 + pc - the preconditioning context 1685 - A - the matrix 1686 1687 Level: intermediate 1688 1689 Note: 1690 A should be obtained by discretizing the vector valued Poisson problem with linear finite elements 1691 1692 Developer Notes: 1693 This automatically converts the matrix to `MATHYPRE` if it is not already of that type 1694 1695 If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSSetAlphaPoissonMatrix()` 1696 1697 .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()` 1698 @*/ 1699 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) 1700 { 1701 PetscFunctionBegin; 1702 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1703 PetscValidHeaderSpecific(A, MAT_CLASSID, 2); 1704 PetscCheckSameComm(pc, 1, A, 2); 1705 PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE)); 1706 PetscFunctionReturn(PETSC_SUCCESS); 1707 } 1708 1709 /*@ 1710 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix for `PCHYPRE` of type ams 1711 1712 Collective 1713 1714 Input Parameters: 1715 + pc - the preconditioning context 1716 - A - the matrix, or NULL to turn it off 1717 1718 Level: intermediate 1719 1720 Note: 1721 A should be obtained by discretizing the Poisson problem with linear finite elements. 1722 1723 Developer Notes: 1724 This automatically converts the matrix to `MATHYPRE` if it is not already of that type 1725 1726 If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSPCHYPRESetBetaPoissonMatrix()` 1727 1728 .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` 1729 @*/ 1730 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) 1731 { 1732 PetscFunctionBegin; 1733 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1734 if (A) { 1735 PetscValidHeaderSpecific(A, MAT_CLASSID, 2); 1736 PetscCheckSameComm(pc, 1, A, 2); 1737 } 1738 PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE)); 1739 PetscFunctionReturn(PETSC_SUCCESS); 1740 } 1741 1742 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc, Vec ozz, Vec zoz, Vec zzo) 1743 { 1744 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1745 1746 PetscFunctionBegin; 1747 /* throw away any vector if already set */ 1748 PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0])); 1749 PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1])); 1750 PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2])); 1751 PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0])); 1752 PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0])); 1753 PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1])); 1754 PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1])); 1755 jac->dim = 2; 1756 if (zzo) { 1757 PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2])); 1758 PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2])); 1759 jac->dim++; 1760 } 1761 PetscFunctionReturn(PETSC_SUCCESS); 1762 } 1763 1764 /*@ 1765 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in the edge element basis for `PCHYPRE` of type ams 1766 1767 Collective 1768 1769 Input Parameters: 1770 + pc - the preconditioning context 1771 . ozz - vector representing (1,0,0) (or (1,0) in 2D) 1772 . zoz - vector representing (0,1,0) (or (0,1) in 2D) 1773 - zzo - vector representing (0,0,1) (use NULL in 2D) 1774 1775 Level: intermediate 1776 1777 Developer Notes: 1778 If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSSetEdgeConstantVectors()` 1779 1780 .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` 1781 @*/ 1782 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) 1783 { 1784 PetscFunctionBegin; 1785 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1786 PetscValidHeaderSpecific(ozz, VEC_CLASSID, 2); 1787 PetscValidHeaderSpecific(zoz, VEC_CLASSID, 3); 1788 if (zzo) PetscValidHeaderSpecific(zzo, VEC_CLASSID, 4); 1789 PetscCheckSameComm(pc, 1, ozz, 2); 1790 PetscCheckSameComm(pc, 1, zoz, 3); 1791 if (zzo) PetscCheckSameComm(pc, 1, zzo, 4); 1792 PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo)); 1793 PetscFunctionReturn(PETSC_SUCCESS); 1794 } 1795 1796 static PetscErrorCode PCHYPREAMSSetInteriorNodes_HYPRE(PC pc, Vec interior) 1797 { 1798 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1799 1800 PetscFunctionBegin; 1801 PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior)); 1802 PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior)); 1803 PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior)); 1804 jac->ams_beta_is_zero_part = PETSC_TRUE; 1805 PetscFunctionReturn(PETSC_SUCCESS); 1806 } 1807 1808 /*@ 1809 PCHYPREAMSSetInteriorNodes - Set the list of interior nodes to a zero-conductivity region for `PCHYPRE` of type ams 1810 1811 Collective 1812 1813 Input Parameters: 1814 + pc - the preconditioning context 1815 - interior - vector. node is interior if its entry in the array is 1.0. 1816 1817 Level: intermediate 1818 1819 Note: 1820 This calls `HYPRE_AMSSetInteriorNodes()` 1821 1822 Developer Notes: 1823 If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSSetInteriorNodes()` 1824 1825 .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` 1826 @*/ 1827 PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior) 1828 { 1829 PetscFunctionBegin; 1830 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 1831 PetscValidHeaderSpecific(interior, VEC_CLASSID, 2); 1832 PetscCheckSameComm(pc, 1, interior, 2); 1833 PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior)); 1834 PetscFunctionReturn(PETSC_SUCCESS); 1835 } 1836 1837 static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 1838 { 1839 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1840 Vec tv; 1841 PetscInt i; 1842 1843 PetscFunctionBegin; 1844 /* throw away any coordinate vector if already set */ 1845 PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0])); 1846 PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1])); 1847 PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2])); 1848 jac->dim = dim; 1849 1850 /* compute IJ vector for coordinates */ 1851 PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &tv)); 1852 PetscCall(VecSetType(tv, VECSTANDARD)); 1853 PetscCall(VecSetSizes(tv, nloc, PETSC_DECIDE)); 1854 for (i = 0; i < dim; i++) { 1855 PetscScalar *array; 1856 PetscInt j; 1857 1858 PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i])); 1859 PetscCall(VecGetArrayWrite(tv, &array)); 1860 for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i]; 1861 PetscCall(VecRestoreArrayWrite(tv, &array)); 1862 PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i])); 1863 } 1864 PetscCall(VecDestroy(&tv)); 1865 PetscFunctionReturn(PETSC_SUCCESS); 1866 } 1867 1868 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc, const char *name[]) 1869 { 1870 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1871 1872 PetscFunctionBegin; 1873 *name = jac->hypre_type; 1874 PetscFunctionReturn(PETSC_SUCCESS); 1875 } 1876 1877 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc, const char name[]) 1878 { 1879 PC_HYPRE *jac = (PC_HYPRE *)pc->data; 1880 PetscBool flag; 1881 1882 PetscFunctionBegin; 1883 if (jac->hypre_type) { 1884 PetscCall(PetscStrcmp(jac->hypre_type, name, &flag)); 1885 PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE preconditioner type once it has been set"); 1886 PetscFunctionReturn(PETSC_SUCCESS); 1887 } else { 1888 PetscCall(PetscStrallocpy(name, &jac->hypre_type)); 1889 } 1890 1891 jac->maxiter = PETSC_DEFAULT; 1892 jac->tol = PETSC_DEFAULT; 1893 jac->printstatistics = PetscLogPrintInfo; 1894 1895 PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag)); 1896 if (flag) { 1897 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); 1898 PetscCallExternal(HYPRE_ParCSRPilutCreate, jac->comm_hypre, &jac->hsolver); 1899 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; 1900 pc->ops->view = PCView_HYPRE_Pilut; 1901 jac->destroy = HYPRE_ParCSRPilutDestroy; 1902 jac->setup = HYPRE_ParCSRPilutSetup; 1903 jac->solve = HYPRE_ParCSRPilutSolve; 1904 jac->factorrowsize = PETSC_DEFAULT; 1905 PetscFunctionReturn(PETSC_SUCCESS); 1906 } 1907 PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag)); 1908 if (flag) { 1909 #if defined(PETSC_USE_64BIT_INDICES) 1910 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices"); 1911 #endif 1912 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); 1913 PetscCallExternal(HYPRE_EuclidCreate, jac->comm_hypre, &jac->hsolver); 1914 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid; 1915 pc->ops->view = PCView_HYPRE_Euclid; 1916 jac->destroy = HYPRE_EuclidDestroy; 1917 jac->setup = HYPRE_EuclidSetup; 1918 jac->solve = HYPRE_EuclidSolve; 1919 jac->factorrowsize = PETSC_DEFAULT; 1920 jac->eu_level = PETSC_DEFAULT; /* default */ 1921 PetscFunctionReturn(PETSC_SUCCESS); 1922 } 1923 PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag)); 1924 if (flag) { 1925 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); 1926 PetscCallExternal(HYPRE_ParaSailsCreate, jac->comm_hypre, &jac->hsolver); 1927 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; 1928 pc->ops->view = PCView_HYPRE_ParaSails; 1929 jac->destroy = HYPRE_ParaSailsDestroy; 1930 jac->setup = HYPRE_ParaSailsSetup; 1931 jac->solve = HYPRE_ParaSailsSolve; 1932 /* initialize */ 1933 jac->nlevels = 1; 1934 jac->threshold = .1; 1935 jac->filter = .1; 1936 jac->loadbal = 0; 1937 if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE; 1938 else jac->logging = (int)PETSC_FALSE; 1939 1940 jac->ruse = (int)PETSC_FALSE; 1941 jac->symt = 0; 1942 PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels); 1943 PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter); 1944 PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal); 1945 PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging); 1946 PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse); 1947 PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt); 1948 PetscFunctionReturn(PETSC_SUCCESS); 1949 } 1950 PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag)); 1951 if (flag) { 1952 PetscCallExternal(HYPRE_BoomerAMGCreate, &jac->hsolver); 1953 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; 1954 pc->ops->view = PCView_HYPRE_BoomerAMG; 1955 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; 1956 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; 1957 pc->ops->matapply = PCMatApply_HYPRE_BoomerAMG; 1958 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG)); 1959 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG)); 1960 jac->destroy = HYPRE_BoomerAMGDestroy; 1961 jac->setup = HYPRE_BoomerAMGSetup; 1962 jac->solve = HYPRE_BoomerAMGSolve; 1963 jac->applyrichardson = PETSC_FALSE; 1964 /* these defaults match the hypre defaults */ 1965 jac->cycletype = 1; 1966 jac->maxlevels = 25; 1967 jac->maxiter = 1; 1968 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ 1969 jac->truncfactor = 0.0; 1970 jac->strongthreshold = .25; 1971 jac->maxrowsum = .9; 1972 jac->coarsentype = 6; 1973 jac->measuretype = 0; 1974 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; 1975 jac->smoothtype = -1; /* Not set by default */ 1976 jac->smoothnumlevels = 25; 1977 jac->eu_level = 0; 1978 jac->eu_droptolerance = 0; 1979 jac->eu_bj = 0; 1980 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */ 1981 jac->relaxtype[2] = 9; /*G.E. */ 1982 jac->relaxweight = 1.0; 1983 jac->outerrelaxweight = 1.0; 1984 jac->relaxorder = 1; 1985 jac->interptype = 0; 1986 jac->Rtype = 0; 1987 jac->Rstrongthreshold = 0.25; 1988 jac->Rfilterthreshold = 0.0; 1989 jac->Adroptype = -1; 1990 jac->Adroptol = 0.0; 1991 jac->agg_nl = 0; 1992 jac->agg_interptype = 4; 1993 jac->pmax = 0; 1994 jac->truncfactor = 0.0; 1995 jac->agg_num_paths = 1; 1996 jac->maxc = 9; 1997 jac->minc = 1; 1998 jac->nodal_coarsening = 0; 1999 jac->nodal_coarsening_diag = 0; 2000 jac->vec_interp_variant = 0; 2001 jac->vec_interp_qmax = 0; 2002 jac->vec_interp_smooth = PETSC_FALSE; 2003 jac->interp_refine = 0; 2004 jac->nodal_relax = PETSC_FALSE; 2005 jac->nodal_relax_levels = 1; 2006 jac->rap2 = 0; 2007 2008 /* GPU defaults 2009 from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options 2010 and /src/parcsr_ls/par_amg.c */ 2011 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2012 jac->keeptranspose = PETSC_TRUE; 2013 jac->mod_rap2 = 1; 2014 jac->coarsentype = 8; 2015 jac->relaxorder = 0; 2016 jac->interptype = 6; 2017 jac->relaxtype[0] = 18; 2018 jac->relaxtype[1] = 18; 2019 jac->agg_interptype = 7; 2020 #else 2021 jac->keeptranspose = PETSC_FALSE; 2022 jac->mod_rap2 = 0; 2023 #endif 2024 PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype); 2025 PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels); 2026 PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter); 2027 PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol); 2028 PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor); 2029 PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold); 2030 PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum); 2031 PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype); 2032 PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype); 2033 PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder); 2034 PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype); 2035 PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl); 2036 PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType, jac->hsolver, jac->agg_interptype); 2037 PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax); 2038 PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths); 2039 PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, jac->relaxtype[0]); /* defaults coarse to 9 */ 2040 PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */ 2041 PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc); 2042 PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc); 2043 /* GPU */ 2044 #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) 2045 PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0); 2046 PetscCallExternal(HYPRE_BoomerAMGSetRAP2, jac->hsolver, jac->rap2); 2047 PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2, jac->hsolver, jac->mod_rap2); 2048 #endif 2049 2050 /* AIR */ 2051 #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) 2052 PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype); 2053 PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold); 2054 PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold); 2055 PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol); 2056 PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype); 2057 #endif 2058 PetscFunctionReturn(PETSC_SUCCESS); 2059 } 2060 PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag)); 2061 if (flag) { 2062 PetscCallExternal(HYPRE_AMSCreate, &jac->hsolver); 2063 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 2064 pc->ops->view = PCView_HYPRE_AMS; 2065 jac->destroy = HYPRE_AMSDestroy; 2066 jac->setup = HYPRE_AMSSetup; 2067 jac->solve = HYPRE_AMSSolve; 2068 jac->coords[0] = NULL; 2069 jac->coords[1] = NULL; 2070 jac->coords[2] = NULL; 2071 jac->interior = NULL; 2072 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 2073 jac->as_print = 0; 2074 jac->as_max_iter = 1; /* used as a preconditioner */ 2075 jac->as_tol = 0.; /* used as a preconditioner */ 2076 jac->ams_cycle_type = 13; 2077 /* Smoothing options */ 2078 jac->as_relax_type = 2; 2079 jac->as_relax_times = 1; 2080 jac->as_relax_weight = 1.0; 2081 jac->as_omega = 1.0; 2082 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2083 jac->as_amg_alpha_opts[0] = 10; 2084 jac->as_amg_alpha_opts[1] = 1; 2085 jac->as_amg_alpha_opts[2] = 6; 2086 jac->as_amg_alpha_opts[3] = 6; 2087 jac->as_amg_alpha_opts[4] = 4; 2088 jac->as_amg_alpha_theta = 0.25; 2089 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2090 jac->as_amg_beta_opts[0] = 10; 2091 jac->as_amg_beta_opts[1] = 1; 2092 jac->as_amg_beta_opts[2] = 6; 2093 jac->as_amg_beta_opts[3] = 6; 2094 jac->as_amg_beta_opts[4] = 4; 2095 jac->as_amg_beta_theta = 0.25; 2096 PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print); 2097 PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter); 2098 PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type); 2099 PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol); 2100 PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega); 2101 PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 2102 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 2103 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 2104 jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */ 2105 jac->as_amg_alpha_opts[4]); /* AMG Pmax */ 2106 PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */ 2107 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 2108 jac->as_amg_beta_opts[2], /* AMG relax_type */ 2109 jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */ 2110 jac->as_amg_beta_opts[4]); /* AMG Pmax */ 2111 /* Zero conductivity */ 2112 jac->ams_beta_is_zero = PETSC_FALSE; 2113 jac->ams_beta_is_zero_part = PETSC_FALSE; 2114 PetscFunctionReturn(PETSC_SUCCESS); 2115 } 2116 PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag)); 2117 if (flag) { 2118 PetscCallExternal(HYPRE_ADSCreate, &jac->hsolver); 2119 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS; 2120 pc->ops->view = PCView_HYPRE_ADS; 2121 jac->destroy = HYPRE_ADSDestroy; 2122 jac->setup = HYPRE_ADSSetup; 2123 jac->solve = HYPRE_ADSSolve; 2124 jac->coords[0] = NULL; 2125 jac->coords[1] = NULL; 2126 jac->coords[2] = NULL; 2127 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */ 2128 jac->as_print = 0; 2129 jac->as_max_iter = 1; /* used as a preconditioner */ 2130 jac->as_tol = 0.; /* used as a preconditioner */ 2131 jac->ads_cycle_type = 13; 2132 /* Smoothing options */ 2133 jac->as_relax_type = 2; 2134 jac->as_relax_times = 1; 2135 jac->as_relax_weight = 1.0; 2136 jac->as_omega = 1.0; 2137 /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2138 jac->ams_cycle_type = 14; 2139 jac->as_amg_alpha_opts[0] = 10; 2140 jac->as_amg_alpha_opts[1] = 1; 2141 jac->as_amg_alpha_opts[2] = 6; 2142 jac->as_amg_alpha_opts[3] = 6; 2143 jac->as_amg_alpha_opts[4] = 4; 2144 jac->as_amg_alpha_theta = 0.25; 2145 /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2146 jac->as_amg_beta_opts[0] = 10; 2147 jac->as_amg_beta_opts[1] = 1; 2148 jac->as_amg_beta_opts[2] = 6; 2149 jac->as_amg_beta_opts[3] = 6; 2150 jac->as_amg_beta_opts[4] = 4; 2151 jac->as_amg_beta_theta = 0.25; 2152 PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print); 2153 PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter); 2154 PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ams_cycle_type); 2155 PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol); 2156 PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega); 2157 PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type, /* AMG coarsen type */ 2158 jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 2159 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 2160 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 2161 jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */ 2162 jac->as_amg_alpha_opts[4]); /* AMG Pmax */ 2163 PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */ 2164 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 2165 jac->as_amg_beta_opts[2], /* AMG relax_type */ 2166 jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */ 2167 jac->as_amg_beta_opts[4]); /* AMG Pmax */ 2168 PetscFunctionReturn(PETSC_SUCCESS); 2169 } 2170 PetscCall(PetscFree(jac->hypre_type)); 2171 2172 jac->hypre_type = NULL; 2173 SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams", name); 2174 } 2175 2176 /* 2177 It only gets here if the HYPRE type has not been set before the call to 2178 ...SetFromOptions() which actually is most of the time 2179 */ 2180 static PetscErrorCode PCSetFromOptions_HYPRE(PC pc, PetscOptionItems *PetscOptionsObject) 2181 { 2182 PetscInt indx; 2183 const char *type[] = {"euclid", "pilut", "parasails", "boomeramg", "ams", "ads"}; 2184 PetscBool flg; 2185 2186 PetscFunctionBegin; 2187 PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options"); 2188 PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg)); 2189 if (flg) { 2190 PetscCall(PCHYPRESetType_HYPRE(pc, type[indx])); 2191 } else { 2192 PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg")); 2193 } 2194 PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject); 2195 PetscOptionsHeadEnd(); 2196 PetscFunctionReturn(PETSC_SUCCESS); 2197 } 2198 2199 /*@C 2200 PCHYPRESetType - Sets which hypre preconditioner you wish to use 2201 2202 Input Parameters: 2203 + pc - the preconditioner context 2204 - name - either euclid, pilut, parasails, boomeramg, ams, ads 2205 2206 Options Database Key: 2207 . pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads 2208 2209 Level: intermediate 2210 2211 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRE` 2212 @*/ 2213 PetscErrorCode PCHYPRESetType(PC pc, const char name[]) 2214 { 2215 PetscFunctionBegin; 2216 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2217 PetscAssertPointer(name, 2); 2218 PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name)); 2219 PetscFunctionReturn(PETSC_SUCCESS); 2220 } 2221 2222 /*@C 2223 PCHYPREGetType - Gets which hypre preconditioner you are using 2224 2225 Input Parameter: 2226 . pc - the preconditioner context 2227 2228 Output Parameter: 2229 . name - either euclid, pilut, parasails, boomeramg, ams, ads 2230 2231 Level: intermediate 2232 2233 .seealso: [](ch_ksp), `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, `PCHYPRE` 2234 @*/ 2235 PetscErrorCode PCHYPREGetType(PC pc, const char *name[]) 2236 { 2237 PetscFunctionBegin; 2238 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2239 PetscAssertPointer(name, 2); 2240 PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name)); 2241 PetscFunctionReturn(PETSC_SUCCESS); 2242 } 2243 2244 /*@C 2245 PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use on GPUs 2246 2247 Logically Collective 2248 2249 Input Parameters: 2250 + pc - the hypre context 2251 - name - one of 'cusparse', 'hypre' 2252 2253 Options Database Key: 2254 . -pc_mg_galerkin_mat_product_algorithm <cusparse,hypre> - Type of SpGEMM to use in hypre 2255 2256 Level: intermediate 2257 2258 Developer Notes: 2259 How the name starts with `PCMG`, should it not be `PCHYPREBoomerAMG`? 2260 2261 .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinGetMatProductAlgorithm()` 2262 @*/ 2263 PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[]) 2264 { 2265 PetscFunctionBegin; 2266 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2267 PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name)); 2268 PetscFunctionReturn(PETSC_SUCCESS); 2269 } 2270 2271 /*@C 2272 PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre to use on GPUs 2273 2274 Not Collective 2275 2276 Input Parameter: 2277 . pc - the multigrid context 2278 2279 Output Parameter: 2280 . name - one of 'cusparse', 'hypre' 2281 2282 Level: intermediate 2283 2284 .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinSetMatProductAlgorithm()` 2285 @*/ 2286 PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[]) 2287 { 2288 PetscFunctionBegin; 2289 PetscValidHeaderSpecific(pc, PC_CLASSID, 1); 2290 PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name)); 2291 PetscFunctionReturn(PETSC_SUCCESS); 2292 } 2293 2294 /*MC 2295 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre as PETSc `PC` 2296 2297 Options Database Keys: 2298 + -pc_hypre_type - One of `euclid`, `pilut`, `parasails`, `boomeramg`, `ams`, or `ads` 2299 . -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see `HYPRE_BOOMERAMGSetNodal()`) 2300 . -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see `HYPRE_BoomerAMGSetInterpVecVariant()`) 2301 - Many others, run with `-pc_type hypre` `-pc_hypre_type XXX` `-help` to see options for the XXX preconditioner 2302 2303 Level: intermediate 2304 2305 Notes: 2306 Apart from `-pc_hypre_type` (for which there is `PCHYPRESetType()`), 2307 the many hypre options can ONLY be set via the options database (e.g. the command line 2308 or with `PetscOptionsSetValue()`, there are no functions to set them) 2309 2310 The options `-pc_hypre_boomeramg_max_iter` and `-pc_hypre_boomeramg_tol` refer to the number of iterations 2311 (V-cycles) and tolerance that boomerAMG does EACH time it is called. So for example, if 2312 `-pc_hypre_boomeramg_max_iter` is set to 2 then 2-V-cycles are being used to define the preconditioner 2313 (`-pc_hypre_boomeramg_tol` should be set to 0.0 - the default - to strictly use a fixed number of 2314 iterations per hypre call). `-ksp_max_it` and `-ksp_rtol` STILL determine the total number of iterations 2315 and tolerance for the Krylov solver. For example, if `-pc_hypre_boomeramg_max_iter` is 2 and `-ksp_max_it` is 10 2316 then AT MOST twenty V-cycles of boomeramg will be used. 2317 2318 Note that the option `-pc_hypre_boomeramg_relax_type_all` defaults to symmetric relaxation 2319 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 2320 Otherwise, you may want to use `-pc_hypre_boomeramg_relax_type_all SOR/Jacobi`. 2321 2322 `MatSetNearNullSpace()` - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use 2323 the following two options: `-pc_hypre_boomeramg_nodal_coarsen <n> -pc_hypre_boomeramg_vec_interp_variant <v>` 2324 2325 See `PCPFMG`, `PCSMG`, and `PCSYSPFMG` for access to hypre's other (nonalgebraic) multigrid solvers 2326 2327 For `PCHYPRE` type of `ams` or `ads` auxiliary data must be provided to the preconditioner with `PCHYPRESetDiscreteGradient()`, 2328 `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`, 2329 `PCHYPREAMSSetInteriorNodes()` 2330 2331 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 2332 since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON` 2333 (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid. 2334 2335 PETSc provides its own geometric and algebraic multigrid solvers `PCMG` and `PCGAMG`, also see `PCHMG` which is useful for certain multicomponent problems 2336 2337 GPU Notes: 2338 To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda 2339 Then pass `VECCUDA` vectors and `MATAIJCUSPARSE` matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers. 2340 2341 To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip 2342 Then pass `VECHIP` vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers. 2343 2344 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`, 2345 `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`, 2346 PCHYPREAMSSetInteriorNodes() 2347 M*/ 2348 2349 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 2350 { 2351 PC_HYPRE *jac; 2352 2353 PetscFunctionBegin; 2354 PetscCall(PetscNew(&jac)); 2355 2356 pc->data = jac; 2357 pc->ops->reset = PCReset_HYPRE; 2358 pc->ops->destroy = PCDestroy_HYPRE; 2359 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 2360 pc->ops->setup = PCSetUp_HYPRE; 2361 pc->ops->apply = PCApply_HYPRE; 2362 jac->comm_hypre = MPI_COMM_NULL; 2363 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE)); 2364 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE)); 2365 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE)); 2366 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE)); 2367 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE)); 2368 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE)); 2369 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE)); 2370 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE)); 2371 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE)); 2372 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG)); 2373 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG)); 2374 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2375 #if defined(HYPRE_USING_HIP) 2376 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); 2377 #endif 2378 #if defined(HYPRE_USING_CUDA) 2379 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); 2380 #endif 2381 #endif 2382 PetscHYPREInitialize(); 2383 PetscFunctionReturn(PETSC_SUCCESS); 2384 } 2385 2386 typedef struct { 2387 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 2388 HYPRE_StructSolver hsolver; 2389 2390 /* keep copy of PFMG options used so may view them */ 2391 PetscInt its; 2392 PetscReal tol; 2393 PetscInt relax_type; 2394 PetscInt rap_type; 2395 PetscInt num_pre_relax, num_post_relax; 2396 PetscInt max_levels; 2397 PetscInt skip_relax; 2398 PetscBool print_statistics; 2399 } PC_PFMG; 2400 2401 static PetscErrorCode PCDestroy_PFMG(PC pc) 2402 { 2403 PC_PFMG *ex = (PC_PFMG *)pc->data; 2404 2405 PetscFunctionBegin; 2406 if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver); 2407 PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 2408 PetscCall(PetscFree(pc->data)); 2409 PetscFunctionReturn(PETSC_SUCCESS); 2410 } 2411 2412 static const char *PFMGRelaxType[] = {"Jacobi", "Weighted-Jacobi", "symmetric-Red/Black-Gauss-Seidel", "Red/Black-Gauss-Seidel"}; 2413 static const char *PFMGRAPType[] = {"Galerkin", "non-Galerkin"}; 2414 2415 static PetscErrorCode PCView_PFMG(PC pc, PetscViewer viewer) 2416 { 2417 PetscBool iascii; 2418 PC_PFMG *ex = (PC_PFMG *)pc->data; 2419 2420 PetscFunctionBegin; 2421 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 2422 if (iascii) { 2423 PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE PFMG preconditioning\n")); 2424 PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); 2425 PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); 2426 PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type])); 2427 PetscCall(PetscViewerASCIIPrintf(viewer, " RAP type %s\n", PFMGRAPType[ex->rap_type])); 2428 PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); 2429 PetscCall(PetscViewerASCIIPrintf(viewer, " max levels %" PetscInt_FMT "\n", ex->max_levels)); 2430 PetscCall(PetscViewerASCIIPrintf(viewer, " skip relax %" PetscInt_FMT "\n", ex->skip_relax)); 2431 } 2432 PetscFunctionReturn(PETSC_SUCCESS); 2433 } 2434 2435 static PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems *PetscOptionsObject) 2436 { 2437 PC_PFMG *ex = (PC_PFMG *)pc->data; 2438 2439 PetscFunctionBegin; 2440 PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options"); 2441 PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL)); 2442 PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL)); 2443 PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its); 2444 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)); 2445 PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax); 2446 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)); 2447 PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax); 2448 2449 PetscCall(PetscOptionsInt("-pc_pfmg_max_levels", "Max Levels for MG hierarchy", "HYPRE_StructPFMGSetMaxLevels", ex->max_levels, &ex->max_levels, NULL)); 2450 PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels); 2451 2452 PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL)); 2453 PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol); 2454 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)); 2455 PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type); 2456 PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL)); 2457 PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type); 2458 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)); 2459 PetscCallExternal(HYPRE_StructPFMGSetSkipRelax, ex->hsolver, ex->skip_relax); 2460 PetscOptionsHeadEnd(); 2461 PetscFunctionReturn(PETSC_SUCCESS); 2462 } 2463 2464 static PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y) 2465 { 2466 PC_PFMG *ex = (PC_PFMG *)pc->data; 2467 PetscScalar *yy; 2468 const PetscScalar *xx; 2469 PetscInt ilower[3], iupper[3]; 2470 HYPRE_Int hlower[3], hupper[3]; 2471 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; 2472 2473 PetscFunctionBegin; 2474 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2475 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 2476 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2477 iupper[0] += ilower[0] - 1; 2478 iupper[1] += ilower[1] - 1; 2479 iupper[2] += ilower[2] - 1; 2480 hlower[0] = (HYPRE_Int)ilower[0]; 2481 hlower[1] = (HYPRE_Int)ilower[1]; 2482 hlower[2] = (HYPRE_Int)ilower[2]; 2483 hupper[0] = (HYPRE_Int)iupper[0]; 2484 hupper[1] = (HYPRE_Int)iupper[1]; 2485 hupper[2] = (HYPRE_Int)iupper[2]; 2486 2487 /* copy x values over to hypre */ 2488 PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0); 2489 PetscCall(VecGetArrayRead(x, &xx)); 2490 PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx); 2491 PetscCall(VecRestoreArrayRead(x, &xx)); 2492 PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb); 2493 PetscCallExternal(HYPRE_StructPFMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx); 2494 2495 /* copy solution values back to PETSc */ 2496 PetscCall(VecGetArray(y, &yy)); 2497 PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy); 2498 PetscCall(VecRestoreArray(y, &yy)); 2499 PetscFunctionReturn(PETSC_SUCCESS); 2500 } 2501 2502 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) 2503 { 2504 PC_PFMG *jac = (PC_PFMG *)pc->data; 2505 HYPRE_Int oits; 2506 2507 PetscFunctionBegin; 2508 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2509 PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, its * jac->its); 2510 PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, rtol); 2511 2512 PetscCall(PCApply_PFMG(pc, b, y)); 2513 PetscCallExternal(HYPRE_StructPFMGGetNumIterations, jac->hsolver, &oits); 2514 *outits = oits; 2515 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2516 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2517 PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, jac->tol); 2518 PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, jac->its); 2519 PetscFunctionReturn(PETSC_SUCCESS); 2520 } 2521 2522 static PetscErrorCode PCSetUp_PFMG(PC pc) 2523 { 2524 PC_PFMG *ex = (PC_PFMG *)pc->data; 2525 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; 2526 PetscBool flg; 2527 2528 PetscFunctionBegin; 2529 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg)); 2530 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner"); 2531 2532 /* create the hypre solver object and set its information */ 2533 if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver); 2534 PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver); 2535 2536 // Print Hypre statistics about the solve process 2537 if (ex->print_statistics) PetscCallExternal(HYPRE_StructPFMGSetPrintLevel, ex->hsolver, 3); 2538 2539 // The hypre options must be repeated here because the StructPFMG was destroyed and recreated 2540 PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its); 2541 PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax); 2542 PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax); 2543 PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels); 2544 PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol); 2545 PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type); 2546 PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type); 2547 2548 PetscCallExternal(HYPRE_StructPFMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx); 2549 PetscCallExternal(HYPRE_StructPFMGSetZeroGuess, ex->hsolver); 2550 PetscFunctionReturn(PETSC_SUCCESS); 2551 } 2552 2553 /*MC 2554 PCPFMG - the hypre PFMG multigrid solver 2555 2556 Options Database Keys: 2557 + -pc_pfmg_its <its> - number of iterations of PFMG to use as preconditioner 2558 . -pc_pfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid solve 2559 . -pc_pfmg_num_post_relax <steps> - number of smoothing steps after coarse grid solve 2560 . -pc_pfmg_tol <tol> - tolerance of PFMG 2561 . -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 2562 . -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 2563 - -pc_pfmg_skip_relax - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations 2564 when the underlying problem is isotropic, one of 0,1 2565 2566 Level: advanced 2567 2568 Notes: 2569 This is for CELL-centered descretizations 2570 2571 See `PCSYSPFMG` for a version suitable for systems of PDEs, and `PCSMG` 2572 2573 See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver 2574 2575 This must be used with the `MATHYPRESTRUCT` matrix type. 2576 2577 This provides only some of the functionality of PFMG, it supports only one block per process defined by a PETSc `DMDA`. 2578 2579 .seealso: [](ch_ksp), `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG` 2580 M*/ 2581 2582 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 2583 { 2584 PC_PFMG *ex; 2585 2586 PetscFunctionBegin; 2587 PetscCall(PetscNew(&ex)); 2588 pc->data = ex; 2589 2590 ex->its = 1; 2591 ex->tol = 1.e-8; 2592 ex->relax_type = 1; 2593 ex->rap_type = 0; 2594 ex->num_pre_relax = 1; 2595 ex->num_post_relax = 1; 2596 ex->max_levels = 0; 2597 ex->skip_relax = 0; 2598 ex->print_statistics = PETSC_FALSE; 2599 2600 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 2601 pc->ops->view = PCView_PFMG; 2602 pc->ops->destroy = PCDestroy_PFMG; 2603 pc->ops->apply = PCApply_PFMG; 2604 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 2605 pc->ops->setup = PCSetUp_PFMG; 2606 2607 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 2608 PetscHYPREInitialize(); 2609 PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver); 2610 PetscFunctionReturn(PETSC_SUCCESS); 2611 } 2612 2613 /* we know we are working with a HYPRE_SStructMatrix */ 2614 typedef struct { 2615 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 2616 HYPRE_SStructSolver ss_solver; 2617 2618 /* keep copy of SYSPFMG options used so may view them */ 2619 PetscInt its; 2620 PetscReal tol; 2621 PetscInt relax_type; 2622 PetscInt num_pre_relax, num_post_relax; 2623 } PC_SysPFMG; 2624 2625 static PetscErrorCode PCDestroy_SysPFMG(PC pc) 2626 { 2627 PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; 2628 2629 PetscFunctionBegin; 2630 if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver); 2631 PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 2632 PetscCall(PetscFree(pc->data)); 2633 PetscFunctionReturn(PETSC_SUCCESS); 2634 } 2635 2636 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi", "Red/Black-Gauss-Seidel"}; 2637 2638 static PetscErrorCode PCView_SysPFMG(PC pc, PetscViewer viewer) 2639 { 2640 PetscBool iascii; 2641 PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; 2642 2643 PetscFunctionBegin; 2644 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 2645 if (iascii) { 2646 PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SysPFMG preconditioning\n")); 2647 PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); 2648 PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); 2649 PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type])); 2650 PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); 2651 } 2652 PetscFunctionReturn(PETSC_SUCCESS); 2653 } 2654 2655 static PetscErrorCode PCSetFromOptions_SysPFMG(PC pc, PetscOptionItems *PetscOptionsObject) 2656 { 2657 PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; 2658 PetscBool flg = PETSC_FALSE; 2659 2660 PetscFunctionBegin; 2661 PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options"); 2662 PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL)); 2663 if (flg) PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel, ex->ss_solver, 3); 2664 PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL)); 2665 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, ex->ss_solver, ex->its); 2666 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)); 2667 PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax, ex->ss_solver, ex->num_pre_relax); 2668 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)); 2669 PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax, ex->ss_solver, ex->num_post_relax); 2670 2671 PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL)); 2672 PetscCallExternal(HYPRE_SStructSysPFMGSetTol, ex->ss_solver, ex->tol); 2673 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)); 2674 PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType, ex->ss_solver, ex->relax_type); 2675 PetscOptionsHeadEnd(); 2676 PetscFunctionReturn(PETSC_SUCCESS); 2677 } 2678 2679 static PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y) 2680 { 2681 PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; 2682 PetscScalar *yy; 2683 const PetscScalar *xx; 2684 PetscInt ilower[3], iupper[3]; 2685 HYPRE_Int hlower[3], hupper[3]; 2686 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data; 2687 PetscInt ordering = mx->dofs_order; 2688 PetscInt nvars = mx->nvars; 2689 PetscInt part = 0; 2690 PetscInt size; 2691 PetscInt i; 2692 2693 PetscFunctionBegin; 2694 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2695 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 2696 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2697 iupper[0] += ilower[0] - 1; 2698 iupper[1] += ilower[1] - 1; 2699 iupper[2] += ilower[2] - 1; 2700 hlower[0] = (HYPRE_Int)ilower[0]; 2701 hlower[1] = (HYPRE_Int)ilower[1]; 2702 hlower[2] = (HYPRE_Int)ilower[2]; 2703 hupper[0] = (HYPRE_Int)iupper[0]; 2704 hupper[1] = (HYPRE_Int)iupper[1]; 2705 hupper[2] = (HYPRE_Int)iupper[2]; 2706 2707 size = 1; 2708 for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1); 2709 2710 /* copy x values over to hypre for variable ordering */ 2711 if (ordering) { 2712 PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0); 2713 PetscCall(VecGetArrayRead(x, &xx)); 2714 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i))); 2715 PetscCall(VecRestoreArrayRead(x, &xx)); 2716 PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b); 2717 PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x); 2718 PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x); 2719 2720 /* copy solution values back to PETSc */ 2721 PetscCall(VecGetArray(y, &yy)); 2722 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i))); 2723 PetscCall(VecRestoreArray(y, &yy)); 2724 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2725 PetscScalar *z; 2726 PetscInt j, k; 2727 2728 PetscCall(PetscMalloc1(nvars * size, &z)); 2729 PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0); 2730 PetscCall(VecGetArrayRead(x, &xx)); 2731 2732 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2733 for (i = 0; i < size; i++) { 2734 k = i * nvars; 2735 for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j]; 2736 } 2737 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i))); 2738 PetscCall(VecRestoreArrayRead(x, &xx)); 2739 PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b); 2740 PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x); 2741 2742 /* copy solution values back to PETSc */ 2743 PetscCall(VecGetArray(y, &yy)); 2744 for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i))); 2745 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2746 for (i = 0; i < size; i++) { 2747 k = i * nvars; 2748 for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i]; 2749 } 2750 PetscCall(VecRestoreArray(y, &yy)); 2751 PetscCall(PetscFree(z)); 2752 } 2753 PetscFunctionReturn(PETSC_SUCCESS); 2754 } 2755 2756 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) 2757 { 2758 PC_SysPFMG *jac = (PC_SysPFMG *)pc->data; 2759 HYPRE_Int oits; 2760 2761 PetscFunctionBegin; 2762 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2763 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, its * jac->its); 2764 PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, rtol); 2765 PetscCall(PCApply_SysPFMG(pc, b, y)); 2766 PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations, jac->ss_solver, &oits); 2767 *outits = oits; 2768 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2769 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2770 PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, jac->tol); 2771 PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, jac->its); 2772 PetscFunctionReturn(PETSC_SUCCESS); 2773 } 2774 2775 static PetscErrorCode PCSetUp_SysPFMG(PC pc) 2776 { 2777 PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; 2778 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data; 2779 PetscBool flg; 2780 2781 PetscFunctionBegin; 2782 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESSTRUCT, &flg)); 2783 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESSTRUCT with this preconditioner"); 2784 2785 /* create the hypre sstruct solver object and set its information */ 2786 if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver); 2787 PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver); 2788 PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess, ex->ss_solver); 2789 PetscCallExternal(HYPRE_SStructSysPFMGSetup, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x); 2790 PetscFunctionReturn(PETSC_SUCCESS); 2791 } 2792 2793 /*MC 2794 PCSYSPFMG - the hypre SysPFMG multigrid solver 2795 2796 Level: advanced 2797 2798 Options Database Keys: 2799 + -pc_syspfmg_its <its> - number of iterations of SysPFMG to use as preconditioner 2800 . -pc_syspfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid 2801 . -pc_syspfmg_num_post_relax <steps> - number of smoothing steps after coarse grid 2802 . -pc_syspfmg_tol <tol> - tolerance of SysPFMG 2803 - -pc_syspfmg_relax_type <Weighted-Jacobi,Red/Black-Gauss-Seidel> - relaxation type for the up and down cycles 2804 2805 Notes: 2806 See `PCPFMG` for hypre's PFMG that works for a scalar PDE and `PCSMG` 2807 2808 See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver 2809 2810 This is for CELL-centered descretizations 2811 2812 This must be used with the `MATHYPRESSTRUCT` matrix type. 2813 2814 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`. 2815 2816 .seealso: [](ch_ksp), `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG` 2817 M*/ 2818 2819 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2820 { 2821 PC_SysPFMG *ex; 2822 2823 PetscFunctionBegin; 2824 PetscCall(PetscNew(&ex)); 2825 pc->data = ex; 2826 2827 ex->its = 1; 2828 ex->tol = 1.e-8; 2829 ex->relax_type = 1; 2830 ex->num_pre_relax = 1; 2831 ex->num_post_relax = 1; 2832 2833 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2834 pc->ops->view = PCView_SysPFMG; 2835 pc->ops->destroy = PCDestroy_SysPFMG; 2836 pc->ops->apply = PCApply_SysPFMG; 2837 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2838 pc->ops->setup = PCSetUp_SysPFMG; 2839 2840 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 2841 PetscHYPREInitialize(); 2842 PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver); 2843 PetscFunctionReturn(PETSC_SUCCESS); 2844 } 2845 2846 /* PC SMG */ 2847 typedef struct { 2848 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 2849 HYPRE_StructSolver hsolver; 2850 PetscInt its; /* keep copy of SMG options used so may view them */ 2851 PetscReal tol; 2852 PetscBool print_statistics; 2853 PetscInt num_pre_relax, num_post_relax; 2854 } PC_SMG; 2855 2856 static PetscErrorCode PCDestroy_SMG(PC pc) 2857 { 2858 PC_SMG *ex = (PC_SMG *)pc->data; 2859 2860 PetscFunctionBegin; 2861 if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver); 2862 PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 2863 PetscCall(PetscFree(pc->data)); 2864 PetscFunctionReturn(PETSC_SUCCESS); 2865 } 2866 2867 static PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer) 2868 { 2869 PetscBool iascii; 2870 PC_SMG *ex = (PC_SMG *)pc->data; 2871 2872 PetscFunctionBegin; 2873 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 2874 if (iascii) { 2875 PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SMG preconditioning\n")); 2876 PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); 2877 PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); 2878 PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); 2879 } 2880 PetscFunctionReturn(PETSC_SUCCESS); 2881 } 2882 2883 static PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems *PetscOptionsObject) 2884 { 2885 PC_SMG *ex = (PC_SMG *)pc->data; 2886 2887 PetscFunctionBegin; 2888 PetscOptionsHeadBegin(PetscOptionsObject, "SMG options"); 2889 2890 PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL)); 2891 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)); 2892 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)); 2893 PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL)); 2894 2895 PetscOptionsHeadEnd(); 2896 PetscFunctionReturn(PETSC_SUCCESS); 2897 } 2898 2899 static PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y) 2900 { 2901 PC_SMG *ex = (PC_SMG *)pc->data; 2902 PetscScalar *yy; 2903 const PetscScalar *xx; 2904 PetscInt ilower[3], iupper[3]; 2905 HYPRE_Int hlower[3], hupper[3]; 2906 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; 2907 2908 PetscFunctionBegin; 2909 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2910 PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); 2911 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2912 iupper[0] += ilower[0] - 1; 2913 iupper[1] += ilower[1] - 1; 2914 iupper[2] += ilower[2] - 1; 2915 hlower[0] = (HYPRE_Int)ilower[0]; 2916 hlower[1] = (HYPRE_Int)ilower[1]; 2917 hlower[2] = (HYPRE_Int)ilower[2]; 2918 hupper[0] = (HYPRE_Int)iupper[0]; 2919 hupper[1] = (HYPRE_Int)iupper[1]; 2920 hupper[2] = (HYPRE_Int)iupper[2]; 2921 2922 /* copy x values over to hypre */ 2923 PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0); 2924 PetscCall(VecGetArrayRead(x, &xx)); 2925 PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx); 2926 PetscCall(VecRestoreArrayRead(x, &xx)); 2927 PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb); 2928 PetscCallExternal(HYPRE_StructSMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx); 2929 2930 /* copy solution values back to PETSc */ 2931 PetscCall(VecGetArray(y, &yy)); 2932 PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy); 2933 PetscCall(VecRestoreArray(y, &yy)); 2934 PetscFunctionReturn(PETSC_SUCCESS); 2935 } 2936 2937 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) 2938 { 2939 PC_SMG *jac = (PC_SMG *)pc->data; 2940 HYPRE_Int oits; 2941 2942 PetscFunctionBegin; 2943 PetscCall(PetscCitationsRegister(hypreCitation, &cite)); 2944 PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, its * jac->its); 2945 PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, rtol); 2946 2947 PetscCall(PCApply_SMG(pc, b, y)); 2948 PetscCallExternal(HYPRE_StructSMGGetNumIterations, jac->hsolver, &oits); 2949 *outits = oits; 2950 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2951 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2952 PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, jac->tol); 2953 PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, jac->its); 2954 PetscFunctionReturn(PETSC_SUCCESS); 2955 } 2956 2957 static PetscErrorCode PCSetUp_SMG(PC pc) 2958 { 2959 PetscInt i, dim; 2960 PC_SMG *ex = (PC_SMG *)pc->data; 2961 Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; 2962 PetscBool flg; 2963 DMBoundaryType p[3]; 2964 PetscInt M[3]; 2965 2966 PetscFunctionBegin; 2967 PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg)); 2968 PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner"); 2969 2970 PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0)); 2971 // Check if power of 2 in periodic directions 2972 for (i = 0; i < dim; i++) { 2973 if (((M[i] & (M[i] - 1)) != 0) && (p[i] == DM_BOUNDARY_PERIODIC)) { 2974 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]); 2975 } 2976 } 2977 2978 /* create the hypre solver object and set its information */ 2979 if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, (ex->hsolver)); 2980 PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver); 2981 // The hypre options must be set here and not in SetFromOptions because it is created here! 2982 PetscCallExternal(HYPRE_StructSMGSetMaxIter, ex->hsolver, ex->its); 2983 PetscCallExternal(HYPRE_StructSMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax); 2984 PetscCallExternal(HYPRE_StructSMGSetNumPostRelax, ex->hsolver, ex->num_post_relax); 2985 PetscCallExternal(HYPRE_StructSMGSetTol, ex->hsolver, ex->tol); 2986 2987 PetscCallExternal(HYPRE_StructSMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx); 2988 PetscCallExternal(HYPRE_StructSMGSetZeroGuess, ex->hsolver); 2989 PetscFunctionReturn(PETSC_SUCCESS); 2990 } 2991 2992 /*MC 2993 PCSMG - the hypre (structured grid) SMG multigrid solver 2994 2995 Level: advanced 2996 2997 Options Database Keys: 2998 + -pc_smg_its <its> - number of iterations of SMG to use as preconditioner 2999 . -pc_smg_num_pre_relax <steps> - number of smoothing steps before coarse grid 3000 . -pc_smg_num_post_relax <steps> - number of smoothing steps after coarse grid 3001 - -pc_smg_tol <tol> - tolerance of SMG 3002 3003 Notes: 3004 This is for CELL-centered descretizations 3005 3006 This must be used with the `MATHYPRESTRUCT` `MatType`. 3007 3008 This does not provide all the functionality of hypre's SMG solver, it supports only one block per process defined by a PETSc `DMDA`. 3009 3010 See `PCSYSPFMG`, `PCSMG`, `PCPFMG`, and `PCHYPRE` for access to hypre's other preconditioners 3011 3012 .seealso: `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG` 3013 M*/ 3014 3015 PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc) 3016 { 3017 PC_SMG *ex; 3018 3019 PetscFunctionBegin; 3020 PetscCall(PetscNew(&ex)); 3021 pc->data = ex; 3022 3023 ex->its = 1; 3024 ex->tol = 1.e-8; 3025 ex->num_pre_relax = 1; 3026 ex->num_post_relax = 1; 3027 3028 pc->ops->setfromoptions = PCSetFromOptions_SMG; 3029 pc->ops->view = PCView_SMG; 3030 pc->ops->destroy = PCDestroy_SMG; 3031 pc->ops->apply = PCApply_SMG; 3032 pc->ops->applyrichardson = PCApplyRichardson_SMG; 3033 pc->ops->setup = PCSetUp_SMG; 3034 3035 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); 3036 PetscHYPREInitialize(); 3037 PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver); 3038 PetscFunctionReturn(PETSC_SUCCESS); 3039 } 3040