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