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