/* Provides an interface to the LLNL package hypre */ #include #include /*I "petscpc.h" I*/ /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */ #include #include #include <../src/vec/vec/impls/hypre/vhyp.h> #include <../src/mat/impls/hypre/mhypre.h> #include <../src/dm/impls/da/hypre/mhyp.h> #include <_hypre_parcsr_ls.h> #include #if defined(PETSC_HAVE_HYPRE_DEVICE) #include #endif static PetscBool cite = PETSC_FALSE; static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = " "{\\url{https://www.llnl.gov/casc/hypre}}\n}\n"; /* Private context (data structure) for the preconditioner. */ typedef struct { HYPRE_Solver hsolver; Mat hpmat; /* MatHYPRE */ HYPRE_Int (*destroy)(HYPRE_Solver); HYPRE_Int (*solve)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector); HYPRE_Int (*setup)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector); MPI_Comm comm_hypre; char *hypre_type; /* options for Pilut and BoomerAMG*/ PetscInt maxiter; PetscReal tol; /* options for Pilut */ PetscInt factorrowsize; /* options for ParaSails */ PetscInt nlevels; PetscReal threshold; PetscReal filter; PetscReal loadbal; PetscInt logging; PetscInt ruse; PetscInt symt; /* options for BoomerAMG */ PetscBool printstatistics; /* options for BoomerAMG */ PetscInt cycletype; PetscInt maxlevels; PetscReal strongthreshold; PetscReal maxrowsum; PetscInt gridsweeps[3]; PetscObjectParameterDeclare(PetscInt, coarsentype); PetscInt measuretype; PetscInt smoothtype; PetscInt smoothsweeps; PetscInt smoothnumlevels; PetscInt eu_level; /* Number of levels for ILU(k) in Euclid */ PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */ PetscInt eu_bj; /* Defines use of Block Jacobi ILU in Euclid */ PetscObjectParameterDeclare(PetscInt, relaxtype[3]); PetscReal relaxweight; PetscReal outerrelaxweight; PetscObjectParameterDeclare(PetscInt, relaxorder); PetscReal truncfactor; PetscBool applyrichardson; PetscInt pmax; PetscObjectParameterDeclare(PetscInt, interptype); PetscInt maxc; PetscInt minc; #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscObjectParameterDeclarePtr(const char, spgemm_type); // this is a global hypre parameter but is closely associated with BoomerAMG #endif /* GPU */ PetscObjectParameterDeclare(PetscBool3, keeptranspose); PetscInt rap2; PetscObjectParameterDeclare(PetscInt, mod_rap2); /* AIR */ PetscInt Rtype; PetscReal Rstrongthreshold; PetscReal Rfilterthreshold; PetscInt Adroptype; PetscReal Adroptol; PetscInt agg_nl; PetscObjectParameterDeclare(PetscInt, agg_interptype); PetscInt agg_num_paths; PetscBool nodal_relax; PetscInt nodal_relax_levels; PetscInt nodal_coarsening; PetscInt nodal_coarsening_diag; PetscInt vec_interp_variant; PetscInt vec_interp_qmax; PetscBool vec_interp_smooth; PetscInt interp_refine; /* NearNullSpace support */ VecHYPRE_IJVector *hmnull; HYPRE_ParVector *phmnull; PetscInt n_hmnull; Vec hmnull_constant; /* options for AS (Auxiliary Space preconditioners) */ PetscInt as_print; PetscInt as_max_iter; PetscReal as_tol; PetscInt as_relax_type; PetscInt as_relax_times; PetscReal as_relax_weight; PetscReal as_omega; PetscInt as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */ PetscReal as_amg_alpha_theta; /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */ PetscInt as_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */ PetscReal as_amg_beta_theta; /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS) */ PetscInt ams_cycle_type; PetscInt ads_cycle_type; /* additional data */ Mat G; /* MatHYPRE */ Mat C; /* MatHYPRE */ Mat alpha_Poisson; /* MatHYPRE */ Mat beta_Poisson; /* MatHYPRE */ /* extra information for AMS */ PetscInt dim; /* geometrical dimension */ VecHYPRE_IJVector coords[3]; VecHYPRE_IJVector constants[3]; VecHYPRE_IJVector interior; Mat RT_PiFull, RT_Pi[3]; Mat ND_PiFull, ND_Pi[3]; PetscBool ams_beta_is_zero; PetscBool ams_beta_is_zero_part; PetscInt ams_proj_freq; } PC_HYPRE; /* Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine. It is used in PCHMG. Other users should avoid using this function. */ static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc, PetscInt *nlevels, Mat *operators[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool same; PetscInt num_levels, l; Mat *mattmp; hypre_ParCSRMatrix **A_array; PetscFunctionBegin; PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same)); PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG"); num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver); PetscCall(PetscMalloc1(num_levels, &mattmp)); A_array = hypre_ParAMGDataAArray((hypre_ParAMGData *)jac->hsolver); for (l = 1; l < num_levels; l++) { PetscCall(MatCreateFromParCSR(A_array[l], MATAIJ, PETSC_OWN_POINTER, &mattmp[num_levels - 1 - l])); /* We want to own the data, and HYPRE can not touch this matrix any more */ A_array[l] = NULL; } *nlevels = num_levels; *operators = mattmp; PetscFunctionReturn(PETSC_SUCCESS); } /* Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine. It is used in PCHMG. Other users should avoid using this function. */ static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc, PetscInt *nlevels, Mat *interpolations[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool same; PetscInt num_levels, l; Mat *mattmp; hypre_ParCSRMatrix **P_array; PetscFunctionBegin; PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same)); PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG"); num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver); PetscCall(PetscMalloc1(num_levels, &mattmp)); P_array = hypre_ParAMGDataPArray((hypre_ParAMGData *)jac->hsolver); for (l = 1; l < num_levels; l++) { PetscCall(MatCreateFromParCSR(P_array[num_levels - 1 - l], MATAIJ, PETSC_OWN_POINTER, &mattmp[l - 1])); /* We want to own the data, and HYPRE can not touch this matrix any more */ P_array[num_levels - 1 - l] = NULL; } *nlevels = num_levels; *interpolations = mattmp; PetscFunctionReturn(PETSC_SUCCESS); } /* Boolean Vecs are created IN PLACE with using data from BoomerAMG. */ static PetscErrorCode PCHYPREGetCFMarkers_BoomerAMG(PC pc, PetscInt *n_per_level[], PetscBT *CFMarkers[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool same; PetscInt num_levels, fine_nodes = 0, coarse_nodes; PetscInt *n_per_temp; PetscBT *markertmp; hypre_IntArray **CF_marker_array; PetscFunctionBegin; PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same)); PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG"); num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver); PetscCall(PetscMalloc1(num_levels, &n_per_temp)); PetscCall(PetscMalloc1(num_levels - 1, &markertmp)); CF_marker_array = hypre_ParAMGDataCFMarkerArray((hypre_ParAMGData *)jac->hsolver); for (PetscInt l = 0, CFMaxIndex = num_levels - 2; CFMaxIndex >= 0; l++, CFMaxIndex--) { fine_nodes = hypre_IntArraySize(CF_marker_array[CFMaxIndex]); coarse_nodes = 0; PetscCall(PetscBTCreate(fine_nodes, &markertmp[l])); for (PetscInt k = 0; k < fine_nodes; k++) { if (hypre_IntArrayDataI(CF_marker_array[CFMaxIndex], k) > 0) { PetscCall(PetscBTSet(markertmp[l], k)); coarse_nodes++; } } n_per_temp[l] = coarse_nodes; } n_per_temp[num_levels - 1] = fine_nodes; *n_per_level = n_per_temp; *CFMarkers = markertmp; PetscFunctionReturn(PETSC_SUCCESS); } /* Resets (frees) Hypre's representation of the near null space */ static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscInt i; PetscFunctionBegin; for (i = 0; i < jac->n_hmnull; i++) PetscCall(VecHYPRE_IJVectorDestroy(&jac->hmnull[i])); PetscCall(PetscFree(jac->hmnull)); PetscCall(PetscFree(jac->phmnull)); PetscCall(VecDestroy(&jac->hmnull_constant)); jac->n_hmnull = 0; PetscFunctionReturn(PETSC_SUCCESS); } static const char *HYPRESpgemmTypes[] = {"cusparse", "hypre"}; static PetscErrorCode PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char name[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscFunctionBegin; jac->spgemm_type = name; PetscFunctionReturn(PETSC_SUCCESS); #endif } static PetscErrorCode PCSetUp_HYPRE(PC pc) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; Mat_HYPRE *hjac; HYPRE_ParCSRMatrix hmat; HYPRE_ParVector bv, xv; PetscBool ishypre; PetscFunctionBegin; /* default type is boomerAMG */ if (!jac->hypre_type) PetscCall(PCHYPRESetType(pc, "boomeramg")); /* get hypre matrix */ if (pc->flag == DIFFERENT_NONZERO_PATTERN) PetscCall(MatDestroy(&jac->hpmat)); PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRE, &ishypre)); if (!ishypre) { #if defined(PETSC_HAVE_HYPRE_DEVICE) && PETSC_PKG_HYPRE_VERSION_LE(2, 30, 0) /* Temporary fix since we do not support MAT_REUSE_MATRIX with HYPRE device */ PetscBool iscuda, iship, iskokkos; PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, "")); PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, "")); PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iskokkos, MATSEQAIJKOKKOS, MATMPIAIJKOKKOS, "")); if (iscuda || iship || iskokkos) PetscCall(MatDestroy(&jac->hpmat)); #endif PetscCall(MatConvert(pc->pmat, MATHYPRE, jac->hpmat ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &jac->hpmat)); } else { PetscCall(PetscObjectReference((PetscObject)pc->pmat)); PetscCall(MatDestroy(&jac->hpmat)); jac->hpmat = pc->pmat; } /* allow debug */ PetscCall(MatViewFromOptions(jac->hpmat, NULL, "-pc_hypre_mat_view")); hjac = (Mat_HYPRE *)jac->hpmat->data; /* special case for BoomerAMG */ if (jac->setup == HYPRE_BoomerAMGSetup) { MatNullSpace mnull; PetscBool has_const; PetscInt bs, nvec, i; PetscMemType memtype; const Vec *vecs; PetscCall(MatGetCurrentMemType(jac->hpmat, &memtype)); if (PetscMemTypeDevice(memtype)) { /* GPU defaults From https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options and /src/parcsr_ls/par_amg.c First handle options which users have interfaces for changing */ PetscObjectParameterSetDefault(jac, coarsentype, 8); PetscObjectParameterSetDefault(jac, relaxorder, 0); PetscObjectParameterSetDefault(jac, interptype, 6); PetscObjectParameterSetDefault(jac, relaxtype[0], 18); PetscObjectParameterSetDefault(jac, relaxtype[1], 18); #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscObjectParameterSetDefault(jac, spgemm_type, HYPRESpgemmTypes[0]); #endif #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_TRUE); PetscObjectParameterSetDefault(jac, mod_rap2, 1); #endif PetscObjectParameterSetDefault(jac, agg_interptype, 7); } else { PetscObjectParameterSetDefault(jac, coarsentype, 6); PetscObjectParameterSetDefault(jac, relaxorder, 1); PetscObjectParameterSetDefault(jac, interptype, 0); PetscObjectParameterSetDefault(jac, relaxtype[0], 6); PetscObjectParameterSetDefault(jac, relaxtype[1], 6); /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */ #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscObjectParameterSetDefault(jac, spgemm_type, "hypre"); #endif #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_FALSE); PetscObjectParameterSetDefault(jac, mod_rap2, 0); #endif PetscObjectParameterSetDefault(jac, agg_interptype, 4); } PetscCallHYPRE(HYPRE_BoomerAMGSetCycleType(jac->hsolver, (HYPRE_Int)jac->cycletype)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxLevels(jac->hsolver, (HYPRE_Int)jac->maxlevels)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter)); PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol)); PetscCallHYPRE(HYPRE_BoomerAMGSetTruncFactor(jac->hsolver, jac->truncfactor)); PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver, jac->strongthreshold)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver, jac->maxrowsum)); PetscCallHYPRE(HYPRE_BoomerAMGSetMeasureType(jac->hsolver, (HYPRE_Int)jac->measuretype)); PetscCallHYPRE(HYPRE_BoomerAMGSetAggNumLevels(jac->hsolver, (HYPRE_Int)jac->agg_nl)); PetscCallHYPRE(HYPRE_BoomerAMGSetPMaxElmts(jac->hsolver, (HYPRE_Int)jac->pmax)); PetscCallHYPRE(HYPRE_BoomerAMGSetNumPaths(jac->hsolver, (HYPRE_Int)jac->agg_num_paths)); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[0], 1)); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[1], 2)); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[2], 3)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxCoarseSize(jac->hsolver, (HYPRE_Int)jac->maxc)); PetscCallHYPRE(HYPRE_BoomerAMGSetMinCoarseSize(jac->hsolver, (HYPRE_Int)jac->minc)); PetscCallHYPRE(HYPRE_BoomerAMGSetCoarsenType(jac->hsolver, (HYPRE_Int)jac->coarsentype)); PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxOrder(jac->hsolver, (HYPRE_Int)jac->relaxorder)); PetscCallHYPRE(HYPRE_BoomerAMGSetInterpType(jac->hsolver, (HYPRE_Int)jac->interptype)); PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[0])); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[0], 1)); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[1], 2)); PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[2], 3)); /* GPU */ #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) { PetscBool flg_cusparse, flg_hypre; PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flg_cusparse)); PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flg_hypre)); if (flg_cusparse) PetscCallHYPRE(HYPRE_SetSpGemmUseCusparse(1)); else if (flg_hypre) PetscCallHYPRE(HYPRE_SetSpGemmUseCusparse(0)); else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEMM type %s; Choices are cusparse, hypre", jac->spgemm_type); } #endif #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) PetscCallHYPRE(HYPRE_BoomerAMGSetKeepTranspose(jac->hsolver, jac->keeptranspose == PETSC_BOOL3_TRUE ? 1 : 0)); PetscCallHYPRE(HYPRE_BoomerAMGSetRAP2(jac->hsolver, (HYPRE_Int)jac->rap2)); PetscCallHYPRE(HYPRE_BoomerAMGSetModuleRAP2(jac->hsolver, (HYPRE_Int)jac->mod_rap2)); #endif PetscCallHYPRE(HYPRE_BoomerAMGSetAggInterpType(jac->hsolver, (HYPRE_Int)jac->agg_interptype)); /* AIR */ #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) PetscCallHYPRE(HYPRE_BoomerAMGSetRestriction(jac->hsolver, (HYPRE_Int)jac->Rtype)); PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThresholdR(jac->hsolver, jac->Rstrongthreshold)); PetscCallHYPRE(HYPRE_BoomerAMGSetFilterThresholdR(jac->hsolver, jac->Rfilterthreshold)); PetscCallHYPRE(HYPRE_BoomerAMGSetADropTol(jac->hsolver, jac->Adroptol)); PetscCallHYPRE(HYPRE_BoomerAMGSetADropType(jac->hsolver, (HYPRE_Int)jac->Adroptype)); #endif PetscCall(MatGetBlockSize(pc->pmat, &bs)); if (bs > 1) PetscCallHYPRE(HYPRE_BoomerAMGSetNumFunctions(jac->hsolver, (HYPRE_Int)bs)); PetscCall(MatGetNearNullSpace(pc->mat, &mnull)); if (mnull) { PetscCall(PCHYPREResetNearNullSpace_Private(pc)); PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs)); PetscCall(PetscMalloc1(nvec + 1, &jac->hmnull)); PetscCall(PetscMalloc1(nvec + 1, &jac->phmnull)); for (i = 0; i < nvec; i++) { PetscCall(VecHYPRE_IJVectorCreate(vecs[i]->map, &jac->hmnull[i])); PetscCall(VecHYPRE_IJVectorCopy(vecs[i], jac->hmnull[i])); PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->hmnull[i]->ij, (void **)&jac->phmnull[i])); } if (has_const) { PetscCall(MatCreateVecs(pc->pmat, &jac->hmnull_constant, NULL)); PetscCall(VecSet(jac->hmnull_constant, 1)); PetscCall(VecNormalize(jac->hmnull_constant, NULL)); PetscCall(VecHYPRE_IJVectorCreate(jac->hmnull_constant->map, &jac->hmnull[nvec])); PetscCall(VecHYPRE_IJVectorCopy(jac->hmnull_constant, jac->hmnull[nvec])); PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->hmnull[nvec]->ij, (void **)&jac->phmnull[nvec])); nvec++; } PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVectors(jac->hsolver, (HYPRE_Int)nvec, jac->phmnull)); jac->n_hmnull = nvec; } } /* special case for AMS */ if (jac->setup == HYPRE_AMSSetup) { Mat_HYPRE *hm; HYPRE_ParCSRMatrix parcsr; PetscCheck(jac->coords[0] || jac->constants[0] || jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]), 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()"); if (jac->dim) PetscCallHYPRE(HYPRE_AMSSetDimension(jac->hsolver, (HYPRE_Int)jac->dim)); if (jac->constants[0]) { HYPRE_ParVector ozz, zoz, zzo = NULL; PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[0]->ij, (void **)(&ozz))); PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[1]->ij, (void **)(&zoz))); if (jac->constants[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[2]->ij, (void **)(&zzo))); PetscCallHYPRE(HYPRE_AMSSetEdgeConstantVectors(jac->hsolver, ozz, zoz, zzo)); } if (jac->coords[0]) { HYPRE_ParVector coords[3]; coords[0] = NULL; coords[1] = NULL; coords[2] = NULL; if (jac->coords[0]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[0]->ij, (void **)(&coords[0]))); if (jac->coords[1]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[1]->ij, (void **)(&coords[1]))); if (jac->coords[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[2]->ij, (void **)(&coords[2]))); PetscCallHYPRE(HYPRE_AMSSetCoordinateVectors(jac->hsolver, coords[0], coords[1], coords[2])); } PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient"); hm = (Mat_HYPRE *)jac->G->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr))); PetscCallHYPRE(HYPRE_AMSSetDiscreteGradient(jac->hsolver, parcsr)); if (jac->alpha_Poisson) { hm = (Mat_HYPRE *)jac->alpha_Poisson->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr))); PetscCallHYPRE(HYPRE_AMSSetAlphaPoissonMatrix(jac->hsolver, parcsr)); } if (jac->ams_beta_is_zero) { PetscCallHYPRE(HYPRE_AMSSetBetaPoissonMatrix(jac->hsolver, NULL)); } else if (jac->beta_Poisson) { hm = (Mat_HYPRE *)jac->beta_Poisson->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr))); PetscCallHYPRE(HYPRE_AMSSetBetaPoissonMatrix(jac->hsolver, parcsr)); } else if (jac->ams_beta_is_zero_part) { if (jac->interior) { HYPRE_ParVector interior = NULL; PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->interior->ij, (void **)(&interior))); PetscCallHYPRE(HYPRE_AMSSetInteriorNodes(jac->hsolver, interior)); } else { jac->ams_beta_is_zero_part = PETSC_FALSE; } } if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) { PetscInt i; HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3]; if (jac->ND_PiFull) { hm = (Mat_HYPRE *)jac->ND_PiFull->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsrfull))); } else { nd_parcsrfull = NULL; } for (i = 0; i < 3; ++i) { if (jac->ND_Pi[i]) { hm = (Mat_HYPRE *)jac->ND_Pi[i]->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsr[i]))); } else { nd_parcsr[i] = NULL; } } PetscCallHYPRE(HYPRE_AMSSetInterpolations(jac->hsolver, nd_parcsrfull, nd_parcsr[0], nd_parcsr[1], nd_parcsr[2])); } } /* special case for ADS */ if (jac->setup == HYPRE_ADSSetup) { Mat_HYPRE *hm; HYPRE_ParCSRMatrix parcsr; 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])))) { SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations"); } 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"); PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient"); PetscCheck(jac->C, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient"); if (jac->coords[0]) { HYPRE_ParVector coords[3]; coords[0] = NULL; coords[1] = NULL; coords[2] = NULL; if (jac->coords[0]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[0]->ij, (void **)(&coords[0]))); if (jac->coords[1]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[1]->ij, (void **)(&coords[1]))); if (jac->coords[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[2]->ij, (void **)(&coords[2]))); PetscCallHYPRE(HYPRE_ADSSetCoordinateVectors(jac->hsolver, coords[0], coords[1], coords[2])); } hm = (Mat_HYPRE *)jac->G->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr))); PetscCallHYPRE(HYPRE_ADSSetDiscreteGradient(jac->hsolver, parcsr)); hm = (Mat_HYPRE *)jac->C->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr))); PetscCallHYPRE(HYPRE_ADSSetDiscreteCurl(jac->hsolver, parcsr)); if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) { PetscInt i; HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3]; HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3]; if (jac->RT_PiFull) { hm = (Mat_HYPRE *)jac->RT_PiFull->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&rt_parcsrfull))); } else { rt_parcsrfull = NULL; } for (i = 0; i < 3; ++i) { if (jac->RT_Pi[i]) { hm = (Mat_HYPRE *)jac->RT_Pi[i]->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&rt_parcsr[i]))); } else { rt_parcsr[i] = NULL; } } if (jac->ND_PiFull) { hm = (Mat_HYPRE *)jac->ND_PiFull->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsrfull))); } else { nd_parcsrfull = NULL; } for (i = 0; i < 3; ++i) { if (jac->ND_Pi[i]) { hm = (Mat_HYPRE *)jac->ND_Pi[i]->data; PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsr[i]))); } else { nd_parcsr[i] = NULL; } } PetscCallHYPRE(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])); } } PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&bv)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&xv)); PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF)); PetscCallHYPRE((*jac->setup)(jac->hsolver, hmat, bv, xv)); PetscCall(PetscFPTrapPop()); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCApply_HYPRE(PC pc, Vec b, Vec x) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data; HYPRE_ParCSRMatrix hmat; HYPRE_ParVector jbv, jxv; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); if (!jac->applyrichardson) PetscCall(VecSet(x, 0.0)); PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b)); if (jac->applyrichardson) PetscCall(VecHYPRE_IJVectorPushVec(hjac->x, x)); else PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x)); PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&jbv)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&jxv)); PetscStackCallExternalVoid( "Hypre solve", do { HYPRE_Int hierr = (*jac->solve)(jac->hsolver, hmat, jbv, jxv); if (hierr) { PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr); HYPRE_ClearAllErrors(); } } while (0)); if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) PetscCallHYPRE(HYPRE_AMSProjectOutGradients(jac->hsolver, jxv)); PetscCall(VecHYPRE_IJVectorPopVec(hjac->x)); PetscCall(VecHYPRE_IJVectorPopVec(hjac->b)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCMatApply_HYPRE_BoomerAMG(PC pc, Mat B, Mat X) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data; hypre_ParCSRMatrix *par_matrix; HYPRE_ParVector hb, hx; const PetscScalar *b; PetscScalar *x; PetscInt m, N, lda; hypre_Vector *x_local; PetscMemType type; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&par_matrix)); PetscCall(MatGetLocalSize(B, &m, NULL)); PetscCall(MatGetSize(B, NULL, &N)); PetscCallHYPRE(HYPRE_ParMultiVectorCreate(hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), (HYPRE_Int)N, &hb)); PetscCallHYPRE(HYPRE_ParMultiVectorCreate(hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), (HYPRE_Int)N, &hx)); PetscCall(MatZeroEntries(X)); PetscCall(MatDenseGetArrayReadAndMemType(B, &b, &type)); PetscCall(MatDenseGetLDA(B, &lda)); 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); PetscCall(MatDenseGetLDA(X, &lda)); 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); x_local = hypre_ParVectorLocalVector(hb); PetscCallHYPRE(hypre_SeqVectorSetDataOwner(x_local, 0)); hypre_VectorData(x_local) = (HYPRE_Complex *)b; PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, NULL)); x_local = hypre_ParVectorLocalVector(hx); PetscCallHYPRE(hypre_SeqVectorSetDataOwner(x_local, 0)); hypre_VectorData(x_local) = (HYPRE_Complex *)x; PetscCallHYPRE(hypre_ParVectorInitialize_v2(hb, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE)); PetscCallHYPRE(hypre_ParVectorInitialize_v2(hx, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE)); PetscStackCallExternalVoid( "Hypre solve", do { HYPRE_Int hierr = (*jac->solve)(jac->hsolver, par_matrix, hb, hx); if (hierr) { PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr); HYPRE_ClearAllErrors(); } } while (0)); PetscCallHYPRE(HYPRE_ParVectorDestroy(hb)); PetscCallHYPRE(HYPRE_ParVectorDestroy(hx)); PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b)); PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCReset_HYPRE(PC pc) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; PetscCall(MatDestroy(&jac->hpmat)); PetscCall(MatDestroy(&jac->G)); PetscCall(MatDestroy(&jac->C)); PetscCall(MatDestroy(&jac->alpha_Poisson)); PetscCall(MatDestroy(&jac->beta_Poisson)); PetscCall(MatDestroy(&jac->RT_PiFull)); PetscCall(MatDestroy(&jac->RT_Pi[0])); PetscCall(MatDestroy(&jac->RT_Pi[1])); PetscCall(MatDestroy(&jac->RT_Pi[2])); PetscCall(MatDestroy(&jac->ND_PiFull)); PetscCall(MatDestroy(&jac->ND_Pi[0])); PetscCall(MatDestroy(&jac->ND_Pi[1])); PetscCall(MatDestroy(&jac->ND_Pi[2])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior)); PetscCall(PCHYPREResetNearNullSpace_Private(pc)); jac->ams_beta_is_zero = PETSC_FALSE; jac->ams_beta_is_zero_part = PETSC_FALSE; jac->dim = 0; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCDestroy_HYPRE(PC pc) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; PetscCall(PCReset_HYPRE(pc)); if (jac->destroy) PetscCallHYPRE((*jac->destroy)(jac->hsolver)); PetscCall(PetscFree(jac->hypre_type)); if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); PetscCall(PetscFree(pc->data)); PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetConstantEdgeVectors_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool flag; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options"); PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter)); PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetDropTolerance(jac->hsolver, jac->tol)); PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetFactorRowSize(jac->hsolver, (HYPRE_Int)jac->factorrowsize)); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_Pilut(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool isascii; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Pilut preconditioning\n")); if (jac->maxiter != PETSC_DEFAULT) { PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " default maximum number of iterations \n")); } if (jac->tol != PETSC_DEFAULT) { PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->tol)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " default drop tolerance \n")); } if (jac->factorrowsize != PETSC_DEFAULT) { PetscCall(PetscViewerASCIIPrintf(viewer, " factor row size %" PetscInt_FMT "\n", jac->factorrowsize)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " default factor row size \n")); } } PetscFunctionReturn(PETSC_SUCCESS); } static const char *HYPREILUType[] = { "Block-Jacobi-ILUk", "Block-Jacobi-ILUT", "", "", "", "", "", "", "", "", /* 0-9 */ "GMRES-ILUk", "GMRES-ILUT", "", "", "", "", "", "", "", "", /* 10-19 */ "NSH-ILUk", "NSH-ILUT", "", "", "", "", "", "", "", "", /* 20-29 */ "RAS-ILUk", "RAS-ILUT", "", "", "", "", "", "", "", "", /* 30-39 */ "ddPQ-GMRES-ILUk", "ddPQ-GMRES-ILUT", "", "", "", "", "", "", "", "", /* 40-49 */ "GMRES-ILU0" /* 50 */ }; static const char *HYPREILUIterSetup[] = {"default", "async-in-place", "async-explicit", "sync-explicit", "semisync-explicit"}; static PetscErrorCode PCSetFromOptions_HYPRE_ILU(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool flg; PetscInt indx; PetscReal tmpdbl; PetscBool tmp_truth; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ILU Options"); /* ILU: ILU Type */ PetscCall(PetscOptionsEList("-pc_hypre_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetType(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup type*/ PetscCall(PetscOptionsEList("-pc_hypre_ilu_iterative_setup_type", "Set ILU iterative setup type", "None", HYPREILUIterSetup, PETSC_STATIC_ARRAY_LENGTH(HYPREILUIterSetup), HYPREILUIterSetup[0], &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupType(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup option*/ PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupOption(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup maxiter */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupMaxIter(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup tolerance */ PetscCall(PetscOptionsReal("-pc_hypre_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupTolerance(jac->hsolver, tmpdbl)); /* ILU: ILU Print Level */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetPrintLevel(jac->hsolver, (HYPRE_Int)indx)); /* ILU: Logging */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetLogging(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU Level */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_level", "Set ILU level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetLevelOfFill(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU Max NNZ per row */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetMaxNnzPerRow(jac->hsolver, (HYPRE_Int)indx)); /* ILU: tolerance */ PetscCall(PetscOptionsReal("-pc_hypre_ilu_tol", "Tolerance for ILU", "None", 0, &tmpdbl, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetTol(jac->hsolver, tmpdbl)); /* ILU: maximum iteration count */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetMaxIter(jac->hsolver, (HYPRE_Int)indx)); /* ILU: drop threshold */ PetscCall(PetscOptionsReal("-pc_hypre_ilu_drop_threshold", "Drop threshold for ILU", "None", 0, &tmpdbl, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetDropThreshold(jac->hsolver, tmpdbl)); /* ILU: Triangular Solve */ PetscCall(PetscOptionsBool("-pc_hypre_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetTriSolve(jac->hsolver, tmp_truth)); /* ILU: Lower Jacobi iteration */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetLowerJacobiIters(jac->hsolver, (HYPRE_Int)indx)); /* ILU: Upper Jacobi iteration */ PetscCall(PetscOptionsInt("-pc_hypre_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetUpperJacobiIters(jac->hsolver, (HYPRE_Int)indx)); /* ILU: local reordering */ PetscCall(PetscOptionsBool("-pc_hypre_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg) PetscCallHYPRE(HYPRE_ILUSetLocalReordering(jac->hsolver, tmp_truth)); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_ILU(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; hypre_ParILUData *ilu_data = (hypre_ParILUData *)jac->hsolver; PetscBool isascii; PetscInt indx; PetscReal tmpdbl; PetscReal *tmpdbl3; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ILU preconditioning\n")); PetscStackCallExternalVoid("hypre_ParILUDataIluType", indx = hypre_ParILUDataIluType(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx)); PetscStackCallExternalVoid("hypre_ParILUDataLfil", indx = hypre_ParILUDataLfil(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataMaxIter", indx = hypre_ParILUDataMaxIter(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataMaxRowNnz", indx = hypre_ParILUDataMaxRowNnz(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataTriSolve", indx = hypre_ParILUDataTriSolve(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataTol", tmpdbl = hypre_ParILUDataTol(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", tmpdbl)); PetscStackCallExternalVoid("hypre_ParILUDataDroptol", tmpdbl3 = hypre_ParILUDataDroptol(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e / %e / %e\n", tmpdbl3[0], tmpdbl3[1], tmpdbl3[2])); PetscStackCallExternalVoid("hypre_ParILUDataReorderingType", indx = hypre_ParILUDataReorderingType(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataLowerJacobiIters", indx = hypre_ParILUDataLowerJacobiIters(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataUpperJacobiIters", indx = hypre_ParILUDataUpperJacobiIters(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataPrintLevel", indx = hypre_ParILUDataPrintLevel(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataLogging", indx = hypre_ParILUDataLogging(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupType", indx = hypre_ParILUDataIterativeSetupType(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx)); PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupOption", indx = hypre_ParILUDataIterativeSetupOption(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupMaxIter", indx = hypre_ParILUDataIterativeSetupMaxIter(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupTolerance", tmpdbl = hypre_ParILUDataIterativeSetupTolerance(ilu_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", tmpdbl)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool flag, eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Euclid Options"); PetscCall(PetscOptionsInt("-pc_hypre_euclid_level", "Factorization levels", "None", jac->eu_level, &jac->eu_level, &flag)); if (flag) PetscCallHYPRE(HYPRE_EuclidSetLevel(jac->hsolver, (HYPRE_Int)jac->eu_level)); PetscCall(PetscOptionsReal("-pc_hypre_euclid_droptolerance", "Drop tolerance for ILU(k) in Euclid", "None", jac->eu_droptolerance, &jac->eu_droptolerance, &flag)); if (flag) { PetscMPIInt size; PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size)); PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "hypre's Euclid does not support a parallel drop tolerance"); PetscCallHYPRE(HYPRE_EuclidSetILUT(jac->hsolver, jac->eu_droptolerance)); } PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag)); if (flag) { jac->eu_bj = eu_bj ? 1 : 0; PetscCallHYPRE(HYPRE_EuclidSetBJ(jac->hsolver, (HYPRE_Int)jac->eu_bj)); } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_Euclid(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool isascii; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Euclid preconditioning\n")); if (jac->eu_level != PETSC_DEFAULT) { PetscCall(PetscViewerASCIIPrintf(viewer, " factorization levels %" PetscInt_FMT "\n", jac->eu_level)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " default factorization levels \n")); } PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->eu_droptolerance)); PetscCall(PetscViewerASCIIPrintf(viewer, " use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc, Vec b, Vec x) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data; HYPRE_ParCSRMatrix hmat; HYPRE_ParVector jbv, jxv; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCall(VecSet(x, 0.0)); PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b)); PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x)); PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&jbv)); PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&jxv)); PetscStackCallExternalVoid( "Hypre Transpose solve", do { HYPRE_Int hierr = HYPRE_BoomerAMGSolveT(jac->hsolver, hmat, jbv, jxv); if (hierr) { /* error code of 1 in BoomerAMG merely means convergence not achieved */ PetscCheck(hierr == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr); HYPRE_ClearAllErrors(); } } while (0)); PetscCall(VecHYPRE_IJVectorPopVec(hjac->x)); PetscCall(VecHYPRE_IJVectorPopVec(hjac->b)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char *spgemm[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) *spgemm = jac->spgemm_type; #endif PetscFunctionReturn(PETSC_SUCCESS); } static const char *HYPREBoomerAMGCycleType[] = {"", "V", "W"}; static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP", "Ruge-Stueben", "", "modifiedRuge-Stueben", "", "", "Falgout", "", "PMIS", "", "HMIS"}; static const char *HYPREBoomerAMGMeasureType[] = {"local", "global"}; /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */ static const char *HYPREBoomerAMGSmoothType[] = {"ILU", "Schwarz-smoothers", "Pilut", "ParaSails", "Euclid"}; 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"}; 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"}; static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscInt bs, n, indx, level; PetscBool flg, tmp_truth; PetscReal tmpdbl, twodbl[2]; const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"}; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options"); PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg)); if (flg) { jac->cycletype = indx + 1; PetscCallHYPRE(HYPRE_BoomerAMGSetCycleType(jac->hsolver, (HYPRE_Int)jac->cycletype)); } PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg, 2)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxLevels(jac->hsolver, (HYPRE_Int)jac->maxlevels)); PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg, 1)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter)); 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol)); bs = 1; if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs)); PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNumFunctions(jac->hsolver, (HYPRE_Int)bs)); PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_truncfactor", "Truncation factor for interpolation (0=no truncation)", "None", jac->truncfactor, &jac->truncfactor, &flg, 0.0)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetTruncFactor(jac->hsolver, jac->truncfactor)); PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_P_max", "Max elements per row for interpolation operator (0=unlimited)", "None", jac->pmax, &jac->pmax, &flg, 0)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetPMaxElmts(jac->hsolver, (HYPRE_Int)jac->pmax)); PetscCall(PetscOptionsRangeInt("-pc_hypre_boomeramg_agg_nl", "Number of levels of aggressive coarsening", "None", jac->agg_nl, &jac->agg_nl, &flg, 0, jac->maxlevels)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetAggNumLevels(jac->hsolver, (HYPRE_Int)jac->agg_nl)); 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNumPaths(jac->hsolver, (HYPRE_Int)jac->agg_num_paths)); PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_strong_threshold", "Threshold for being strongly connected", "None", jac->strongthreshold, &jac->strongthreshold, &flg, 0.0)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver, jac->strongthreshold)); PetscCall(PetscOptionsRangeReal("-pc_hypre_boomeramg_max_row_sum", "Maximum row sum", "None", jac->maxrowsum, &jac->maxrowsum, &flg, 0.0, 1.0)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver, jac->maxrowsum)); /* Grid sweeps */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg)); if (flg) { /* modify the jac structure so we can view the updated options with PC_View */ jac->gridsweeps[0] = indx; jac->gridsweeps[1] = indx; /*defaults coarse to 1 */ jac->gridsweeps[2] = 1; } PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNodal(jac->hsolver, (HYPRE_Int)jac->nodal_coarsening)); 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNodalDiag(jac->hsolver, (HYPRE_Int)jac->nodal_coarsening_diag)); PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVecVariant(jac->hsolver, (HYPRE_Int)jac->vec_interp_variant)); 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVecQMax(jac->hsolver, (HYPRE_Int)jac->vec_interp_qmax)); 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothInterpVectors(jac->hsolver, jac->vec_interp_smooth)); PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpRefine(jac->hsolver, (HYPRE_Int)jac->interp_refine)); PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg)); if (flg) { PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 1)); jac->gridsweeps[0] = indx; } PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg)); if (flg) { PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 2)); jac->gridsweeps[1] = indx; } PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg)); if (flg) { PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 3)); jac->gridsweeps[2] = indx; } /* Smooth type */ PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_smooth_type", "Enable more complex smoothers", "None", HYPREBoomerAMGSmoothType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGSmoothType), HYPREBoomerAMGSmoothType[0], &indx, &flg)); if (flg) { jac->smoothtype = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothType(jac->hsolver, (HYPRE_Int)indx + 5)); jac->smoothnumlevels = 25; PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, 25)); } /* Number of smoothing levels */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels", "Number of levels on which more complex smoothers are used", "None", 25, &indx, &flg)); if (flg && (jac->smoothtype != -1)) { jac->smoothnumlevels = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, (HYPRE_Int)indx)); } /* Smooth num sweeps */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_sweeps", "Set number of smoother sweeps", "None", 1, &indx, &flg)); if (flg && indx > 0) { jac->smoothsweeps = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumSweeps(jac->hsolver, (HYPRE_Int)indx)); } /* ILU: ILU Type */ PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUType(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup type*/ 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)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupType(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup option*/ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupOption(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup maxiter */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupMaxIter(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU iterative setup tolerance */ PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg)); if (flg) PetscCallHYPRE(hypre_BoomerAMGSetILUIterSetupTolerance(jac->hsolver, tmpdbl)); /* ILU: ILU Print Level */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetPrintLevel(jac->hsolver, (HYPRE_Int)indx)); /* ILU: Logging */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetLogging(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU Level */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_level", "Set ILU level", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULevel(jac->hsolver, (HYPRE_Int)indx)); /* ILU: ILU Max NNZ per row */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUMaxRowNnz(jac->hsolver, (HYPRE_Int)indx)); /* ILU: maximum iteration count */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUMaxIter(jac->hsolver, (HYPRE_Int)indx)); /* ILU: drop threshold */ PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_drop_tol", "Drop tolerance for ILU", "None", 0, &tmpdbl, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUDroptol(jac->hsolver, tmpdbl)); /* ILU: Triangular Solve */ PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUTriSolve(jac->hsolver, tmp_truth)); /* ILU: Lower Jacobi iteration */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULowerJacobiIters(jac->hsolver, (HYPRE_Int)indx)); /* ILU: Upper Jacobi iteration */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUUpperJacobiIters(jac->hsolver, (HYPRE_Int)indx)); /* ILU: local reordering */ PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULocalReordering(jac->hsolver, tmp_truth)); /* Number of levels for ILU(k) for Euclid */ PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_eu_level", "Number of levels for ILU(k) in Euclid smoother", "None", 0, &indx, &flg)); if (flg && (jac->smoothtype == 4)) { jac->eu_level = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetEuLevel(jac->hsolver, (HYPRE_Int)indx)); } /* Filter for ILU(k) for Euclid */ PetscReal droptolerance; PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance", "Drop tolerance for ILU(k) in Euclid smoother", "None", 0, &droptolerance, &flg)); if (flg && (jac->smoothtype == 4)) { jac->eu_droptolerance = droptolerance; PetscCallHYPRE(HYPRE_BoomerAMGSetEuLevel(jac->hsolver, droptolerance)); } /* Use Block Jacobi ILUT for Euclid */ PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg && (jac->smoothtype == 4)) { jac->eu_bj = tmp_truth; PetscCallHYPRE(HYPRE_BoomerAMGSetEuBJ(jac->hsolver, (HYPRE_Int)jac->eu_bj)); } /* Relax type */ PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all", "Relax type for the up and down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]], &indx, &flg)); if (flg) jac->relaxtype[0] = jac->relaxtype[1] = indx; PetscCall( PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down", "Relax type for the down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]], &indx, &flg)); if (flg) jac->relaxtype[0] = indx; PetscCall( PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up", "Relax type for the up cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), jac->relaxtype[1] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[1]], &indx, &flg)); if (flg) jac->relaxtype[1] = indx; PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[jac->relaxtype[2]], &indx, &flg)); if (flg) jac->relaxtype[2] = indx; /* Relaxation Weight */ 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)); if (flg) { PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxWt(jac->hsolver, tmpdbl)); jac->relaxweight = tmpdbl; } n = 2; twodbl[0] = twodbl[1] = 1.0; PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level", "Set the relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg)); if (flg) { 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); indx = (int)PetscAbsReal(twodbl[1]); PetscCallHYPRE(HYPRE_BoomerAMGSetLevelRelaxWt(jac->hsolver, twodbl[0], (HYPRE_Int)indx)); } /* Outer relaxation Weight */ 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)); if (flg) { PetscCallHYPRE(HYPRE_BoomerAMGSetOuterWt(jac->hsolver, tmpdbl)); jac->outerrelaxweight = tmpdbl; } n = 2; twodbl[0] = twodbl[1] = 1.0; PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg)); if (flg) { 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); indx = (int)PetscAbsReal(twodbl[1]); PetscCallHYPRE(HYPRE_BoomerAMGSetLevelOuterWt(jac->hsolver, twodbl[0], (HYPRE_Int)indx)); } /* the Relax Order */ PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg) jac->relaxorder = !tmp_truth; PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg)); if (flg) { jac->measuretype = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetMeasureType(jac->hsolver, (HYPRE_Int)jac->measuretype)); } PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), jac->coarsentype < 0 ? "unknown" : HYPREBoomerAMGCoarsenType[jac->coarsentype], &indx, &flg)); if (flg) jac->coarsentype = indx; PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxCoarseSize(jac->hsolver, (HYPRE_Int)jac->maxc)); PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg)); if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMinCoarseSize(jac->hsolver, (HYPRE_Int)jac->minc)); #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) // global parameter but is closely associated with BoomerAMG 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), jac->spgemm_type, &indx, &flg)); if (flg) PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, HYPRESpgemmTypes[indx])); #endif /* AIR */ #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) 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)); PetscCallHYPRE(HYPRE_BoomerAMGSetRestriction(jac->hsolver, (HYPRE_Int)jac->Rtype)); if (jac->Rtype) { HYPRE_Int **grid_relax_points = hypre_TAlloc(HYPRE_Int *, 4, HYPRE_MEMORY_HOST); char *prerelax[256]; char *postrelax[256]; char stringF[2] = "F", stringC[2] = "C", stringA[2] = "A"; PetscInt ns_down = 256, ns_up = 256; PetscBool matchF, matchC, matchA; 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 */ PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR", "Threshold for R", "None", jac->Rstrongthreshold, &jac->Rstrongthreshold, NULL)); PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThresholdR(jac->hsolver, jac->Rstrongthreshold)); PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR", "Filter threshold for R", "None", jac->Rfilterthreshold, &jac->Rfilterthreshold, NULL)); PetscCallHYPRE(HYPRE_BoomerAMGSetFilterThresholdR(jac->hsolver, jac->Rfilterthreshold)); 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)); PetscCallHYPRE(HYPRE_BoomerAMGSetADropTol(jac->hsolver, (HYPRE_Int)jac->Adroptol)); 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)); PetscCallHYPRE(HYPRE_BoomerAMGSetADropType(jac->hsolver, (HYPRE_Int)jac->Adroptype)); PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_prerelax", "Defines prerelax scheme", "None", prerelax, &ns_down, NULL)); PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_postrelax", "Defines postrelax scheme", "None", postrelax, &ns_up, NULL)); 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"); 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"); grid_relax_points[0] = NULL; grid_relax_points[1] = hypre_TAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST); grid_relax_points[2] = hypre_TAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST); grid_relax_points[3] = hypre_TAlloc(HYPRE_Int, jac->gridsweeps[2], HYPRE_MEMORY_HOST); grid_relax_points[3][0] = 0; // set down relax scheme for (PetscInt i = 0; i < ns_down; i++) { PetscCall(PetscStrcasecmp(prerelax[i], stringF, &matchF)); PetscCall(PetscStrcasecmp(prerelax[i], stringC, &matchC)); PetscCall(PetscStrcasecmp(prerelax[i], stringA, &matchA)); PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_prerelax are C, F, and A"); if (matchF) grid_relax_points[1][i] = -1; else if (matchC) grid_relax_points[1][i] = 1; else if (matchA) grid_relax_points[1][i] = 0; } // set up relax scheme for (PetscInt i = 0; i < ns_up; i++) { PetscCall(PetscStrcasecmp(postrelax[i], stringF, &matchF)); PetscCall(PetscStrcasecmp(postrelax[i], stringC, &matchC)); PetscCall(PetscStrcasecmp(postrelax[i], stringA, &matchA)); PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_postrelax are C, F, and A"); if (matchF) grid_relax_points[2][i] = -1; else if (matchC) grid_relax_points[2][i] = 1; else if (matchA) grid_relax_points[2][i] = 0; } // set coarse relax scheme for (PetscInt i = 0; i < jac->gridsweeps[2]; i++) grid_relax_points[3][i] = 0; // Pass relax schemes to hypre PetscCallHYPRE(HYPRE_BoomerAMGSetGridRelaxPoints(jac->hsolver, grid_relax_points)); // cleanup memory for (PetscInt i = 0; i < ns_down; i++) PetscCall(PetscFree(prerelax[i])); for (PetscInt i = 0; i < ns_up; i++) PetscCall(PetscFree(postrelax[i])); } #endif #if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9) 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); #endif PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_interp_type", "Interpolation type", "None", HYPREBoomerAMGInterpType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGInterpType), jac->interptype < 0 ? "unknown" : HYPREBoomerAMGInterpType[jac->interptype], &indx, &flg)); if (flg) jac->interptype = indx; PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", &flg)); if (flg) { level = 3; PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", level, &level, NULL)); jac->printstatistics = PETSC_TRUE; PetscCallHYPRE(HYPRE_BoomerAMGSetPrintLevel(jac->hsolver, (HYPRE_Int)level)); } PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", &flg)); if (flg) { level = 3; PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", level, &level, NULL)); jac->printstatistics = PETSC_TRUE; PetscCallHYPRE(HYPRE_BoomerAMGSetDebugFlag(jac->hsolver, (HYPRE_Int)level)); } PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg)); if (flg && tmp_truth) { PetscInt tmp_int; PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", (HYPRE_Int)jac->nodal_relax_levels, &tmp_int, &flg)); if (flg) jac->nodal_relax_levels = tmp_int; PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothType(jac->hsolver, 6)); PetscCallHYPRE(HYPRE_BoomerAMGSetDomainType(jac->hsolver, 1)); PetscCallHYPRE(HYPRE_BoomerAMGSetOverlap(jac->hsolver, 0)); PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, (HYPRE_Int)jac->nodal_relax_levels)); } PetscCall(PetscOptionsBool3("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL)); /* options for ParaSails solvers */ PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flg)); if (flg) { jac->symt = indx; PetscCallHYPRE(HYPRE_BoomerAMGSetSym(jac->hsolver, (HYPRE_Int)jac->symt)); } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } 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) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; HYPRE_Int oits; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->maxiter))); PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, rtol)); jac->applyrichardson = PETSC_TRUE; PetscCall(PCApply_HYPRE(pc, b, y)); jac->applyrichardson = PETSC_FALSE; PetscCallHYPRE(HYPRE_BoomerAMGGetNumIterations(jac->hsolver, &oits)); *outits = oits; if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; else *reason = PCRICHARDSON_CONVERGED_RTOL; PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol)); PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; hypre_ParAMGData *amg_data = (hypre_ParAMGData *)jac->hsolver; PetscBool isascii; PetscInt indx; PetscReal val; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE BoomerAMG preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " Cycle type %s\n", HYPREBoomerAMGCycleType[jac->cycletype])); PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of levels %" PetscInt_FMT "\n", jac->maxlevels)); PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of iterations PER hypre call %" PetscInt_FMT "\n", jac->maxiter)); PetscCall(PetscViewerASCIIPrintf(viewer, " Convergence tolerance PER hypre call %g\n", (double)jac->tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for strong coupling %g\n", (double)jac->strongthreshold)); PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation truncation factor %g\n", (double)jac->truncfactor)); PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: max elements per row %" PetscInt_FMT "\n", jac->pmax)); if (jac->interp_refine) PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: number of steps of weighted refinement %" PetscInt_FMT "\n", jac->interp_refine)); PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", jac->agg_nl)); PetscCall(PetscViewerASCIIPrintf(viewer, " Number of paths for aggressive coarsening %" PetscInt_FMT "\n", jac->agg_num_paths)); PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum row sums %g\n", (double)jac->maxrowsum)); PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps down %" PetscInt_FMT "\n", jac->gridsweeps[0])); PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps up %" PetscInt_FMT "\n", jac->gridsweeps[1])); PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps on coarse %" PetscInt_FMT "\n", jac->gridsweeps[2])); PetscCall(PetscViewerASCIIPrintf(viewer, " Relax down %s\n", jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]])); PetscCall(PetscViewerASCIIPrintf(viewer, " Relax up %s\n", jac->relaxtype[1] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[1]])); PetscCall(PetscViewerASCIIPrintf(viewer, " Relax on coarse %s\n", HYPREBoomerAMGRelaxType[jac->relaxtype[2]])); PetscCall(PetscViewerASCIIPrintf(viewer, " Relax weight (all) %g\n", (double)jac->relaxweight)); PetscCall(PetscViewerASCIIPrintf(viewer, " Outer relax weight (all) %g\n", (double)jac->outerrelaxweight)); PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum size of coarsest grid %" PetscInt_FMT "\n", jac->maxc)); PetscCall(PetscViewerASCIIPrintf(viewer, " Minimum size of coarsest grid %" PetscInt_FMT "\n", jac->minc)); if (jac->relaxorder == PETSC_DECIDE) { PetscCall(PetscViewerASCIIPrintf(viewer, " CF-relaxation option not yet determined\n")); } else if (jac->relaxorder) { PetscCall(PetscViewerASCIIPrintf(viewer, " Using CF-relaxation\n")); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " Not using CF-relaxation\n")); } if (jac->smoothtype != -1) { PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth type %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype])); PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num levels %" PetscInt_FMT "\n", jac->smoothnumlevels)); PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num sweeps %" PetscInt_FMT "\n", jac->smoothsweeps)); if (jac->smoothtype == 0) { PetscStackCallExternalVoid("hypre_ParAMGDataILUType", indx = hypre_ParAMGDataILUType(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILULevel", indx = hypre_ParAMGDataILULevel(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUMaxIter", indx = hypre_ParAMGDataILUMaxIter(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUMaxRowNnz", indx = hypre_ParAMGDataILUMaxRowNnz(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUTriSolve", indx = hypre_ParAMGDataILUTriSolve(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataTol", val = hypre_ParAMGDataTol(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", val)); PetscStackCallExternalVoid("hypre_ParAMGDataILUDroptol", val = hypre_ParAMGDataILUDroptol(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e\n", val)); PetscStackCallExternalVoid("hypre_ParAMGDataILULocalReordering", indx = hypre_ParAMGDataILULocalReordering(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILULowerJacobiIters", indx = hypre_ParAMGDataILULowerJacobiIters(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUUpperJacobiIters", indx = hypre_ParAMGDataILUUpperJacobiIters(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataPrintLevel", indx = hypre_ParAMGDataPrintLevel(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataLogging", indx = hypre_ParAMGDataLogging(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupType", indx = hypre_ParAMGDataILUIterSetupType(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupOption", indx = hypre_ParAMGDataILUIterSetupOption(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupMaxIter", indx = hypre_ParAMGDataILUIterSetupMaxIter(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx)); PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupTolerance", val = hypre_ParAMGDataILUIterSetupTolerance(amg_data)); PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", val)); } } else { PetscCall(PetscViewerASCIIPrintf(viewer, " Not using more complex smoothers.\n")); } if (jac->smoothtype == 3) { PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level)); PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance)); PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj)); } PetscCall(PetscViewerASCIIPrintf(viewer, " Measure type %s\n", HYPREBoomerAMGMeasureType[jac->measuretype])); PetscCall(PetscViewerASCIIPrintf(viewer, " Coarsen type %s\n", jac->coarsentype < 0 ? "not yet set" : HYPREBoomerAMGCoarsenType[jac->coarsentype])); PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation type %s\n", jac->interptype != 100 ? (jac->interptype < 0 ? "not yet set" : HYPREBoomerAMGInterpType[jac->interptype]) : "1pt")); if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening)); if (jac->vec_interp_variant) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant)); PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax)); PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth)); } if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels)); #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", jac->spgemm_type)); #else PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", "hypre")); #endif /* AIR */ if (jac->Rtype) { PetscCall(PetscViewerASCIIPrintf(viewer, " Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype)); PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for R %g\n", (double)jac->Rstrongthreshold)); PetscCall(PetscViewerASCIIPrintf(viewer, " Filter for R %g\n", (double)jac->Rfilterthreshold)); PetscCall(PetscViewerASCIIPrintf(viewer, " A drop tolerance %g\n", (double)jac->Adroptol)); PetscCall(PetscViewerASCIIPrintf(viewer, " A drop type %" PetscInt_FMT "\n", jac->Adroptype)); } } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscInt indx; PetscBool flag; const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"}; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ParaSails Options"); PetscCall(PetscOptionsInt("-pc_hypre_parasails_nlevels", "Number of number of levels", "None", jac->nlevels, &jac->nlevels, 0)); PetscCall(PetscOptionsReal("-pc_hypre_parasails_thresh", "Threshold", "None", jac->threshold, &jac->threshold, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetParams(jac->hsolver, jac->threshold, (HYPRE_Int)jac->nlevels)); PetscCall(PetscOptionsReal("-pc_hypre_parasails_filter", "filter", "None", jac->filter, &jac->filter, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetFilter(jac->hsolver, jac->filter)); PetscCall(PetscOptionsReal("-pc_hypre_parasails_loadbal", "Load balance", "None", jac->loadbal, &jac->loadbal, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetLoadbal(jac->hsolver, (HYPRE_Int)jac->loadbal)); PetscCall(PetscOptionsBool("-pc_hypre_parasails_logging", "Print info to screen", "None", (PetscBool)jac->logging, (PetscBool *)&jac->logging, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetLogging(jac->hsolver, (HYPRE_Int)jac->logging)); PetscCall(PetscOptionsBool("-pc_hypre_parasails_reuse", "Reuse nonzero pattern in preconditioner", "None", (PetscBool)jac->ruse, (PetscBool *)&jac->ruse, &flag)); if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetReuse(jac->hsolver, (HYPRE_Int)jac->ruse)); PetscCall(PetscOptionsEList("-pc_hypre_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flag)); if (flag) { jac->symt = indx; PetscCallHYPRE(HYPRE_ParaSailsSetSym(jac->hsolver, (HYPRE_Int)jac->symt)); } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_ParaSails(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool isascii; const char *symt = 0; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ParaSails preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " nlevels %" PetscInt_FMT "\n", jac->nlevels)); PetscCall(PetscViewerASCIIPrintf(viewer, " threshold %g\n", (double)jac->threshold)); PetscCall(PetscViewerASCIIPrintf(viewer, " filter %g\n", (double)jac->filter)); PetscCall(PetscViewerASCIIPrintf(viewer, " load balance %g\n", (double)jac->loadbal)); PetscCall(PetscViewerASCIIPrintf(viewer, " reuse nonzero structure %s\n", PetscBools[jac->ruse])); PetscCall(PetscViewerASCIIPrintf(viewer, " print info to screen %s\n", PetscBools[jac->logging])); if (!jac->symt) symt = "nonsymmetric matrix and preconditioner"; else if (jac->symt == 1) symt = "SPD matrix and preconditioner"; else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner"; else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt); PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", symt)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscInt n; PetscBool flag, flag2, flag3, flag4; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE AMS Options"); PetscCall(PetscOptionsInt("-pc_hypre_ams_print_level", "Debugging output level for AMS", "None", jac->as_print, &jac->as_print, &flag)); if (flag) PetscCallHYPRE(HYPRE_AMSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print)); 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)); if (flag) PetscCallHYPRE(HYPRE_AMSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter)); PetscCall(PetscOptionsInt("-pc_hypre_ams_cycle_type", "Cycle type for AMS multigrid", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag)); if (flag) PetscCallHYPRE(HYPRE_AMSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type)); PetscCall(PetscOptionsReal("-pc_hypre_ams_tol", "Error tolerance for AMS multigrid", "None", jac->as_tol, &jac->as_tol, &flag)); if (flag) PetscCallHYPRE(HYPRE_AMSSetTol(jac->hsolver, jac->as_tol)); PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_type", "Relaxation type for AMS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag)); PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_times", "Number of relaxation steps for AMS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2)); PetscCall(PetscOptionsReal("-pc_hypre_ams_relax_weight", "Relaxation weight for AMS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3)); PetscCall(PetscOptionsReal("-pc_hypre_ams_omega", "SSOR coefficient for AMS smoother", "None", jac->as_omega, &jac->as_omega, &flag4)); if (flag || flag2 || flag3 || flag4) PetscCallHYPRE(HYPRE_AMSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega)); 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)); n = 5; PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options", "AMG options for vector Poisson", "None", jac->as_amg_alpha_opts, &n, &flag2)); if (flag || flag2) { PetscCallHYPRE(HYPRE_AMSSetAlphaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */ jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */ } 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)); n = 5; PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options", "AMG options for scalar Poisson solver", "None", jac->as_amg_beta_opts, &n, &flag2)); if (flag || flag2) { PetscCallHYPRE(HYPRE_AMSSetBetaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */ jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */ } 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)); if (flag) { /* override HYPRE's default only if the options is used */ PetscCallHYPRE(HYPRE_AMSSetProjectionFrequency(jac->hsolver, (HYPRE_Int)jac->ams_proj_freq)); } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_AMS(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool isascii; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE AMS preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter)); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type)); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type)); PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega)); if (jac->alpha_Poisson) { PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (passed in by user)\n")); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (computed) \n")); } PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_alpha_theta)); if (!jac->ams_beta_is_zero) { if (jac->beta_Poisson) { PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (passed in by user)\n")); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (computed) \n")); } PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4])); PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_beta_theta)); 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)); } else { PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver not used (zero-conductivity everywhere) \n")); } } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PC pc, PetscOptionItems PetscOptionsObject) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscInt n; PetscBool flag, flag2, flag3, flag4; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ADS Options"); PetscCall(PetscOptionsInt("-pc_hypre_ads_print_level", "Debugging output level for ADS", "None", jac->as_print, &jac->as_print, &flag)); if (flag) PetscCallHYPRE(HYPRE_ADSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print)); 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)); if (flag) PetscCallHYPRE(HYPRE_ADSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter)); PetscCall(PetscOptionsInt("-pc_hypre_ads_cycle_type", "Cycle type for ADS multigrid", "None", jac->ads_cycle_type, &jac->ads_cycle_type, &flag)); if (flag) PetscCallHYPRE(HYPRE_ADSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ads_cycle_type)); PetscCall(PetscOptionsReal("-pc_hypre_ads_tol", "Error tolerance for ADS multigrid", "None", jac->as_tol, &jac->as_tol, &flag)); if (flag) PetscCallHYPRE(HYPRE_ADSSetTol(jac->hsolver, jac->as_tol)); PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_type", "Relaxation type for ADS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag)); PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_times", "Number of relaxation steps for ADS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2)); PetscCall(PetscOptionsReal("-pc_hypre_ads_relax_weight", "Relaxation weight for ADS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3)); PetscCall(PetscOptionsReal("-pc_hypre_ads_omega", "SSOR coefficient for ADS smoother", "None", jac->as_omega, &jac->as_omega, &flag4)); if (flag || flag2 || flag3 || flag4) PetscCallHYPRE(HYPRE_ADSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega)); 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)); n = 5; PetscCall(PetscOptionsIntArray("-pc_hypre_ads_ams_options", "AMG options for AMS solver inside ADS", "None", jac->as_amg_alpha_opts, &n, &flag2)); 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)); if (flag || flag2 || flag3) { PetscCallHYPRE(HYPRE_ADSSetAMSOptions(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type, /* AMS cycle type */ (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */ jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */ } 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)); n = 5; PetscCall(PetscOptionsIntArray("-pc_hypre_ads_amg_options", "AMG options for vector AMG solver inside ADS", "None", jac->as_amg_beta_opts, &n, &flag2)); if (flag || flag2) { PetscCallHYPRE(HYPRE_ADSSetAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */ jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */ } PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_HYPRE_ADS(PC pc, PetscViewer viewer) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool isascii; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ADS preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter)); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type)); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type)); PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight)); PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega)); PetscCall(PetscViewerASCIIPrintf(viewer, " AMS solver using boomerAMG\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type)); PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0])); PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1])); PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2])); PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3])); PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4])); PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_alpha_theta)); PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver using boomerAMG\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0])); PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1])); PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2])); PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3])); PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4])); PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_beta_theta)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool ishypre; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)G)); PetscCall(MatDestroy(&jac->G)); jac->G = G; } else { PetscCall(MatDestroy(&jac->G)); PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetDiscreteGradient - Set the discrete gradient matrix for `PCHYPRE` type of AMS or ADS Collective Input Parameters: + pc - the preconditioning context - G - the discrete gradient Level: intermediate Notes: `G` should have as many rows as the number of edges and as many columns as the number of vertices in the mesh 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 Developer Note: This automatically converts the matrix to `MATHYPRE` if it is not already of that type .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteCurl()` @*/ PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscValidHeaderSpecific(G, MAT_CLASSID, 2); PetscCheckSameComm(pc, 1, G, 2); PetscTryMethod(pc, "PCHYPRESetDiscreteGradient_C", (PC, Mat), (pc, G)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool ishypre; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)C)); PetscCall(MatDestroy(&jac->C)); jac->C = C; } else { PetscCall(MatDestroy(&jac->C)); PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C)); } PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetDiscreteCurl - Set the discrete curl matrix for `PCHYPRE` type of ADS Collective Input Parameters: + pc - the preconditioning context - C - the discrete curl Level: intermediate Notes: `C` should have as many rows as the number of faces and as many columns as the number of edges in the mesh 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 Developer Notes: This automatically converts the matrix to `MATHYPRE` if it is not already of that type If this is only for `PCHYPRE` type of ADS it should be called `PCHYPREADSSetDiscreteCurl()` .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()` @*/ PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscValidHeaderSpecific(C, MAT_CLASSID, 2); PetscCheckSameComm(pc, 1, C, 2); PetscTryMethod(pc, "PCHYPRESetDiscreteCurl_C", (PC, Mat), (pc, C)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool ishypre; PetscInt i; PetscFunctionBegin; PetscCall(MatDestroy(&jac->RT_PiFull)); PetscCall(MatDestroy(&jac->ND_PiFull)); for (i = 0; i < 3; ++i) { PetscCall(MatDestroy(&jac->RT_Pi[i])); PetscCall(MatDestroy(&jac->ND_Pi[i])); } jac->dim = dim; if (RT_PiFull) { PetscCall(PetscObjectTypeCompare((PetscObject)RT_PiFull, MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)RT_PiFull)); jac->RT_PiFull = RT_PiFull; } else { PetscCall(MatConvert(RT_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_PiFull)); } } if (RT_Pi) { for (i = 0; i < dim; ++i) { if (RT_Pi[i]) { PetscCall(PetscObjectTypeCompare((PetscObject)RT_Pi[i], MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)RT_Pi[i])); jac->RT_Pi[i] = RT_Pi[i]; } else { PetscCall(MatConvert(RT_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_Pi[i])); } } } } if (ND_PiFull) { PetscCall(PetscObjectTypeCompare((PetscObject)ND_PiFull, MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)ND_PiFull)); jac->ND_PiFull = ND_PiFull; } else { PetscCall(MatConvert(ND_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_PiFull)); } } if (ND_Pi) { for (i = 0; i < dim; ++i) { if (ND_Pi[i]) { PetscCall(PetscObjectTypeCompare((PetscObject)ND_Pi[i], MATHYPRE, &ishypre)); if (ishypre) { PetscCall(PetscObjectReference((PetscObject)ND_Pi[i])); jac->ND_Pi[i] = ND_Pi[i]; } else { PetscCall(MatConvert(ND_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_Pi[i])); } } } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetInterpolations - Set the interpolation matrices for `PCHYPRE` type of AMS or ADS Collective Input Parameters: + pc - the preconditioning context . dim - the dimension of the problem, only used in AMS . RT_PiFull - Raviart-Thomas interpolation matrix . RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix . ND_PiFull - Nedelec interpolation matrix - ND_Pi - x/y/z component of Nedelec interpolation matrix Level: intermediate Notes: For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to `NULL`. For ADS, both type of interpolation matrices are needed. Developer Note: This automatically converts the matrix to `MATHYPRE` if it is not already of that type .seealso: [](ch_ksp), `PCHYPRE` @*/ PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) { PetscInt i; PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); if (RT_PiFull) { PetscValidHeaderSpecific(RT_PiFull, MAT_CLASSID, 3); PetscCheckSameComm(pc, 1, RT_PiFull, 3); } if (RT_Pi) { PetscAssertPointer(RT_Pi, 4); for (i = 0; i < dim; ++i) { if (RT_Pi[i]) { PetscValidHeaderSpecific(RT_Pi[i], MAT_CLASSID, 4); PetscCheckSameComm(pc, 1, RT_Pi[i], 4); } } } if (ND_PiFull) { PetscValidHeaderSpecific(ND_PiFull, MAT_CLASSID, 5); PetscCheckSameComm(pc, 1, ND_PiFull, 5); } if (ND_Pi) { PetscAssertPointer(ND_Pi, 6); for (i = 0; i < dim; ++i) { if (ND_Pi[i]) { PetscValidHeaderSpecific(ND_Pi[i], MAT_CLASSID, 6); PetscCheckSameComm(pc, 1, ND_Pi[i], 6); } } } PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool ishypre; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre)); if (ishypre) { if (isalpha) { PetscCall(PetscObjectReference((PetscObject)A)); PetscCall(MatDestroy(&jac->alpha_Poisson)); jac->alpha_Poisson = A; } else { if (A) { PetscCall(PetscObjectReference((PetscObject)A)); } else { jac->ams_beta_is_zero = PETSC_TRUE; } PetscCall(MatDestroy(&jac->beta_Poisson)); jac->beta_Poisson = A; } } else { if (isalpha) { PetscCall(MatDestroy(&jac->alpha_Poisson)); PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson)); } else { if (A) { PetscCall(MatDestroy(&jac->beta_Poisson)); PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson)); } else { PetscCall(MatDestroy(&jac->beta_Poisson)); jac->ams_beta_is_zero = PETSC_TRUE; } } } PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetAlphaPoissonMatrix - Set the vector Poisson matrix for `PCHYPRE` of type AMS Collective Input Parameters: + pc - the preconditioning context - A - the matrix Level: intermediate Note: `A` should be obtained by discretizing the vector valued Poisson problem with linear finite elements Developer Notes: This automatically converts the matrix to `MATHYPRE` if it is not already of that type If this is only for `PCHYPRE` type of AMS it should be called `PCHYPREAMSSetAlphaPoissonMatrix()` .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()` @*/ PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscValidHeaderSpecific(A, MAT_CLASSID, 2); PetscCheckSameComm(pc, 1, A, 2); PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetBetaPoissonMatrix - Set the Poisson matrix for `PCHYPRE` of type AMS Collective Input Parameters: + pc - the preconditioning context - A - the matrix, or `NULL` to turn it off Level: intermediate Note: `A` should be obtained by discretizing the Poisson problem with linear finite elements. Developer Notes: This automatically converts the matrix to `MATHYPRE` if it is not already of that type If this is only for `PCHYPRE` type of AMS it should be called `PCHYPREAMSPCHYPRESetBetaPoissonMatrix()` .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` @*/ PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); if (A) { PetscValidHeaderSpecific(A, MAT_CLASSID, 2); PetscCheckSameComm(pc, 1, A, 2); } PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc, Vec ozz, Vec zoz, Vec zzo) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; /* throw away any vector if already set */ PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2])); PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0])); PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0])); PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1])); PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1])); jac->dim = 2; if (zzo) { PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2])); PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2])); jac->dim++; } PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in the edge element basis for `PCHYPRE` of type AMS Collective Input Parameters: + pc - the preconditioning context . ozz - vector representing (1,0,0) (or (1,0) in 2D) . zoz - vector representing (0,1,0) (or (0,1) in 2D) - zzo - vector representing (0,0,1) (use NULL in 2D) Level: intermediate Developer Note: If this is only for `PCHYPRE` type of AMS it should be called `PCHYPREAMSSetEdgeConstantVectors()` .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` @*/ PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscValidHeaderSpecific(ozz, VEC_CLASSID, 2); PetscValidHeaderSpecific(zoz, VEC_CLASSID, 3); if (zzo) PetscValidHeaderSpecific(zzo, VEC_CLASSID, 4); PetscCheckSameComm(pc, 1, ozz, 2); PetscCheckSameComm(pc, 1, zoz, 3); if (zzo) PetscCheckSameComm(pc, 1, zzo, 4); PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPREAMSSetInteriorNodes_HYPRE(PC pc, Vec interior) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior)); PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior)); PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior)); jac->ams_beta_is_zero_part = PETSC_TRUE; PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPREAMSSetInteriorNodes - Set the list of interior nodes to a zero-conductivity region for `PCHYPRE` of type AMS Collective Input Parameters: + pc - the preconditioning context - interior - vector. node is interior if its entry in the array is 1.0. Level: intermediate Note: This calls `HYPRE_AMSSetInteriorNodes()` .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()` @*/ PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscValidHeaderSpecific(interior, VEC_CLASSID, 2); PetscCheckSameComm(pc, 1, interior, 2); PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; Vec tv; PetscInt i; PetscFunctionBegin; /* throw away any coordinate vector if already set */ PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1])); PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2])); jac->dim = dim; /* compute IJ vector for coordinates */ PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &tv)); PetscCall(VecSetType(tv, VECSTANDARD)); PetscCall(VecSetSizes(tv, nloc, PETSC_DECIDE)); for (i = 0; i < dim; i++) { PetscScalar *array; PetscInt j; PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i])); PetscCall(VecGetArrayWrite(tv, &array)); for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i]; PetscCall(VecRestoreArrayWrite(tv, &array)); PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i])); } PetscCall(VecDestroy(&tv)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPREGetType_HYPRE(PC pc, const char *name[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; *name = jac->hypre_type; PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCHYPRESetType_HYPRE(PC pc, const char name[]) { PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscBool flag; PetscFunctionBegin; if (jac->hypre_type) { PetscCall(PetscStrcmp(jac->hypre_type, name, &flag)); if (flag) PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PCReset_HYPRE(pc)); PetscCall(PetscFree(jac->hypre_type)); PetscCall(PetscStrallocpy(name, &jac->hypre_type)); jac->maxiter = PETSC_DEFAULT; jac->tol = PETSC_DEFAULT; jac->printstatistics = PetscLogPrintInfo; PetscCall(PetscStrcmp("ilu", jac->hypre_type, &flag)); if (flag) { PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); PetscCallHYPRE(HYPRE_ILUCreate(&jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ILU; pc->ops->view = PCView_HYPRE_ILU; jac->destroy = HYPRE_ILUDestroy; jac->setup = HYPRE_ILUSetup; jac->solve = HYPRE_ILUSolve; jac->factorrowsize = PETSC_DEFAULT; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag)); if (flag) { PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); PetscCallHYPRE(HYPRE_ParCSRPilutCreate(jac->comm_hypre, &jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; pc->ops->view = PCView_HYPRE_Pilut; jac->destroy = HYPRE_ParCSRPilutDestroy; jac->setup = HYPRE_ParCSRPilutSetup; jac->solve = HYPRE_ParCSRPilutSolve; jac->factorrowsize = PETSC_DEFAULT; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag)); if (flag) { #if defined(PETSC_USE_64BIT_INDICES) SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices"); #endif PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); PetscCallHYPRE(HYPRE_EuclidCreate(jac->comm_hypre, &jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid; pc->ops->view = PCView_HYPRE_Euclid; jac->destroy = HYPRE_EuclidDestroy; jac->setup = HYPRE_EuclidSetup; jac->solve = HYPRE_EuclidSolve; jac->factorrowsize = PETSC_DEFAULT; jac->eu_level = PETSC_DEFAULT; /* default */ PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag)); if (flag) { PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre)); PetscCallHYPRE(HYPRE_ParaSailsCreate(jac->comm_hypre, &jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; pc->ops->view = PCView_HYPRE_ParaSails; jac->destroy = HYPRE_ParaSailsDestroy; jac->setup = HYPRE_ParaSailsSetup; jac->solve = HYPRE_ParaSailsSolve; /* initialize */ jac->nlevels = 1; jac->threshold = .1; jac->filter = .1; jac->loadbal = 0; if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE; else jac->logging = (int)PETSC_FALSE; jac->ruse = (int)PETSC_FALSE; jac->symt = 0; PetscCallHYPRE(HYPRE_ParaSailsSetParams(jac->hsolver, jac->threshold, (HYPRE_Int)jac->nlevels)); PetscCallHYPRE(HYPRE_ParaSailsSetFilter(jac->hsolver, jac->filter)); PetscCallHYPRE(HYPRE_ParaSailsSetLoadbal(jac->hsolver, (HYPRE_Int)jac->loadbal)); PetscCallHYPRE(HYPRE_ParaSailsSetLogging(jac->hsolver, (HYPRE_Int)jac->logging)); PetscCallHYPRE(HYPRE_ParaSailsSetReuse(jac->hsolver, (HYPRE_Int)jac->ruse)); PetscCallHYPRE(HYPRE_ParaSailsSetSym(jac->hsolver, (HYPRE_Int)jac->symt)); PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag)); if (flag) { PetscCallHYPRE(HYPRE_BoomerAMGCreate(&jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; pc->ops->view = PCView_HYPRE_BoomerAMG; pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; pc->ops->matapply = PCMatApply_HYPRE_BoomerAMG; PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", PCHYPREGetCFMarkers_BoomerAMG)); jac->destroy = HYPRE_BoomerAMGDestroy; jac->setup = HYPRE_BoomerAMGSetup; jac->solve = HYPRE_BoomerAMGSolve; jac->applyrichardson = PETSC_FALSE; /* these defaults match the hypre defaults */ jac->cycletype = 1; jac->maxlevels = 25; jac->maxiter = 1; jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ jac->truncfactor = 0.0; jac->strongthreshold = .25; jac->maxrowsum = .9; jac->measuretype = 0; jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; jac->smoothtype = -1; /* Not set by default */ jac->smoothnumlevels = 25; jac->eu_level = 0; jac->eu_droptolerance = 0; jac->eu_bj = 0; jac->relaxweight = 1.0; jac->outerrelaxweight = 1.0; jac->Rtype = 0; jac->Rstrongthreshold = 0.25; jac->Rfilterthreshold = 0.0; jac->Adroptype = -1; jac->Adroptol = 0.0; jac->agg_nl = 0; jac->pmax = 0; jac->truncfactor = 0.0; jac->agg_num_paths = 1; jac->maxc = 9; jac->minc = 1; jac->nodal_coarsening = 0; jac->nodal_coarsening_diag = 0; jac->vec_interp_variant = 0; jac->vec_interp_qmax = 0; jac->vec_interp_smooth = PETSC_FALSE; jac->interp_refine = 0; jac->nodal_relax = PETSC_FALSE; jac->nodal_relax_levels = 1; jac->rap2 = 0; PetscObjectParameterSetDefault(jac, relaxtype[2], 9); /* G.E. */ /* Initialize the following parameters with invalid value so we can recognize user input that sets the parameter. If there is no user input they are overwritten in PCSetUp_HYPRE() depending on if the matrix is on the CPU or the GPU */ PetscObjectParameterSetDefault(jac, relaxorder, PETSC_DECIDE); PetscObjectParameterSetDefault(jac, coarsentype, PETSC_DECIDE); PetscObjectParameterSetDefault(jac, interptype, PETSC_DECIDE); PetscObjectParameterSetDefault(jac, relaxtype[0], PETSC_DECIDE); PetscObjectParameterSetDefault(jac, relaxtype[1], PETSC_DECIDE); #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0) PetscObjectParameterSetDefault(jac, spgemm_type, "not yet set"); #endif #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0) PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_UNKNOWN); PetscObjectParameterSetDefault(jac, mod_rap2, PETSC_DECIDE); #endif PetscObjectParameterSetDefault(jac, agg_interptype, PETSC_DECIDE); PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag)); if (flag) { PetscCallHYPRE(HYPRE_AMSCreate(&jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; pc->ops->view = PCView_HYPRE_AMS; jac->destroy = HYPRE_AMSDestroy; jac->setup = HYPRE_AMSSetup; jac->solve = HYPRE_AMSSolve; jac->coords[0] = NULL; jac->coords[1] = NULL; jac->coords[2] = NULL; jac->interior = NULL; /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ jac->as_print = 0; jac->as_max_iter = 1; /* used as a preconditioner */ jac->as_tol = 0.; /* used as a preconditioner */ jac->ams_cycle_type = 13; /* Smoothing options */ jac->as_relax_type = 2; jac->as_relax_times = 1; jac->as_relax_weight = 1.0; jac->as_omega = 1.0; /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ jac->as_amg_alpha_opts[0] = 10; jac->as_amg_alpha_opts[1] = 1; jac->as_amg_alpha_opts[2] = 6; jac->as_amg_alpha_opts[3] = 6; jac->as_amg_alpha_opts[4] = 4; jac->as_amg_alpha_theta = 0.25; /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ jac->as_amg_beta_opts[0] = 10; jac->as_amg_beta_opts[1] = 1; jac->as_amg_beta_opts[2] = 6; jac->as_amg_beta_opts[3] = 6; jac->as_amg_beta_opts[4] = 4; jac->as_amg_beta_theta = 0.25; PetscCallHYPRE(HYPRE_AMSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print)); PetscCallHYPRE(HYPRE_AMSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter)); PetscCallHYPRE(HYPRE_AMSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type)); PetscCallHYPRE(HYPRE_AMSSetTol(jac->hsolver, jac->as_tol)); PetscCallHYPRE(HYPRE_AMSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega)); PetscCallHYPRE(HYPRE_AMSSetAlphaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */ jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */ PetscCallHYPRE(HYPRE_AMSSetBetaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */ jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */ /* Zero conductivity */ jac->ams_beta_is_zero = PETSC_FALSE; jac->ams_beta_is_zero_part = PETSC_FALSE; PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag)); if (flag) { PetscCallHYPRE(HYPRE_ADSCreate(&jac->hsolver)); pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS; pc->ops->view = PCView_HYPRE_ADS; jac->destroy = HYPRE_ADSDestroy; jac->setup = HYPRE_ADSSetup; jac->solve = HYPRE_ADSSolve; jac->coords[0] = NULL; jac->coords[1] = NULL; jac->coords[2] = NULL; /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */ jac->as_print = 0; jac->as_max_iter = 1; /* used as a preconditioner */ jac->as_tol = 0.; /* used as a preconditioner */ jac->ads_cycle_type = 13; /* Smoothing options */ jac->as_relax_type = 2; jac->as_relax_times = 1; jac->as_relax_weight = 1.0; jac->as_omega = 1.0; /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */ jac->ams_cycle_type = 14; jac->as_amg_alpha_opts[0] = 10; jac->as_amg_alpha_opts[1] = 1; jac->as_amg_alpha_opts[2] = 6; jac->as_amg_alpha_opts[3] = 6; jac->as_amg_alpha_opts[4] = 4; jac->as_amg_alpha_theta = 0.25; /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ jac->as_amg_beta_opts[0] = 10; jac->as_amg_beta_opts[1] = 1; jac->as_amg_beta_opts[2] = 6; jac->as_amg_beta_opts[3] = 6; jac->as_amg_beta_opts[4] = 4; jac->as_amg_beta_theta = 0.25; PetscCallHYPRE(HYPRE_ADSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print)); PetscCallHYPRE(HYPRE_ADSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter)); PetscCallHYPRE(HYPRE_ADSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type)); PetscCallHYPRE(HYPRE_ADSSetTol(jac->hsolver, jac->as_tol)); PetscCallHYPRE(HYPRE_ADSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega)); PetscCallHYPRE(HYPRE_ADSSetAMSOptions(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type, /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */ jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */ PetscCallHYPRE(HYPRE_ADSSetAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */ (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */ (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */ jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */ (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */ PetscFunctionReturn(PETSC_SUCCESS); } PetscCall(PetscFree(jac->hypre_type)); jac->hypre_type = NULL; SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE preconditioner %s; Choices are euclid, ilu, pilut, parasails, boomeramg, ams, ads", name); } /* It only gets here if the HYPRE type has not been set before the call to ...SetFromOptions() which actually is most of the time */ static PetscErrorCode PCSetFromOptions_HYPRE(PC pc, PetscOptionItems PetscOptionsObject) { PetscInt indx; const char *type[] = {"ilu", "euclid", "pilut", "parasails", "boomeramg", "ams", "ads"}; PetscBool flg; PC_HYPRE *jac = (PC_HYPRE *)pc->data; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options"); PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg)); if (flg) PetscCall(PCHYPRESetType_HYPRE(pc, type[indx])); /* Set the type if it was never set. */ if (!jac->hypre_type) PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg")); PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPRESetType - Sets which hypre preconditioner you wish to use Input Parameters: + pc - the preconditioner context - name - either euclid, ilu, pilut, parasails, boomeramg, ams, or ads Options Database Key: . pc_hypre_type - One of euclid, ilu, pilut, parasails, boomeramg, ams, or ads Level: intermediate .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRE` @*/ PetscErrorCode PCHYPRESetType(PC pc, const char name[]) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscAssertPointer(name, 2); PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name)); PetscFunctionReturn(PETSC_SUCCESS); } /*@C PCHYPREGetCFMarkers - Gets CF marker arrays for all levels (except the finest level) Logically Collective Input Parameter: . pc - the preconditioner context Output Parameters: + n_per_level - the number of nodes per level (size of `num_levels`) - CFMarkers - the Coarse/Fine Boolean arrays (size of `num_levels` - 1) Level: advanced Note: 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. .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()` @*/ PetscErrorCode PCHYPREGetCFMarkers(PC pc, PetscInt *n_per_level[], PetscBT *CFMarkers[]) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscAssertPointer(n_per_level, 2); PetscAssertPointer(CFMarkers, 3); PetscUseMethod(pc, "PCHYPREGetCFMarkers_C", (PC, PetscInt *[], PetscBT *[]), (pc, n_per_level, CFMarkers)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCHYPREGetType - Gets which hypre preconditioner you are using Input Parameter: . pc - the preconditioner context Output Parameter: . name - either euclid, ilu, pilut, parasails, boomeramg, ams, or ads Level: intermediate .seealso: [](ch_ksp), `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, `PCHYPRE` @*/ PetscErrorCode PCHYPREGetType(PC pc, const char *name[]) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscAssertPointer(name, 2); PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCMGGalerkinSetMatProductAlgorithm - Set type of sparse matrix-matrix product for hypre's BoomerAMG to use on GPUs Logically Collective Input Parameters: + pc - the hypre context - name - one of 'cusparse', 'hypre' Options Database Key: . -pc_mg_galerkin_mat_product_algorithm - Type of sparse matrix-matrix product to use in hypre Level: intermediate Developer Note: How the name starts with `PCMG`, should it not be `PCHYPREBoomerAMG`? .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinGetMatProductAlgorithm()` @*/ PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[]) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name)); PetscFunctionReturn(PETSC_SUCCESS); } /*@ PCMGGalerkinGetMatProductAlgorithm - Get type of sparse matrix-matrix product for hypre's BoomerAMG to use on GPUs Not Collective Input Parameter: . pc - the multigrid context Output Parameter: . name - one of 'cusparse', 'hypre' Level: intermediate .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinSetMatProductAlgorithm()` @*/ PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[]) { PetscFunctionBegin; PetscValidHeaderSpecific(pc, PC_CLASSID, 1); PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name)); PetscFunctionReturn(PETSC_SUCCESS); } /*MC PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre as PETSc `PC` Options Database Keys: + -pc_hypre_type - One of `euclid`, `ilu`, `pilut`, `parasails`, `boomeramg`, `ams`, or `ads` . -pc_hypre_boomeramg_nodal_coarsen - where `n` is from 1 to 6 (see `HYPRE_BoomerAMGSetNodal()`) . -pc_hypre_boomeramg_vec_interp_variant - where `v` is from 1 to 3 (see `HYPRE_BoomerAMGSetInterpVecVariant()`) - Many others - run with `-pc_type hypre` `-pc_hypre_type XXX` `-help` to see options for the XXX preconditioner Level: intermediate Notes: Apart from `-pc_hypre_type` (for which there is `PCHYPRESetType()`), the many hypre options can ONLY be set via the options database (e.g. the command line or with `PetscOptionsSetValue()`, there are no functions to set them) The options `-pc_hypre_boomeramg_max_iter` and `-pc_hypre_boomeramg_tol` refer to the number of iterations (V-cycles) and tolerance that boomerAMG does EACH time it is called. So for example, if `-pc_hypre_boomeramg_max_iter` is set to 2 then 2-V-cycles are being used to define the preconditioner (`-pc_hypre_boomeramg_tol` should be set to 0.0 - the default - to strictly use a fixed number of iterations per hypre call). `-ksp_max_it` and `-ksp_rtol` STILL determine the total number of iterations and tolerance for the Krylov solver. For example, if `-pc_hypre_boomeramg_max_iter` is 2 and `-ksp_max_it` is 10 then AT MOST twenty V-cycles of BoomerAMG will be used. Note that the option `-pc_hypre_boomeramg_relax_type_all` defaults to symmetric relaxation (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. Otherwise, you may want to use `-pc_hypre_boomeramg_relax_type_all SOR/Jacobi`. If you provide a near null space to your matrix with `MatSetNearNullSpace()` it is ignored by hypre's BoomerAMG UNLESS you also use the following two options: `-pc_hypre_boomeramg_nodal_coarsen -pc_hypre_boomeramg_vec_interp_variant ` See `PCPFMG`, `PCSMG`, and `PCSYSPFMG` for access to hypre's other (nonalgebraic) multigrid solvers For `PCHYPRE` type of `ams` or `ads` auxiliary data must be provided to the preconditioner with `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`, `PCHYPREAMSSetInteriorNodes()` 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 since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON` (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid. PETSc provides its own geometric and algebraic multigrid solvers `PCMG` and `PCGAMG`, also see `PCHMG` which is useful for certain multicomponent problems. hypre supports performance logging via the `Caliper` library. With `--download-hypre --download-caliper`, hypre will be automatically configured with the support. Enabling Caliper logging requires setting the `CALI_CONFIG` environment variable before running your hypre code. For example, .vb export CALI_CONFIG=runtime-report,max_column_width=200,calc.inclusive,mpi-report,output=stdout .ve Then run a hypre code, and you will see profiling results on stdout. See https://software.llnl.gov/Caliper/#guides for more options. GPU Notes: To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda Then pass `VECCUDA` vectors and `MATAIJCUSPARSE` matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers. To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip Then pass `VECHIP` vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers. .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`, PCHYPREAMSSetInteriorNodes() M*/ PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) { PC_HYPRE *jac; PetscFunctionBegin; PetscCall(PetscNew(&jac)); pc->data = jac; pc->ops->reset = PCReset_HYPRE; pc->ops->destroy = PCDestroy_HYPRE; pc->ops->setfromoptions = PCSetFromOptions_HYPRE; pc->ops->setup = PCSetUp_HYPRE; pc->ops->apply = PCApply_HYPRE; jac->hypre_type = NULL; jac->comm_hypre = MPI_COMM_NULL; PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG)); PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG)); #if defined(PETSC_HAVE_HYPRE_DEVICE) #if defined(HYPRE_USING_HIP) PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP)); #endif #if defined(HYPRE_USING_CUDA) PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA)); #endif #endif PetscCall(PetscHYPREInitialize()); PetscFunctionReturn(PETSC_SUCCESS); } typedef struct { MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ HYPRE_StructSolver hsolver; /* keep copy of PFMG options used so may view them */ PetscInt its; PetscReal tol; PetscInt relax_type; PetscInt rap_type; PetscInt num_pre_relax, num_post_relax; PetscInt max_levels; PetscInt skip_relax; PetscBool print_statistics; } PC_PFMG; static PetscErrorCode PCDestroy_PFMG(PC pc) { PC_PFMG *ex = (PC_PFMG *)pc->data; PetscFunctionBegin; if (ex->hsolver) PetscCallHYPRE(HYPRE_StructPFMGDestroy(ex->hsolver)); PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscFree(pc->data)); PetscFunctionReturn(PETSC_SUCCESS); } static const char *PFMGRelaxType[] = {"Jacobi", "Weighted-Jacobi", "symmetric-Red/Black-Gauss-Seidel", "Red/Black-Gauss-Seidel"}; static const char *PFMGRAPType[] = {"Galerkin", "non-Galerkin"}; static PetscErrorCode PCView_PFMG(PC pc, PetscViewer viewer) { PetscBool isascii; PC_PFMG *ex = (PC_PFMG *)pc->data; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE PFMG preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type])); PetscCall(PetscViewerASCIIPrintf(viewer, " RAP type %s\n", PFMGRAPType[ex->rap_type])); PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); PetscCall(PetscViewerASCIIPrintf(viewer, " max levels %" PetscInt_FMT "\n", ex->max_levels)); PetscCall(PetscViewerASCIIPrintf(viewer, " skip relax %" PetscInt_FMT "\n", ex->skip_relax)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems PetscOptionsObject) { PC_PFMG *ex = (PC_PFMG *)pc->data; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options"); PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL)); PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL)); PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its)); 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)); PetscCallHYPRE(HYPRE_StructPFMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax)); 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)); PetscCallHYPRE(HYPRE_StructPFMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax)); PetscCall(PetscOptionsInt("-pc_pfmg_max_levels", "Max Levels for MG hierarchy", "HYPRE_StructPFMGSetMaxLevels", ex->max_levels, &ex->max_levels, NULL)); PetscCallHYPRE(HYPRE_StructPFMGSetMaxLevels(ex->hsolver, (HYPRE_Int)ex->max_levels)); PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL)); PetscCallHYPRE(HYPRE_StructPFMGSetTol(ex->hsolver, ex->tol)); 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)); PetscCallHYPRE(HYPRE_StructPFMGSetRelaxType(ex->hsolver, (HYPRE_Int)ex->relax_type)); PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL)); PetscCallHYPRE(HYPRE_StructPFMGSetRAPType(ex->hsolver, (HYPRE_Int)ex->rap_type)); 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)); PetscCallHYPRE(HYPRE_StructPFMGSetSkipRelax(ex->hsolver, (HYPRE_Int)ex->skip_relax)); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y) { PC_PFMG *ex = (PC_PFMG *)pc->data; PetscScalar *yy; const PetscScalar *xx; PetscInt ilower[3], iupper[3]; HYPRE_Int hlower[3], hupper[3]; Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ iupper[0] += ilower[0] - 1; iupper[1] += ilower[1] - 1; iupper[2] += ilower[2] - 1; hlower[0] = (HYPRE_Int)ilower[0]; hlower[1] = (HYPRE_Int)ilower[1]; hlower[2] = (HYPRE_Int)ilower[2]; hupper[0] = (HYPRE_Int)iupper[0]; hupper[1] = (HYPRE_Int)iupper[1]; hupper[2] = (HYPRE_Int)iupper[2]; /* copy x values over to hypre */ PetscCallHYPRE(HYPRE_StructVectorSetConstantValues(mx->hb, 0.0)); PetscCall(VecGetArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_StructVectorSetBoxValues(mx->hb, hlower, hupper, (HYPRE_Complex *)xx)); PetscCall(VecRestoreArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_StructVectorAssemble(mx->hb)); PetscCallHYPRE(HYPRE_StructPFMGSolve(ex->hsolver, mx->hmat, mx->hb, mx->hx)); /* copy solution values back to PETSc */ PetscCall(VecGetArray(y, &yy)); PetscCallHYPRE(HYPRE_StructVectorGetBoxValues(mx->hx, hlower, hupper, (HYPRE_Complex *)yy)); PetscCall(VecRestoreArray(y, &yy)); PetscFunctionReturn(PETSC_SUCCESS); } 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) { PC_PFMG *jac = (PC_PFMG *)pc->data; HYPRE_Int oits; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->its))); PetscCallHYPRE(HYPRE_StructPFMGSetTol(jac->hsolver, rtol)); PetscCall(PCApply_PFMG(pc, b, y)); PetscCallHYPRE(HYPRE_StructPFMGGetNumIterations(jac->hsolver, &oits)); *outits = oits; if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; else *reason = PCRICHARDSON_CONVERGED_RTOL; PetscCallHYPRE(HYPRE_StructPFMGSetTol(jac->hsolver, jac->tol)); PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->its)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetUp_PFMG(PC pc) { PC_PFMG *ex = (PC_PFMG *)pc->data; Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; PetscBool flg; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg)); PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner"); /* create the hypre solver object and set its information */ if (ex->hsolver) PetscCallHYPRE(HYPRE_StructPFMGDestroy(ex->hsolver)); PetscCallHYPRE(HYPRE_StructPFMGCreate(ex->hcomm, &ex->hsolver)); // Print Hypre statistics about the solve process if (ex->print_statistics) PetscCallHYPRE(HYPRE_StructPFMGSetPrintLevel(ex->hsolver, 3)); // The hypre options must be repeated here because the StructPFMG was destroyed and recreated PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its)); PetscCallHYPRE(HYPRE_StructPFMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax)); PetscCallHYPRE(HYPRE_StructPFMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax)); PetscCallHYPRE(HYPRE_StructPFMGSetMaxLevels(ex->hsolver, (HYPRE_Int)ex->max_levels)); PetscCallHYPRE(HYPRE_StructPFMGSetTol(ex->hsolver, ex->tol)); PetscCallHYPRE(HYPRE_StructPFMGSetRelaxType(ex->hsolver, (HYPRE_Int)ex->relax_type)); PetscCallHYPRE(HYPRE_StructPFMGSetRAPType(ex->hsolver, (HYPRE_Int)ex->rap_type)); PetscCallHYPRE(HYPRE_StructPFMGSetup(ex->hsolver, mx->hmat, mx->hb, mx->hx)); PetscCallHYPRE(HYPRE_StructPFMGSetZeroGuess(ex->hsolver)); PetscFunctionReturn(PETSC_SUCCESS); } /*MC PCPFMG - the hypre PFMG multigrid solver Options Database Keys: + -pc_pfmg_its - number of iterations of PFMG to use as preconditioner . -pc_pfmg_num_pre_relax - number of smoothing steps before coarse grid solve . -pc_pfmg_num_post_relax - number of smoothing steps after coarse grid solve . -pc_pfmg_tol - tolerance of PFMG . -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 . -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin - -pc_pfmg_skip_relax - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations when the underlying problem is isotropic, one of 0,1 Level: advanced Notes: This is for CELL-centered descretizations See `PCSYSPFMG` for a version suitable for systems of PDEs, and `PCSMG` See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver This must be used with the `MATHYPRESTRUCT` matrix type. This provides only some of the functionality of PFMG, it supports only one block per process defined by a PETSc `DMDA`. .seealso: [](ch_ksp), `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG` M*/ PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) { PC_PFMG *ex; PetscFunctionBegin; PetscCall(PetscNew(&ex)); pc->data = ex; ex->its = 1; ex->tol = 1.e-8; ex->relax_type = 1; ex->rap_type = 0; ex->num_pre_relax = 1; ex->num_post_relax = 1; ex->max_levels = 0; ex->skip_relax = 0; ex->print_statistics = PETSC_FALSE; pc->ops->setfromoptions = PCSetFromOptions_PFMG; pc->ops->view = PCView_PFMG; pc->ops->destroy = PCDestroy_PFMG; pc->ops->apply = PCApply_PFMG; pc->ops->applyrichardson = PCApplyRichardson_PFMG; pc->ops->setup = PCSetUp_PFMG; PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscHYPREInitialize()); PetscCallHYPRE(HYPRE_StructPFMGCreate(ex->hcomm, &ex->hsolver)); PetscFunctionReturn(PETSC_SUCCESS); } /* we know we are working with a HYPRE_SStructMatrix */ typedef struct { MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ HYPRE_SStructSolver ss_solver; /* keep copy of SYSPFMG options used so may view them */ PetscInt its; PetscReal tol; PetscInt relax_type; PetscInt num_pre_relax, num_post_relax; } PC_SysPFMG; static PetscErrorCode PCDestroy_SysPFMG(PC pc) { PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; PetscFunctionBegin; if (ex->ss_solver) PetscCallHYPRE(HYPRE_SStructSysPFMGDestroy(ex->ss_solver)); PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscFree(pc->data)); PetscFunctionReturn(PETSC_SUCCESS); } static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi", "Red/Black-Gauss-Seidel"}; static PetscErrorCode PCView_SysPFMG(PC pc, PetscViewer viewer) { PetscBool isascii; PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SysPFMG preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type])); PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_SysPFMG(PC pc, PetscOptionItems PetscOptionsObject) { PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; PetscBool flg = PETSC_FALSE; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options"); PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL)); if (flg) PetscCallHYPRE(HYPRE_SStructSysPFMGSetPrintLevel(ex->ss_solver, 3)); PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(ex->ss_solver, (HYPRE_Int)ex->its)); 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)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetNumPreRelax(ex->ss_solver, (HYPRE_Int)ex->num_pre_relax)); 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)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetNumPostRelax(ex->ss_solver, (HYPRE_Int)ex->num_post_relax)); PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(ex->ss_solver, ex->tol)); 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)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetRelaxType(ex->ss_solver, (HYPRE_Int)ex->relax_type)); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y) { PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; PetscScalar *yy; const PetscScalar *xx; PetscInt ilower[3], iupper[3]; HYPRE_Int hlower[3], hupper[3]; Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data; PetscInt ordering = mx->dofs_order; PetscInt nvars = mx->nvars; HYPRE_Int part = 0; PetscInt size; PetscInt i; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ iupper[0] += ilower[0] - 1; iupper[1] += ilower[1] - 1; iupper[2] += ilower[2] - 1; hlower[0] = (HYPRE_Int)ilower[0]; hlower[1] = (HYPRE_Int)ilower[1]; hlower[2] = (HYPRE_Int)ilower[2]; hupper[0] = (HYPRE_Int)iupper[0]; hupper[1] = (HYPRE_Int)iupper[1]; hupper[2] = (HYPRE_Int)iupper[2]; size = 1; for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1); /* copy x values over to hypre for variable ordering */ if (ordering) { PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0)); PetscCall(VecGetArrayRead(x, &xx)); for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(xx + (size * i)))); PetscCall(VecRestoreArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b)); PetscCallHYPRE(HYPRE_SStructMatrixMatvec(1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x)); PetscCallHYPRE(HYPRE_SStructSysPFMGSolve(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x)); /* copy solution values back to PETSc */ PetscCall(VecGetArray(y, &yy)); for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(yy + (size * i)))); PetscCall(VecRestoreArray(y, &yy)); } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ PetscScalar *z; PetscInt j, k; PetscCall(PetscMalloc1(nvars * size, &z)); PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0)); PetscCall(VecGetArrayRead(x, &xx)); /* transform nodal to hypre's variable ordering for sys_pfmg */ for (i = 0; i < size; i++) { k = i * nvars; for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j]; } for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i)))); PetscCall(VecRestoreArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b)); PetscCallHYPRE(HYPRE_SStructSysPFMGSolve(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x)); /* copy solution values back to PETSc */ PetscCall(VecGetArray(y, &yy)); for (i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i)))); /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ for (i = 0; i < size; i++) { k = i * nvars; for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i]; } PetscCall(VecRestoreArray(y, &yy)); PetscCall(PetscFree(z)); } PetscFunctionReturn(PETSC_SUCCESS); } 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) { PC_SysPFMG *jac = (PC_SysPFMG *)pc->data; HYPRE_Int oits; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(jac->ss_solver, (HYPRE_Int)(its * jac->its))); PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(jac->ss_solver, rtol)); PetscCall(PCApply_SysPFMG(pc, b, y)); PetscCallHYPRE(HYPRE_SStructSysPFMGGetNumIterations(jac->ss_solver, &oits)); *outits = oits; if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; else *reason = PCRICHARDSON_CONVERGED_RTOL; PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(jac->ss_solver, jac->tol)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(jac->ss_solver, (HYPRE_Int)jac->its)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetUp_SysPFMG(PC pc) { PC_SysPFMG *ex = (PC_SysPFMG *)pc->data; Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data; PetscBool flg; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESSTRUCT, &flg)); PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESSTRUCT with this preconditioner"); /* create the hypre sstruct solver object and set its information */ if (ex->ss_solver) PetscCallHYPRE(HYPRE_SStructSysPFMGDestroy(ex->ss_solver)); PetscCallHYPRE(HYPRE_SStructSysPFMGCreate(ex->hcomm, &ex->ss_solver)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetZeroGuess(ex->ss_solver)); PetscCallHYPRE(HYPRE_SStructSysPFMGSetup(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x)); PetscFunctionReturn(PETSC_SUCCESS); } /*MC PCSYSPFMG - the hypre SysPFMG multigrid solver Level: advanced Options Database Keys: + -pc_syspfmg_its - number of iterations of SysPFMG to use as preconditioner . -pc_syspfmg_num_pre_relax - number of smoothing steps before coarse grid . -pc_syspfmg_num_post_relax - number of smoothing steps after coarse grid . -pc_syspfmg_tol - tolerance of SysPFMG - -pc_syspfmg_relax_type - relaxation type for the up and down cycles Notes: See `PCPFMG` for hypre's PFMG that works for a scalar PDE and `PCSMG` See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver This is for CELL-centered descretizations This must be used with the `MATHYPRESSTRUCT` matrix type. 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`. .seealso: [](ch_ksp), `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG` M*/ PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) { PC_SysPFMG *ex; PetscFunctionBegin; PetscCall(PetscNew(&ex)); pc->data = ex; ex->its = 1; ex->tol = 1.e-8; ex->relax_type = 1; ex->num_pre_relax = 1; ex->num_post_relax = 1; pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; pc->ops->view = PCView_SysPFMG; pc->ops->destroy = PCDestroy_SysPFMG; pc->ops->apply = PCApply_SysPFMG; pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; pc->ops->setup = PCSetUp_SysPFMG; PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscHYPREInitialize()); PetscCallHYPRE(HYPRE_SStructSysPFMGCreate(ex->hcomm, &ex->ss_solver)); PetscFunctionReturn(PETSC_SUCCESS); } /* PC SMG */ typedef struct { MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ HYPRE_StructSolver hsolver; PetscInt its; /* keep copy of SMG options used so may view them */ PetscReal tol; PetscBool print_statistics; PetscInt num_pre_relax, num_post_relax; } PC_SMG; static PetscErrorCode PCDestroy_SMG(PC pc) { PC_SMG *ex = (PC_SMG *)pc->data; PetscFunctionBegin; if (ex->hsolver) PetscCallHYPRE(HYPRE_StructSMGDestroy(ex->hsolver)); PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscFree(pc->data)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer) { PetscBool isascii; PC_SMG *ex = (PC_SMG *)pc->data; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); if (isascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SMG preconditioning\n")); PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its)); PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol)); PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax)); } PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems PetscOptionsObject) { PC_SMG *ex = (PC_SMG *)pc->data; PetscFunctionBegin; PetscOptionsHeadBegin(PetscOptionsObject, "SMG options"); PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL)); 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)); 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)); PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL)); PetscOptionsHeadEnd(); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y) { PC_SMG *ex = (PC_SMG *)pc->data; PetscScalar *yy; const PetscScalar *xx; PetscInt ilower[3], iupper[3]; HYPRE_Int hlower[3], hupper[3]; Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2])); /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ iupper[0] += ilower[0] - 1; iupper[1] += ilower[1] - 1; iupper[2] += ilower[2] - 1; hlower[0] = (HYPRE_Int)ilower[0]; hlower[1] = (HYPRE_Int)ilower[1]; hlower[2] = (HYPRE_Int)ilower[2]; hupper[0] = (HYPRE_Int)iupper[0]; hupper[1] = (HYPRE_Int)iupper[1]; hupper[2] = (HYPRE_Int)iupper[2]; /* copy x values over to hypre */ PetscCallHYPRE(HYPRE_StructVectorSetConstantValues(mx->hb, 0.0)); PetscCall(VecGetArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_StructVectorSetBoxValues(mx->hb, hlower, hupper, (HYPRE_Complex *)xx)); PetscCall(VecRestoreArrayRead(x, &xx)); PetscCallHYPRE(HYPRE_StructVectorAssemble(mx->hb)); PetscCallHYPRE(HYPRE_StructSMGSolve(ex->hsolver, mx->hmat, mx->hb, mx->hx)); /* copy solution values back to PETSc */ PetscCall(VecGetArray(y, &yy)); PetscCallHYPRE(HYPRE_StructVectorGetBoxValues(mx->hx, hlower, hupper, (HYPRE_Complex *)yy)); PetscCall(VecRestoreArray(y, &yy)); PetscFunctionReturn(PETSC_SUCCESS); } 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) { PC_SMG *jac = (PC_SMG *)pc->data; HYPRE_Int oits; PetscFunctionBegin; PetscCall(PetscCitationsRegister(hypreCitation, &cite)); PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->its))); PetscCallHYPRE(HYPRE_StructSMGSetTol(jac->hsolver, rtol)); PetscCall(PCApply_SMG(pc, b, y)); PetscCallHYPRE(HYPRE_StructSMGGetNumIterations(jac->hsolver, &oits)); *outits = oits; if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; else *reason = PCRICHARDSON_CONVERGED_RTOL; PetscCallHYPRE(HYPRE_StructSMGSetTol(jac->hsolver, jac->tol)); PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->its)); PetscFunctionReturn(PETSC_SUCCESS); } static PetscErrorCode PCSetUp_SMG(PC pc) { PetscInt i, dim; PC_SMG *ex = (PC_SMG *)pc->data; Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data; PetscBool flg; DMBoundaryType p[3]; PetscInt M[3]; PetscFunctionBegin; PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg)); PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner"); PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0)); // Check if power of 2 in periodic directions for (i = 0; i < dim; i++) { PetscCheck((M[i] & (M[i] - 1)) == 0 || p[i] != DM_BOUNDARY_PERIODIC, 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]); } /* create the hypre solver object and set its information */ if (ex->hsolver) PetscCallHYPRE(HYPRE_StructSMGDestroy(ex->hsolver)); PetscCallHYPRE(HYPRE_StructSMGCreate(ex->hcomm, &ex->hsolver)); // The hypre options must be set here and not in SetFromOptions because it is created here! PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its)); PetscCallHYPRE(HYPRE_StructSMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax)); PetscCallHYPRE(HYPRE_StructSMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax)); PetscCallHYPRE(HYPRE_StructSMGSetTol(ex->hsolver, ex->tol)); PetscCallHYPRE(HYPRE_StructSMGSetup(ex->hsolver, mx->hmat, mx->hb, mx->hx)); PetscCallHYPRE(HYPRE_StructSMGSetZeroGuess(ex->hsolver)); PetscFunctionReturn(PETSC_SUCCESS); } /*MC PCSMG - the hypre (structured grid) SMG multigrid solver Level: advanced Options Database Keys: + -pc_smg_its - number of iterations of SMG to use as preconditioner . -pc_smg_num_pre_relax - number of smoothing steps before coarse grid . -pc_smg_num_post_relax - number of smoothing steps after coarse grid - -pc_smg_tol - tolerance of SMG Notes: This is for CELL-centered descretizations This must be used with the `MATHYPRESTRUCT` `MatType`. This does not provide all the functionality of hypre's SMG solver, it supports only one block per process defined by a PETSc `DMDA`. See `PCSYSPFMG`, `PCSMG`, `PCPFMG`, and `PCHYPRE` for access to hypre's other preconditioners .seealso: `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG` M*/ PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc) { PC_SMG *ex; PetscFunctionBegin; PetscCall(PetscNew(&ex)); pc->data = ex; ex->its = 1; ex->tol = 1.e-8; ex->num_pre_relax = 1; ex->num_post_relax = 1; pc->ops->setfromoptions = PCSetFromOptions_SMG; pc->ops->view = PCView_SMG; pc->ops->destroy = PCDestroy_SMG; pc->ops->apply = PCApply_SMG; pc->ops->applyrichardson = PCApplyRichardson_SMG; pc->ops->setup = PCSetUp_SMG; PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm)); PetscCall(PetscHYPREInitialize()); PetscCallHYPRE(HYPRE_StructSMGCreate(ex->hcomm, &ex->hsolver)); PetscFunctionReturn(PETSC_SUCCESS); }