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