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 %g must be greater than or equal to zero",(double)jac->pmax); 704 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 705 } 706 707 ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);CHKERRQ(ierr); 708 if (flg) { 709 if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl); 710 711 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 712 } 713 714 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); 715 if (flg) { 716 if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths); 717 718 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 719 } 720 721 ierr = PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);CHKERRQ(ierr); 722 if (flg) { 723 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); 724 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 725 } 726 ierr = PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);CHKERRQ(ierr); 727 if (flg) { 728 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); 729 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); 730 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 731 } 732 733 /* Grid sweeps */ 734 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); 735 if (flg) { 736 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx)); 737 /* modify the jac structure so we can view the updated options with PC_View */ 738 jac->gridsweeps[0] = indx; 739 jac->gridsweeps[1] = indx; 740 /*defaults coarse to 1 */ 741 jac->gridsweeps[2] = 1; 742 } 743 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); 744 if (flg) { 745 PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening)); 746 } 747 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); 748 if (flg) { 749 PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag)); 750 } 751 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); 752 if (flg) { 753 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant)); 754 } 755 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); 756 if (flg) { 757 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax)); 758 } 759 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); 760 if (flg) { 761 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth)); 762 } 763 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); 764 if (flg) { 765 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine)); 766 } 767 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);CHKERRQ(ierr); 768 if (flg) { 769 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1)); 770 jac->gridsweeps[0] = indx; 771 } 772 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);CHKERRQ(ierr); 773 if (flg) { 774 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2)); 775 jac->gridsweeps[1] = indx; 776 } 777 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);CHKERRQ(ierr); 778 if (flg) { 779 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3)); 780 jac->gridsweeps[2] = indx; 781 } 782 783 /* Smooth type */ 784 ierr = PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);CHKERRQ(ierr); 785 if (flg) { 786 jac->smoothtype = indx; 787 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6)); 788 jac->smoothnumlevels = 25; 789 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25)); 790 } 791 792 /* Number of smoothing levels */ 793 ierr = PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);CHKERRQ(ierr); 794 if (flg && (jac->smoothtype != -1)) { 795 jac->smoothnumlevels = indx; 796 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx)); 797 } 798 799 /* Number of levels for ILU(k) for Euclid */ 800 ierr = PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);CHKERRQ(ierr); 801 if (flg && (jac->smoothtype == 3)) { 802 jac->eu_level = indx; 803 PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx)); 804 } 805 806 /* Filter for ILU(k) for Euclid */ 807 double droptolerance; 808 ierr = PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);CHKERRQ(ierr); 809 if (flg && (jac->smoothtype == 3)) { 810 jac->eu_droptolerance = droptolerance; 811 PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance)); 812 } 813 814 /* Use Block Jacobi ILUT for Euclid */ 815 ierr = PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 816 if (flg && (jac->smoothtype == 3)) { 817 jac->eu_bj = tmp_truth; 818 PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj)); 819 } 820 821 /* Relax type */ 822 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); 823 if (flg) { 824 jac->relaxtype[0] = jac->relaxtype[1] = indx; 825 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx)); 826 /* by default, coarse type set to 9 */ 827 jac->relaxtype[2] = 9; 828 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3)); 829 } 830 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr); 831 if (flg) { 832 jac->relaxtype[0] = indx; 833 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1)); 834 } 835 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr); 836 if (flg) { 837 jac->relaxtype[1] = indx; 838 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2)); 839 } 840 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);CHKERRQ(ierr); 841 if (flg) { 842 jac->relaxtype[2] = indx; 843 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3)); 844 } 845 846 /* Relaxation Weight */ 847 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); 848 if (flg) { 849 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl)); 850 jac->relaxweight = tmpdbl; 851 } 852 853 n = 2; 854 twodbl[0] = twodbl[1] = 1.0; 855 ierr = PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr); 856 if (flg) { 857 if (n == 2) { 858 indx = (int)PetscAbsReal(twodbl[1]); 859 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx)); 860 } 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); 861 } 862 863 /* Outer relaxation Weight */ 864 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); 865 if (flg) { 866 PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl)); 867 jac->outerrelaxweight = tmpdbl; 868 } 869 870 n = 2; 871 twodbl[0] = twodbl[1] = 1.0; 872 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); 873 if (flg) { 874 if (n == 2) { 875 indx = (int)PetscAbsReal(twodbl[1]); 876 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx)); 877 } 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); 878 } 879 880 /* the Relax Order */ 881 ierr = PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 882 883 if (flg && tmp_truth) { 884 jac->relaxorder = 0; 885 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 886 } 887 ierr = PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);CHKERRQ(ierr); 888 if (flg) { 889 jac->measuretype = indx; 890 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 891 } 892 /* update list length 3/07 */ 893 ierr = PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);CHKERRQ(ierr); 894 if (flg) { 895 jac->coarsentype = indx; 896 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 897 } 898 899 ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg);CHKERRQ(ierr); 900 if (flg) { 901 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc)); 902 } 903 ierr = PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg);CHKERRQ(ierr); 904 if (flg) { 905 PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc)); 906 } 907 908 /* AIR */ 909 #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0) 910 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); 911 PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype)); 912 if (jac->Rtype) { 913 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 */ 914 915 ierr = PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR","Threshold for R","None",jac->Rstrongthreshold,&jac->Rstrongthreshold,NULL);CHKERRQ(ierr); 916 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold)); 917 918 ierr = PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR","Filter threshold for R","None",jac->Rfilterthreshold,&jac->Rfilterthreshold,NULL);CHKERRQ(ierr); 919 PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold)); 920 921 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); 922 PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol)); 923 924 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); 925 PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype)); 926 } 927 #endif 928 929 /* new 3/07 */ 930 ierr = PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);CHKERRQ(ierr); 931 if (flg || jac->Rtype) { 932 if (flg) jac->interptype = indx; 933 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 934 } 935 936 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);CHKERRQ(ierr); 937 if (flg) { 938 level = 3; 939 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);CHKERRQ(ierr); 940 941 jac->printstatistics = PETSC_TRUE; 942 PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level)); 943 } 944 945 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);CHKERRQ(ierr); 946 if (flg) { 947 level = 3; 948 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);CHKERRQ(ierr); 949 950 jac->printstatistics = PETSC_TRUE; 951 PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level)); 952 } 953 954 ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr); 955 if (flg && tmp_truth) { 956 PetscInt tmp_int; 957 ierr = PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);CHKERRQ(ierr); 958 if (flg) jac->nodal_relax_levels = tmp_int; 959 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6)); 960 PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1)); 961 PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0)); 962 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels)); 963 } 964 965 ierr = PetscOptionsBool("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL);CHKERRQ(ierr); 966 PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,(jac->hsolver,jac->keeptranspose ? 1 : 0)); 967 968 /* options for ParaSails solvers */ 969 ierr = PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flg);CHKERRQ(ierr); 970 if (flg) { 971 jac->symt = indx; 972 PetscStackCallStandard(HYPRE_BoomerAMGSetSym,(jac->hsolver,jac->symt)); 973 } 974 975 ierr = PetscOptionsTail();CHKERRQ(ierr); 976 PetscFunctionReturn(0); 977 } 978 979 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) 980 { 981 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 982 PetscErrorCode ierr; 983 HYPRE_Int oits; 984 985 PetscFunctionBegin; 986 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 987 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter)); 988 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol)); 989 jac->applyrichardson = PETSC_TRUE; 990 ierr = PCApply_HYPRE(pc,b,y);CHKERRQ(ierr); 991 jac->applyrichardson = PETSC_FALSE; 992 PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits)); 993 *outits = oits; 994 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 995 else *reason = PCRICHARDSON_CONVERGED_RTOL; 996 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 997 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 998 PetscFunctionReturn(0); 999 } 1000 1001 static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer) 1002 { 1003 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1004 PetscErrorCode ierr; 1005 PetscBool iascii; 1006 1007 PetscFunctionBegin; 1008 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1009 if (iascii) { 1010 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");CHKERRQ(ierr); 1011 ierr = PetscViewerASCIIPrintf(viewer," Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);CHKERRQ(ierr); 1012 ierr = PetscViewerASCIIPrintf(viewer," Maximum number of levels %D\n",jac->maxlevels);CHKERRQ(ierr); 1013 ierr = PetscViewerASCIIPrintf(viewer," Maximum number of iterations PER hypre call %D\n",jac->maxiter);CHKERRQ(ierr); 1014 ierr = PetscViewerASCIIPrintf(viewer," Convergence tolerance PER hypre call %g\n",(double)jac->tol);CHKERRQ(ierr); 1015 ierr = PetscViewerASCIIPrintf(viewer," Threshold for strong coupling %g\n",(double)jac->strongthreshold);CHKERRQ(ierr); 1016 ierr = PetscViewerASCIIPrintf(viewer," Interpolation truncation factor %g\n",(double)jac->truncfactor);CHKERRQ(ierr); 1017 ierr = PetscViewerASCIIPrintf(viewer," Interpolation: max elements per row %D\n",jac->pmax);CHKERRQ(ierr); 1018 if (jac->interp_refine) { 1019 ierr = PetscViewerASCIIPrintf(viewer," Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);CHKERRQ(ierr); 1020 } 1021 ierr = PetscViewerASCIIPrintf(viewer," Number of levels of aggressive coarsening %D\n",jac->agg_nl);CHKERRQ(ierr); 1022 ierr = PetscViewerASCIIPrintf(viewer," Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);CHKERRQ(ierr); 1023 1024 ierr = PetscViewerASCIIPrintf(viewer," Maximum row sums %g\n",(double)jac->maxrowsum);CHKERRQ(ierr); 1025 1026 ierr = PetscViewerASCIIPrintf(viewer," Sweeps down %D\n",jac->gridsweeps[0]);CHKERRQ(ierr); 1027 ierr = PetscViewerASCIIPrintf(viewer," Sweeps up %D\n",jac->gridsweeps[1]);CHKERRQ(ierr); 1028 ierr = PetscViewerASCIIPrintf(viewer," Sweeps on coarse %D\n",jac->gridsweeps[2]);CHKERRQ(ierr); 1029 1030 ierr = PetscViewerASCIIPrintf(viewer," Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);CHKERRQ(ierr); 1031 ierr = PetscViewerASCIIPrintf(viewer," Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);CHKERRQ(ierr); 1032 ierr = PetscViewerASCIIPrintf(viewer," Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);CHKERRQ(ierr); 1033 1034 ierr = PetscViewerASCIIPrintf(viewer," Relax weight (all) %g\n",(double)jac->relaxweight);CHKERRQ(ierr); 1035 ierr = PetscViewerASCIIPrintf(viewer," Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);CHKERRQ(ierr); 1036 1037 if (jac->relaxorder) { 1038 ierr = PetscViewerASCIIPrintf(viewer," Using CF-relaxation\n");CHKERRQ(ierr); 1039 } else { 1040 ierr = PetscViewerASCIIPrintf(viewer," Not using CF-relaxation\n");CHKERRQ(ierr); 1041 } 1042 if (jac->smoothtype!=-1) { 1043 ierr = PetscViewerASCIIPrintf(viewer," Smooth type %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);CHKERRQ(ierr); 1044 ierr = PetscViewerASCIIPrintf(viewer," Smooth num levels %D\n",jac->smoothnumlevels);CHKERRQ(ierr); 1045 } else { 1046 ierr = PetscViewerASCIIPrintf(viewer," Not using more complex smoothers.\n");CHKERRQ(ierr); 1047 } 1048 if (jac->smoothtype==3) { 1049 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU(k) levels %D\n",jac->eu_level);CHKERRQ(ierr); 1050 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);CHKERRQ(ierr); 1051 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);CHKERRQ(ierr); 1052 } 1053 ierr = PetscViewerASCIIPrintf(viewer," Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);CHKERRQ(ierr); 1054 ierr = PetscViewerASCIIPrintf(viewer," Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);CHKERRQ(ierr); 1055 ierr = PetscViewerASCIIPrintf(viewer," Interpolation type %s\n",jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt");CHKERRQ(ierr); 1056 if (jac->nodal_coarsening) { 1057 ierr = PetscViewerASCIIPrintf(viewer," Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);CHKERRQ(ierr); 1058 } 1059 if (jac->vec_interp_variant) { 1060 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);CHKERRQ(ierr); 1061 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);CHKERRQ(ierr); 1062 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);CHKERRQ(ierr); 1063 } 1064 if (jac->nodal_relax) { 1065 ierr = PetscViewerASCIIPrintf(viewer," Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);CHKERRQ(ierr); 1066 } 1067 1068 /* AIR */ 1069 if (jac->Rtype) { 1070 ierr = PetscViewerASCIIPrintf(viewer," Using approximate ideal restriction type %D\n",jac->Rtype);CHKERRQ(ierr); 1071 ierr = PetscViewerASCIIPrintf(viewer," Threshold for R %g\n",(double)jac->Rstrongthreshold);CHKERRQ(ierr); 1072 ierr = PetscViewerASCIIPrintf(viewer," Filter for R %g\n",(double)jac->Rfilterthreshold);CHKERRQ(ierr); 1073 ierr = PetscViewerASCIIPrintf(viewer," A drop tolerance %g\n",(double)jac->Adroptol);CHKERRQ(ierr); 1074 ierr = PetscViewerASCIIPrintf(viewer," A drop type %D\n",jac->Adroptype);CHKERRQ(ierr); 1075 } 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 /* --------------------------------------------------------------------------------------------*/ 1081 static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc) 1082 { 1083 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1084 PetscErrorCode ierr; 1085 PetscInt indx; 1086 PetscBool flag; 1087 const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"}; 1088 1089 PetscFunctionBegin; 1090 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");CHKERRQ(ierr); 1091 ierr = PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);CHKERRQ(ierr); 1092 ierr = PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshold,&jac->threshold,&flag);CHKERRQ(ierr); 1093 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels)); 1094 1095 ierr = PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);CHKERRQ(ierr); 1096 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 1097 1098 ierr = PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);CHKERRQ(ierr); 1099 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 1100 1101 ierr = PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);CHKERRQ(ierr); 1102 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 1103 1104 ierr = PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);CHKERRQ(ierr); 1105 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 1106 1107 ierr = PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);CHKERRQ(ierr); 1108 if (flag) { 1109 jac->symt = indx; 1110 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 1111 } 1112 1113 ierr = PetscOptionsTail();CHKERRQ(ierr); 1114 PetscFunctionReturn(0); 1115 } 1116 1117 static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer) 1118 { 1119 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1120 PetscErrorCode ierr; 1121 PetscBool iascii; 1122 const char *symt = 0; 1123 1124 PetscFunctionBegin; 1125 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1126 if (iascii) { 1127 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");CHKERRQ(ierr); 1128 ierr = PetscViewerASCIIPrintf(viewer," nlevels %d\n",jac->nlevels);CHKERRQ(ierr); 1129 ierr = PetscViewerASCIIPrintf(viewer," threshold %g\n",(double)jac->threshold);CHKERRQ(ierr); 1130 ierr = PetscViewerASCIIPrintf(viewer," filter %g\n",(double)jac->filter);CHKERRQ(ierr); 1131 ierr = PetscViewerASCIIPrintf(viewer," load balance %g\n",(double)jac->loadbal);CHKERRQ(ierr); 1132 ierr = PetscViewerASCIIPrintf(viewer," reuse nonzero structure %s\n",PetscBools[jac->ruse]);CHKERRQ(ierr); 1133 ierr = PetscViewerASCIIPrintf(viewer," print info to screen %s\n",PetscBools[jac->logging]);CHKERRQ(ierr); 1134 if (!jac->symt) symt = "nonsymmetric matrix and preconditioner"; 1135 else if (jac->symt == 1) symt = "SPD matrix and preconditioner"; 1136 else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner"; 1137 else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt); 1138 ierr = PetscViewerASCIIPrintf(viewer," %s\n",symt);CHKERRQ(ierr); 1139 } 1140 PetscFunctionReturn(0); 1141 } 1142 /* --------------------------------------------------------------------------------------------*/ 1143 static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc) 1144 { 1145 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1146 PetscErrorCode ierr; 1147 PetscInt n; 1148 PetscBool flag,flag2,flag3,flag4; 1149 1150 PetscFunctionBegin; 1151 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");CHKERRQ(ierr); 1152 ierr = PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr); 1153 if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print)); 1154 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); 1155 if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1156 ierr = PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);CHKERRQ(ierr); 1157 if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1158 ierr = PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr); 1159 if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol)); 1160 ierr = PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr); 1161 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); 1162 ierr = PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr); 1163 ierr = PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr); 1164 if (flag || flag2 || flag3 || flag4) { 1165 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1166 jac->as_relax_times, 1167 jac->as_relax_weight, 1168 jac->as_omega)); 1169 } 1170 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); 1171 n = 5; 1172 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr); 1173 if (flag || flag2) { 1174 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1175 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1176 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1177 jac->as_amg_alpha_theta, 1178 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1179 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1180 } 1181 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); 1182 n = 5; 1183 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);CHKERRQ(ierr); 1184 if (flag || flag2) { 1185 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1186 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1187 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1188 jac->as_amg_beta_theta, 1189 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1190 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1191 } 1192 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); 1193 if (flag) { /* override HYPRE's default only if the options is used */ 1194 PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq)); 1195 } 1196 ierr = PetscOptionsTail();CHKERRQ(ierr); 1197 PetscFunctionReturn(0); 1198 } 1199 1200 static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer) 1201 { 1202 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1203 PetscErrorCode ierr; 1204 PetscBool iascii; 1205 1206 PetscFunctionBegin; 1207 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1208 if (iascii) { 1209 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS preconditioning\n");CHKERRQ(ierr); 1210 ierr = PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr); 1211 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr); 1212 ierr = PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr); 1213 ierr = PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr); 1214 ierr = PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr); 1215 ierr = PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr); 1216 ierr = PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);CHKERRQ(ierr); 1217 if (jac->alpha_Poisson) { 1218 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver (passed in by user)\n");CHKERRQ(ierr); 1219 } else { 1220 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver (computed) \n");CHKERRQ(ierr); 1221 } 1222 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr); 1223 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr); 1224 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr); 1225 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr); 1226 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr); 1227 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr); 1228 if (!jac->ams_beta_is_zero) { 1229 if (jac->beta_Poisson) { 1230 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver (passed in by user)\n");CHKERRQ(ierr); 1231 } else { 1232 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver (computed) \n");CHKERRQ(ierr); 1233 } 1234 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr); 1235 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr); 1236 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr); 1237 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr); 1238 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr); 1239 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr); 1240 if (jac->ams_beta_is_zero_part) { 1241 ierr = PetscViewerASCIIPrintf(viewer," compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);CHKERRQ(ierr); 1242 } 1243 } else { 1244 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver not used (zero-conductivity everywhere) \n");CHKERRQ(ierr); 1245 } 1246 } 1247 PetscFunctionReturn(0); 1248 } 1249 1250 static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc) 1251 { 1252 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1253 PetscErrorCode ierr; 1254 PetscInt n; 1255 PetscBool flag,flag2,flag3,flag4; 1256 1257 PetscFunctionBegin; 1258 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");CHKERRQ(ierr); 1259 ierr = PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr); 1260 if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print)); 1261 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); 1262 if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1263 ierr = PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);CHKERRQ(ierr); 1264 if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type)); 1265 ierr = PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr); 1266 if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol)); 1267 ierr = PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr); 1268 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); 1269 ierr = PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr); 1270 ierr = PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr); 1271 if (flag || flag2 || flag3 || flag4) { 1272 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1273 jac->as_relax_times, 1274 jac->as_relax_weight, 1275 jac->as_omega)); 1276 } 1277 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); 1278 n = 5; 1279 ierr = PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr); 1280 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); 1281 if (flag || flag2 || flag3) { 1282 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMS cycle type */ 1283 jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1284 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1285 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1286 jac->as_amg_alpha_theta, 1287 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1288 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1289 } 1290 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); 1291 n = 5; 1292 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); 1293 if (flag || flag2) { 1294 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1295 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1296 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1297 jac->as_amg_beta_theta, 1298 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1299 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1300 } 1301 ierr = PetscOptionsTail();CHKERRQ(ierr); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer) 1306 { 1307 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1308 PetscErrorCode ierr; 1309 PetscBool iascii; 1310 1311 PetscFunctionBegin; 1312 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1313 if (iascii) { 1314 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ADS preconditioning\n");CHKERRQ(ierr); 1315 ierr = PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr); 1316 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ads_cycle_type);CHKERRQ(ierr); 1317 ierr = PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr); 1318 ierr = PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr); 1319 ierr = PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr); 1320 ierr = PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr); 1321 ierr = PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);CHKERRQ(ierr); 1322 ierr = PetscViewerASCIIPrintf(viewer," AMS solver using boomerAMG\n");CHKERRQ(ierr); 1323 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr); 1324 ierr = PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr); 1325 ierr = PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr); 1326 ierr = PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr); 1327 ierr = PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr); 1328 ierr = PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr); 1329 ierr = PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr); 1330 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver using boomerAMG\n");CHKERRQ(ierr); 1331 ierr = PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr); 1332 ierr = PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr); 1333 ierr = PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr); 1334 ierr = PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr); 1335 ierr = PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr); 1336 ierr = PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr); 1337 } 1338 PetscFunctionReturn(0); 1339 } 1340 1341 static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G) 1342 { 1343 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1344 PetscBool ishypre; 1345 PetscErrorCode ierr; 1346 1347 PetscFunctionBegin; 1348 ierr = PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);CHKERRQ(ierr); 1349 if (ishypre) { 1350 ierr = PetscObjectReference((PetscObject)G);CHKERRQ(ierr); 1351 ierr = MatDestroy(&jac->G);CHKERRQ(ierr); 1352 jac->G = G; 1353 } else { 1354 ierr = MatDestroy(&jac->G);CHKERRQ(ierr); 1355 ierr = MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);CHKERRQ(ierr); 1356 } 1357 PetscFunctionReturn(0); 1358 } 1359 1360 /*@ 1361 PCHYPRESetDiscreteGradient - Set discrete gradient matrix 1362 1363 Collective on PC 1364 1365 Input Parameters: 1366 + pc - the preconditioning context 1367 - G - the discrete gradient 1368 1369 Level: intermediate 1370 1371 Notes: 1372 G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh 1373 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 1374 1375 .seealso: 1376 @*/ 1377 PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G) 1378 { 1379 PetscErrorCode ierr; 1380 1381 PetscFunctionBegin; 1382 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1383 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 1384 PetscCheckSameComm(pc,1,G,2); 1385 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));CHKERRQ(ierr); 1386 PetscFunctionReturn(0); 1387 } 1388 1389 static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C) 1390 { 1391 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1392 PetscBool ishypre; 1393 PetscErrorCode ierr; 1394 1395 PetscFunctionBegin; 1396 ierr = PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);CHKERRQ(ierr); 1397 if (ishypre) { 1398 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 1399 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1400 jac->C = C; 1401 } else { 1402 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1403 ierr = MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 1404 } 1405 PetscFunctionReturn(0); 1406 } 1407 1408 /*@ 1409 PCHYPRESetDiscreteCurl - Set discrete curl matrix 1410 1411 Collective on PC 1412 1413 Input Parameters: 1414 + pc - the preconditioning context 1415 - C - the discrete curl 1416 1417 Level: intermediate 1418 1419 Notes: 1420 C should have as many rows as the number of faces and as many columns as the number of edges in the mesh 1421 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 1422 1423 .seealso: 1424 @*/ 1425 PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C) 1426 { 1427 PetscErrorCode ierr; 1428 1429 PetscFunctionBegin; 1430 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1431 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 1432 PetscCheckSameComm(pc,1,C,2); 1433 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));CHKERRQ(ierr); 1434 PetscFunctionReturn(0); 1435 } 1436 1437 static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) 1438 { 1439 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1440 PetscBool ishypre; 1441 PetscErrorCode ierr; 1442 PetscInt i; 1443 PetscFunctionBegin; 1444 1445 ierr = MatDestroy(&jac->RT_PiFull);CHKERRQ(ierr); 1446 ierr = MatDestroy(&jac->ND_PiFull);CHKERRQ(ierr); 1447 for (i=0;i<3;++i) { 1448 ierr = MatDestroy(&jac->RT_Pi[i]);CHKERRQ(ierr); 1449 ierr = MatDestroy(&jac->ND_Pi[i]);CHKERRQ(ierr); 1450 } 1451 1452 jac->dim = dim; 1453 if (RT_PiFull) { 1454 ierr = PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr); 1455 if (ishypre) { 1456 ierr = PetscObjectReference((PetscObject)RT_PiFull);CHKERRQ(ierr); 1457 jac->RT_PiFull = RT_PiFull; 1458 } else { 1459 ierr = MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);CHKERRQ(ierr); 1460 } 1461 } 1462 if (RT_Pi) { 1463 for (i=0;i<dim;++i) { 1464 if (RT_Pi[i]) { 1465 ierr = PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr); 1466 if (ishypre) { 1467 ierr = PetscObjectReference((PetscObject)RT_Pi[i]);CHKERRQ(ierr); 1468 jac->RT_Pi[i] = RT_Pi[i]; 1469 } else { 1470 ierr = MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);CHKERRQ(ierr); 1471 } 1472 } 1473 } 1474 } 1475 if (ND_PiFull) { 1476 ierr = PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr); 1477 if (ishypre) { 1478 ierr = PetscObjectReference((PetscObject)ND_PiFull);CHKERRQ(ierr); 1479 jac->ND_PiFull = ND_PiFull; 1480 } else { 1481 ierr = MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);CHKERRQ(ierr); 1482 } 1483 } 1484 if (ND_Pi) { 1485 for (i=0;i<dim;++i) { 1486 if (ND_Pi[i]) { 1487 ierr = PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr); 1488 if (ishypre) { 1489 ierr = PetscObjectReference((PetscObject)ND_Pi[i]);CHKERRQ(ierr); 1490 jac->ND_Pi[i] = ND_Pi[i]; 1491 } else { 1492 ierr = MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);CHKERRQ(ierr); 1493 } 1494 } 1495 } 1496 } 1497 1498 PetscFunctionReturn(0); 1499 } 1500 1501 /*@ 1502 PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner 1503 1504 Collective on PC 1505 1506 Input Parameters: 1507 + pc - the preconditioning context 1508 - dim - the dimension of the problem, only used in AMS 1509 - RT_PiFull - Raviart-Thomas interpolation matrix 1510 - RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix 1511 - ND_PiFull - Nedelec interpolation matrix 1512 - ND_Pi - x/y/z component of Nedelec interpolation matrix 1513 1514 Notes: 1515 For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL. 1516 For ADS, both type of interpolation matrices are needed. 1517 Level: intermediate 1518 1519 .seealso: 1520 @*/ 1521 PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) 1522 { 1523 PetscErrorCode ierr; 1524 PetscInt i; 1525 1526 PetscFunctionBegin; 1527 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1528 if (RT_PiFull) { 1529 PetscValidHeaderSpecific(RT_PiFull,MAT_CLASSID,3); 1530 PetscCheckSameComm(pc,1,RT_PiFull,3); 1531 } 1532 if (RT_Pi) { 1533 PetscValidPointer(RT_Pi,4); 1534 for (i=0;i<dim;++i) { 1535 if (RT_Pi[i]) { 1536 PetscValidHeaderSpecific(RT_Pi[i],MAT_CLASSID,4); 1537 PetscCheckSameComm(pc,1,RT_Pi[i],4); 1538 } 1539 } 1540 } 1541 if (ND_PiFull) { 1542 PetscValidHeaderSpecific(ND_PiFull,MAT_CLASSID,5); 1543 PetscCheckSameComm(pc,1,ND_PiFull,5); 1544 } 1545 if (ND_Pi) { 1546 PetscValidPointer(ND_Pi,6); 1547 for (i=0;i<dim;++i) { 1548 if (ND_Pi[i]) { 1549 PetscValidHeaderSpecific(ND_Pi[i],MAT_CLASSID,6); 1550 PetscCheckSameComm(pc,1,ND_Pi[i],6); 1551 } 1552 } 1553 } 1554 ierr = PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));CHKERRQ(ierr); 1555 PetscFunctionReturn(0); 1556 } 1557 1558 static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha) 1559 { 1560 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1561 PetscBool ishypre; 1562 PetscErrorCode ierr; 1563 1564 PetscFunctionBegin; 1565 ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr); 1566 if (ishypre) { 1567 if (isalpha) { 1568 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 1569 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr); 1570 jac->alpha_Poisson = A; 1571 } else { 1572 if (A) { 1573 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 1574 } else { 1575 jac->ams_beta_is_zero = PETSC_TRUE; 1576 } 1577 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1578 jac->beta_Poisson = A; 1579 } 1580 } else { 1581 if (isalpha) { 1582 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr); 1583 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);CHKERRQ(ierr); 1584 } else { 1585 if (A) { 1586 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1587 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);CHKERRQ(ierr); 1588 } else { 1589 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1590 jac->ams_beta_is_zero = PETSC_TRUE; 1591 } 1592 } 1593 } 1594 PetscFunctionReturn(0); 1595 } 1596 1597 /*@ 1598 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix 1599 1600 Collective on PC 1601 1602 Input Parameters: 1603 + pc - the preconditioning context 1604 - A - the matrix 1605 1606 Level: intermediate 1607 1608 Notes: 1609 A should be obtained by discretizing the vector valued Poisson problem with linear finite elements 1610 1611 .seealso: 1612 @*/ 1613 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) 1614 { 1615 PetscErrorCode ierr; 1616 1617 PetscFunctionBegin; 1618 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1619 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 1620 PetscCheckSameComm(pc,1,A,2); 1621 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));CHKERRQ(ierr); 1622 PetscFunctionReturn(0); 1623 } 1624 1625 /*@ 1626 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix 1627 1628 Collective on PC 1629 1630 Input Parameters: 1631 + pc - the preconditioning context 1632 - A - the matrix 1633 1634 Level: intermediate 1635 1636 Notes: 1637 A should be obtained by discretizing the Poisson problem with linear finite elements. 1638 Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL. 1639 1640 .seealso: 1641 @*/ 1642 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) 1643 { 1644 PetscErrorCode ierr; 1645 1646 PetscFunctionBegin; 1647 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1648 if (A) { 1649 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 1650 PetscCheckSameComm(pc,1,A,2); 1651 } 1652 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));CHKERRQ(ierr); 1653 PetscFunctionReturn(0); 1654 } 1655 1656 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo) 1657 { 1658 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1659 PetscErrorCode ierr; 1660 1661 PetscFunctionBegin; 1662 /* throw away any vector if already set */ 1663 ierr = VecHYPRE_IJVectorDestroy(&jac->constants[0]);CHKERRQ(ierr); 1664 ierr = VecHYPRE_IJVectorDestroy(&jac->constants[1]);CHKERRQ(ierr); 1665 ierr = VecHYPRE_IJVectorDestroy(&jac->constants[2]);CHKERRQ(ierr); 1666 ierr = VecHYPRE_IJVectorCreate(ozz->map,&jac->constants[0]);CHKERRQ(ierr); 1667 ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr); 1668 ierr = VecHYPRE_IJVectorCreate(zoz->map,&jac->constants[1]);CHKERRQ(ierr); 1669 ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr); 1670 jac->dim = 2; 1671 if (zzo) { 1672 ierr = VecHYPRE_IJVectorCreate(zzo->map,&jac->constants[2]);CHKERRQ(ierr); 1673 ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr); 1674 jac->dim++; 1675 } 1676 PetscFunctionReturn(0); 1677 } 1678 1679 /*@ 1680 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis 1681 1682 Collective on PC 1683 1684 Input Parameters: 1685 + pc - the preconditioning context 1686 - ozz - vector representing (1,0,0) (or (1,0) in 2D) 1687 - zoz - vector representing (0,1,0) (or (0,1) in 2D) 1688 - zzo - vector representing (0,0,1) (use NULL in 2D) 1689 1690 Level: intermediate 1691 1692 Notes: 1693 1694 .seealso: 1695 @*/ 1696 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) 1697 { 1698 PetscErrorCode ierr; 1699 1700 PetscFunctionBegin; 1701 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1702 PetscValidHeaderSpecific(ozz,VEC_CLASSID,2); 1703 PetscValidHeaderSpecific(zoz,VEC_CLASSID,3); 1704 if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4); 1705 PetscCheckSameComm(pc,1,ozz,2); 1706 PetscCheckSameComm(pc,1,zoz,3); 1707 if (zzo) PetscCheckSameComm(pc,1,zzo,4); 1708 ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr); 1709 PetscFunctionReturn(0); 1710 } 1711 1712 static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 1713 { 1714 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1715 Vec tv; 1716 PetscInt i; 1717 PetscErrorCode ierr; 1718 1719 PetscFunctionBegin; 1720 /* throw away any coordinate vector if already set */ 1721 ierr = VecHYPRE_IJVectorDestroy(&jac->coords[0]);CHKERRQ(ierr); 1722 ierr = VecHYPRE_IJVectorDestroy(&jac->coords[1]);CHKERRQ(ierr); 1723 ierr = VecHYPRE_IJVectorDestroy(&jac->coords[2]);CHKERRQ(ierr); 1724 jac->dim = dim; 1725 1726 /* compute IJ vector for coordinates */ 1727 ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr); 1728 ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr); 1729 ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr); 1730 for (i=0;i<dim;i++) { 1731 PetscScalar *array; 1732 PetscInt j; 1733 1734 ierr = VecHYPRE_IJVectorCreate(tv->map,&jac->coords[i]);CHKERRQ(ierr); 1735 ierr = VecGetArrayWrite(tv,&array);CHKERRQ(ierr); 1736 for (j=0;j<nloc;j++) array[j] = coords[j*dim+i]; 1737 ierr = VecRestoreArrayWrite(tv,&array);CHKERRQ(ierr); 1738 ierr = VecHYPRE_IJVectorCopy(tv,jac->coords[i]);CHKERRQ(ierr); 1739 } 1740 ierr = VecDestroy(&tv);CHKERRQ(ierr); 1741 PetscFunctionReturn(0); 1742 } 1743 1744 /* ---------------------------------------------------------------------------------*/ 1745 1746 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[]) 1747 { 1748 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1749 1750 PetscFunctionBegin; 1751 *name = jac->hypre_type; 1752 PetscFunctionReturn(0); 1753 } 1754 1755 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[]) 1756 { 1757 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1758 PetscErrorCode ierr; 1759 PetscBool flag; 1760 1761 PetscFunctionBegin; 1762 if (jac->hypre_type) { 1763 ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr); 1764 if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set"); 1765 PetscFunctionReturn(0); 1766 } else { 1767 ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr); 1768 } 1769 1770 jac->maxiter = PETSC_DEFAULT; 1771 jac->tol = PETSC_DEFAULT; 1772 jac->printstatistics = PetscLogPrintInfo; 1773 1774 ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr); 1775 if (flag) { 1776 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRMPI(ierr); 1777 PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver)); 1778 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; 1779 pc->ops->view = PCView_HYPRE_Pilut; 1780 jac->destroy = HYPRE_ParCSRPilutDestroy; 1781 jac->setup = HYPRE_ParCSRPilutSetup; 1782 jac->solve = HYPRE_ParCSRPilutSolve; 1783 jac->factorrowsize = PETSC_DEFAULT; 1784 PetscFunctionReturn(0); 1785 } 1786 ierr = PetscStrcmp("euclid",jac->hypre_type,&flag);CHKERRQ(ierr); 1787 if (flag) { 1788 #if defined(PETSC_HAVE_64BIT_INDICES) 1789 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Hypre Euclid not support with 64 bit indices"); 1790 #endif 1791 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRMPI(ierr); 1792 PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver)); 1793 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid; 1794 pc->ops->view = PCView_HYPRE_Euclid; 1795 jac->destroy = HYPRE_EuclidDestroy; 1796 jac->setup = HYPRE_EuclidSetup; 1797 jac->solve = HYPRE_EuclidSolve; 1798 jac->factorrowsize = PETSC_DEFAULT; 1799 jac->eu_level = PETSC_DEFAULT; /* default */ 1800 PetscFunctionReturn(0); 1801 } 1802 ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr); 1803 if (flag) { 1804 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRMPI(ierr); 1805 PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver)); 1806 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; 1807 pc->ops->view = PCView_HYPRE_ParaSails; 1808 jac->destroy = HYPRE_ParaSailsDestroy; 1809 jac->setup = HYPRE_ParaSailsSetup; 1810 jac->solve = HYPRE_ParaSailsSolve; 1811 /* initialize */ 1812 jac->nlevels = 1; 1813 jac->threshold = .1; 1814 jac->filter = .1; 1815 jac->loadbal = 0; 1816 if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE; 1817 else jac->logging = (int) PETSC_FALSE; 1818 1819 jac->ruse = (int) PETSC_FALSE; 1820 jac->symt = 0; 1821 PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels)); 1822 PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 1823 PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 1824 PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 1825 PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 1826 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 1827 PetscFunctionReturn(0); 1828 } 1829 ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr); 1830 if (flag) { 1831 ierr = HYPRE_BoomerAMGCreate(&jac->hsolver); 1832 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; 1833 pc->ops->view = PCView_HYPRE_BoomerAMG; 1834 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; 1835 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; 1836 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);CHKERRQ(ierr); 1837 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);CHKERRQ(ierr); 1838 jac->destroy = HYPRE_BoomerAMGDestroy; 1839 jac->setup = HYPRE_BoomerAMGSetup; 1840 jac->solve = HYPRE_BoomerAMGSolve; 1841 jac->applyrichardson = PETSC_FALSE; 1842 /* these defaults match the hypre defaults */ 1843 jac->cycletype = 1; 1844 jac->maxlevels = 25; 1845 jac->maxiter = 1; 1846 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ 1847 jac->truncfactor = 0.0; 1848 jac->strongthreshold = .25; 1849 jac->maxrowsum = .9; 1850 jac->coarsentype = 6; 1851 jac->measuretype = 0; 1852 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; 1853 jac->smoothtype = -1; /* Not set by default */ 1854 jac->smoothnumlevels = 25; 1855 jac->eu_level = 0; 1856 jac->eu_droptolerance = 0; 1857 jac->eu_bj = 0; 1858 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */ 1859 jac->relaxtype[2] = 9; /*G.E. */ 1860 jac->relaxweight = 1.0; 1861 jac->outerrelaxweight = 1.0; 1862 jac->relaxorder = 1; 1863 jac->interptype = 0; 1864 jac->Rtype = 0; 1865 jac->Rstrongthreshold = 0.25; 1866 jac->Rfilterthreshold = 0.0; 1867 jac->Adroptype = -1; 1868 jac->Adroptol = 0.0; 1869 jac->agg_nl = 0; 1870 jac->agg_interptype = 4; 1871 jac->pmax = 0; 1872 jac->truncfactor = 0.0; 1873 jac->agg_num_paths = 1; 1874 jac->maxc = 9; 1875 jac->minc = 1; 1876 1877 jac->nodal_coarsening = 0; 1878 jac->nodal_coarsening_diag = 0; 1879 jac->vec_interp_variant = 0; 1880 jac->vec_interp_qmax = 0; 1881 jac->vec_interp_smooth = PETSC_FALSE; 1882 jac->interp_refine = 0; 1883 jac->nodal_relax = PETSC_FALSE; 1884 jac->nodal_relax_levels = 1; 1885 jac->rap2 = 0; 1886 1887 /* GPU defaults 1888 from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options 1889 and /src/parcsr_ls/par_amg.c */ 1890 #if defined(PETSC_HAVE_HYPRE_DEVICE) 1891 jac->keeptranspose = PETSC_TRUE; 1892 jac->mod_rap2 = 1; 1893 jac->coarsentype = 8; 1894 jac->relaxorder = 0; 1895 jac->interptype = 6; 1896 jac->relaxtype[0] = 18; 1897 jac->relaxtype[1] = 18; 1898 jac->agg_interptype = 7; 1899 #else 1900 jac->keeptranspose = PETSC_FALSE; 1901 jac->mod_rap2 = 0; 1902 #endif 1903 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 1904 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 1905 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 1906 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 1907 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 1908 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 1909 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 1910 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 1911 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 1912 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 1913 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 1914 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 1915 PetscStackCallStandard(HYPRE_BoomerAMGSetAggInterpType,(jac->hsolver,jac->agg_interptype)); 1916 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 1917 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 1918 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /* defaults coarse to 9 */ 1919 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /* defaults coarse to 1 */ 1920 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc)); 1921 PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc)); 1922 1923 /* GPU */ 1924 #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0) 1925 PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,(jac->hsolver,jac->keeptranspose ? 1 : 0)); 1926 PetscStackCallStandard(HYPRE_BoomerAMGSetRAP2,(jac->hsolver, jac->rap2)); 1927 PetscStackCallStandard(HYPRE_BoomerAMGSetModuleRAP2,(jac->hsolver, jac->mod_rap2)); 1928 #endif 1929 1930 /* AIR */ 1931 #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0) 1932 PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype)); 1933 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold)); 1934 PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold)); 1935 PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol)); 1936 PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype)); 1937 #endif 1938 PetscFunctionReturn(0); 1939 } 1940 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr); 1941 if (flag) { 1942 ierr = HYPRE_AMSCreate(&jac->hsolver); 1943 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 1944 pc->ops->view = PCView_HYPRE_AMS; 1945 jac->destroy = HYPRE_AMSDestroy; 1946 jac->setup = HYPRE_AMSSetup; 1947 jac->solve = HYPRE_AMSSolve; 1948 jac->coords[0] = NULL; 1949 jac->coords[1] = NULL; 1950 jac->coords[2] = NULL; 1951 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 1952 jac->as_print = 0; 1953 jac->as_max_iter = 1; /* used as a preconditioner */ 1954 jac->as_tol = 0.; /* used as a preconditioner */ 1955 jac->ams_cycle_type = 13; 1956 /* Smoothing options */ 1957 jac->as_relax_type = 2; 1958 jac->as_relax_times = 1; 1959 jac->as_relax_weight = 1.0; 1960 jac->as_omega = 1.0; 1961 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1962 jac->as_amg_alpha_opts[0] = 10; 1963 jac->as_amg_alpha_opts[1] = 1; 1964 jac->as_amg_alpha_opts[2] = 6; 1965 jac->as_amg_alpha_opts[3] = 6; 1966 jac->as_amg_alpha_opts[4] = 4; 1967 jac->as_amg_alpha_theta = 0.25; 1968 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1969 jac->as_amg_beta_opts[0] = 10; 1970 jac->as_amg_beta_opts[1] = 1; 1971 jac->as_amg_beta_opts[2] = 6; 1972 jac->as_amg_beta_opts[3] = 6; 1973 jac->as_amg_beta_opts[4] = 4; 1974 jac->as_amg_beta_theta = 0.25; 1975 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print)); 1976 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1977 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1978 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol)); 1979 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1980 jac->as_relax_times, 1981 jac->as_relax_weight, 1982 jac->as_omega)); 1983 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1984 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1985 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1986 jac->as_amg_alpha_theta, 1987 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1988 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1989 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1990 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1991 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1992 jac->as_amg_beta_theta, 1993 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1994 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1995 /* Zero conductivity */ 1996 jac->ams_beta_is_zero = PETSC_FALSE; 1997 jac->ams_beta_is_zero_part = PETSC_FALSE; 1998 PetscFunctionReturn(0); 1999 } 2000 ierr = PetscStrcmp("ads",jac->hypre_type,&flag);CHKERRQ(ierr); 2001 if (flag) { 2002 ierr = HYPRE_ADSCreate(&jac->hsolver); 2003 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS; 2004 pc->ops->view = PCView_HYPRE_ADS; 2005 jac->destroy = HYPRE_ADSDestroy; 2006 jac->setup = HYPRE_ADSSetup; 2007 jac->solve = HYPRE_ADSSolve; 2008 jac->coords[0] = NULL; 2009 jac->coords[1] = NULL; 2010 jac->coords[2] = NULL; 2011 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */ 2012 jac->as_print = 0; 2013 jac->as_max_iter = 1; /* used as a preconditioner */ 2014 jac->as_tol = 0.; /* used as a preconditioner */ 2015 jac->ads_cycle_type = 13; 2016 /* Smoothing options */ 2017 jac->as_relax_type = 2; 2018 jac->as_relax_times = 1; 2019 jac->as_relax_weight = 1.0; 2020 jac->as_omega = 1.0; 2021 /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2022 jac->ams_cycle_type = 14; 2023 jac->as_amg_alpha_opts[0] = 10; 2024 jac->as_amg_alpha_opts[1] = 1; 2025 jac->as_amg_alpha_opts[2] = 6; 2026 jac->as_amg_alpha_opts[3] = 6; 2027 jac->as_amg_alpha_opts[4] = 4; 2028 jac->as_amg_alpha_theta = 0.25; 2029 /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 2030 jac->as_amg_beta_opts[0] = 10; 2031 jac->as_amg_beta_opts[1] = 1; 2032 jac->as_amg_beta_opts[2] = 6; 2033 jac->as_amg_beta_opts[3] = 6; 2034 jac->as_amg_beta_opts[4] = 4; 2035 jac->as_amg_beta_theta = 0.25; 2036 PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print)); 2037 PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 2038 PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 2039 PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol)); 2040 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 2041 jac->as_relax_times, 2042 jac->as_relax_weight, 2043 jac->as_omega)); 2044 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */ 2045 jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 2046 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 2047 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 2048 jac->as_amg_alpha_theta, 2049 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 2050 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 2051 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 2052 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 2053 jac->as_amg_beta_opts[2], /* AMG relax_type */ 2054 jac->as_amg_beta_theta, 2055 jac->as_amg_beta_opts[3], /* AMG interp_type */ 2056 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 2057 PetscFunctionReturn(0); 2058 } 2059 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 2060 2061 jac->hypre_type = NULL; 2062 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name); 2063 } 2064 2065 /* 2066 It only gets here if the HYPRE type has not been set before the call to 2067 ...SetFromOptions() which actually is most of the time 2068 */ 2069 PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc) 2070 { 2071 PetscErrorCode ierr; 2072 PetscInt indx; 2073 const char *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"}; 2074 PetscBool flg; 2075 2076 PetscFunctionBegin; 2077 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr); 2078 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,ALEN(type),"boomeramg",&indx,&flg);CHKERRQ(ierr); 2079 if (flg) { 2080 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr); 2081 } else { 2082 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr); 2083 } 2084 if (pc->ops->setfromoptions) { 2085 ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr); 2086 } 2087 ierr = PetscOptionsTail();CHKERRQ(ierr); 2088 PetscFunctionReturn(0); 2089 } 2090 2091 /*@C 2092 PCHYPRESetType - Sets which hypre preconditioner you wish to use 2093 2094 Input Parameters: 2095 + pc - the preconditioner context 2096 - name - either euclid, pilut, parasails, boomeramg, ams, ads 2097 2098 Options Database Keys: 2099 -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads 2100 2101 Level: intermediate 2102 2103 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 2104 PCHYPRE 2105 2106 @*/ 2107 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 2108 { 2109 PetscErrorCode ierr; 2110 2111 PetscFunctionBegin; 2112 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2113 PetscValidCharPointer(name,2); 2114 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr); 2115 PetscFunctionReturn(0); 2116 } 2117 2118 /*@C 2119 PCHYPREGetType - Gets which hypre preconditioner you are using 2120 2121 Input Parameter: 2122 . pc - the preconditioner context 2123 2124 Output Parameter: 2125 . name - either euclid, pilut, parasails, boomeramg, ams, ads 2126 2127 Level: intermediate 2128 2129 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC, 2130 PCHYPRE 2131 2132 @*/ 2133 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 2134 { 2135 PetscErrorCode ierr; 2136 2137 PetscFunctionBegin; 2138 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 2139 PetscValidPointer(name,2); 2140 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr); 2141 PetscFunctionReturn(0); 2142 } 2143 2144 /*MC 2145 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 2146 2147 Options Database Keys: 2148 + -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads 2149 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX 2150 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 Options Database Keys: 2181 + -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal()) 2182 - -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant()) 2183 2184 Depending on the linear system you may see the same or different convergence depending on the values you use. 2185 2186 See PCPFMG for access to the hypre Struct PFMG solver 2187 2188 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 2189 PCHYPRESetType(), PCPFMG 2190 2191 M*/ 2192 2193 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 2194 { 2195 PC_HYPRE *jac; 2196 PetscErrorCode ierr; 2197 2198 PetscFunctionBegin; 2199 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 2200 2201 pc->data = jac; 2202 pc->ops->reset = PCReset_HYPRE; 2203 pc->ops->destroy = PCDestroy_HYPRE; 2204 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 2205 pc->ops->setup = PCSetUp_HYPRE; 2206 pc->ops->apply = PCApply_HYPRE; 2207 jac->comm_hypre = MPI_COMM_NULL; 2208 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 2209 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 2210 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr); 2211 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr); 2212 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr); 2213 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);CHKERRQ(ierr); 2214 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);CHKERRQ(ierr); 2215 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);CHKERRQ(ierr); 2216 #if defined(PETSC_HAVE_HYPRE_DEVICE) 2217 #if defined(HYPRE_USING_HIP) 2218 ierr = PetscHIPInitializeCheck();CHKERRQ(ierr); 2219 #endif 2220 #if defined(HYPRE_USING_CUDA) 2221 ierr = PetscCUDAInitializeCheck();CHKERRQ(ierr); 2222 #endif 2223 #endif 2224 PetscFunctionReturn(0); 2225 } 2226 2227 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 2228 2229 typedef struct { 2230 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 2231 HYPRE_StructSolver hsolver; 2232 2233 /* keep copy of PFMG options used so may view them */ 2234 PetscInt its; 2235 double tol; 2236 PetscInt relax_type; 2237 PetscInt rap_type; 2238 PetscInt num_pre_relax,num_post_relax; 2239 PetscInt max_levels; 2240 } PC_PFMG; 2241 2242 PetscErrorCode PCDestroy_PFMG(PC pc) 2243 { 2244 PetscErrorCode ierr; 2245 PC_PFMG *ex = (PC_PFMG*) pc->data; 2246 2247 PetscFunctionBegin; 2248 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 2249 ierr = MPI_Comm_free(&ex->hcomm);CHKERRMPI(ierr); 2250 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2251 PetscFunctionReturn(0); 2252 } 2253 2254 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 2255 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 2256 2257 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 2258 { 2259 PetscErrorCode ierr; 2260 PetscBool iascii; 2261 PC_PFMG *ex = (PC_PFMG*) pc->data; 2262 2263 PetscFunctionBegin; 2264 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2265 if (iascii) { 2266 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 2267 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 2268 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 2269 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 2270 ierr = PetscViewerASCIIPrintf(viewer," RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 2271 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 2272 ierr = PetscViewerASCIIPrintf(viewer," max levels %d\n",ex->max_levels);CHKERRQ(ierr); 2273 } 2274 PetscFunctionReturn(0); 2275 } 2276 2277 PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2278 { 2279 PetscErrorCode ierr; 2280 PC_PFMG *ex = (PC_PFMG*) pc->data; 2281 PetscBool flg = PETSC_FALSE; 2282 2283 PetscFunctionBegin; 2284 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr); 2285 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2286 if (flg) { 2287 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,3)); 2288 } 2289 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2290 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 2291 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); 2292 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 2293 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); 2294 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 2295 2296 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 2297 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 2298 2299 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2300 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 2301 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); 2302 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 2303 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 2304 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 2305 ierr = PetscOptionsTail();CHKERRQ(ierr); 2306 PetscFunctionReturn(0); 2307 } 2308 2309 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 2310 { 2311 PetscErrorCode ierr; 2312 PC_PFMG *ex = (PC_PFMG*) pc->data; 2313 PetscScalar *yy; 2314 const PetscScalar *xx; 2315 PetscInt ilower[3],iupper[3]; 2316 HYPRE_Int hlower[3],hupper[3]; 2317 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2318 2319 PetscFunctionBegin; 2320 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2321 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2322 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2323 iupper[0] += ilower[0] - 1; 2324 iupper[1] += ilower[1] - 1; 2325 iupper[2] += ilower[2] - 1; 2326 hlower[0] = (HYPRE_Int)ilower[0]; 2327 hlower[1] = (HYPRE_Int)ilower[1]; 2328 hlower[2] = (HYPRE_Int)ilower[2]; 2329 hupper[0] = (HYPRE_Int)iupper[0]; 2330 hupper[1] = (HYPRE_Int)iupper[1]; 2331 hupper[2] = (HYPRE_Int)iupper[2]; 2332 2333 /* copy x values over to hypre */ 2334 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 2335 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2336 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx)); 2337 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2338 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 2339 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2340 2341 /* copy solution values back to PETSc */ 2342 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2343 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy)); 2344 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2345 PetscFunctionReturn(0); 2346 } 2347 2348 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) 2349 { 2350 PC_PFMG *jac = (PC_PFMG*)pc->data; 2351 PetscErrorCode ierr; 2352 HYPRE_Int oits; 2353 2354 PetscFunctionBegin; 2355 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2356 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 2357 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 2358 2359 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 2360 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits)); 2361 *outits = oits; 2362 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2363 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2364 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 2365 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 2366 PetscFunctionReturn(0); 2367 } 2368 2369 PetscErrorCode PCSetUp_PFMG(PC pc) 2370 { 2371 PetscErrorCode ierr; 2372 PC_PFMG *ex = (PC_PFMG*) pc->data; 2373 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2374 PetscBool flg; 2375 2376 PetscFunctionBegin; 2377 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 2378 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 2379 2380 /* create the hypre solver object and set its information */ 2381 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 2382 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2383 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2384 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 2385 PetscFunctionReturn(0); 2386 } 2387 2388 /*MC 2389 PCPFMG - the hypre PFMG multigrid solver 2390 2391 Level: advanced 2392 2393 Options Database: 2394 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 2395 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2396 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2397 . -pc_pfmg_tol <tol> tolerance of PFMG 2398 . -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 2399 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 2400 2401 Notes: 2402 This is for CELL-centered descretizations 2403 2404 This must be used with the MATHYPRESTRUCT matrix type. 2405 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 2406 2407 .seealso: PCMG, MATHYPRESTRUCT 2408 M*/ 2409 2410 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 2411 { 2412 PetscErrorCode ierr; 2413 PC_PFMG *ex; 2414 2415 PetscFunctionBegin; 2416 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2417 pc->data = ex; 2418 2419 ex->its = 1; 2420 ex->tol = 1.e-8; 2421 ex->relax_type = 1; 2422 ex->rap_type = 0; 2423 ex->num_pre_relax = 1; 2424 ex->num_post_relax = 1; 2425 ex->max_levels = 0; 2426 2427 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 2428 pc->ops->view = PCView_PFMG; 2429 pc->ops->destroy = PCDestroy_PFMG; 2430 pc->ops->apply = PCApply_PFMG; 2431 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 2432 pc->ops->setup = PCSetUp_PFMG; 2433 2434 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRMPI(ierr); 2435 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2436 PetscFunctionReturn(0); 2437 } 2438 2439 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 2440 2441 /* we know we are working with a HYPRE_SStructMatrix */ 2442 typedef struct { 2443 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 2444 HYPRE_SStructSolver ss_solver; 2445 2446 /* keep copy of SYSPFMG options used so may view them */ 2447 PetscInt its; 2448 double tol; 2449 PetscInt relax_type; 2450 PetscInt num_pre_relax,num_post_relax; 2451 } PC_SysPFMG; 2452 2453 PetscErrorCode PCDestroy_SysPFMG(PC pc) 2454 { 2455 PetscErrorCode ierr; 2456 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2457 2458 PetscFunctionBegin; 2459 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2460 ierr = MPI_Comm_free(&ex->hcomm);CHKERRMPI(ierr); 2461 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2462 PetscFunctionReturn(0); 2463 } 2464 2465 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 2466 2467 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 2468 { 2469 PetscErrorCode ierr; 2470 PetscBool iascii; 2471 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2472 2473 PetscFunctionBegin; 2474 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2475 if (iascii) { 2476 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 2477 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 2478 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 2479 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 2480 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 2481 } 2482 PetscFunctionReturn(0); 2483 } 2484 2485 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2486 { 2487 PetscErrorCode ierr; 2488 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2489 PetscBool flg = PETSC_FALSE; 2490 2491 PetscFunctionBegin; 2492 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr); 2493 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2494 if (flg) { 2495 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,3)); 2496 } 2497 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2498 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 2499 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); 2500 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 2501 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); 2502 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 2503 2504 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2505 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 2506 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); 2507 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 2508 ierr = PetscOptionsTail();CHKERRQ(ierr); 2509 PetscFunctionReturn(0); 2510 } 2511 2512 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 2513 { 2514 PetscErrorCode ierr; 2515 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2516 PetscScalar *yy; 2517 const PetscScalar *xx; 2518 PetscInt ilower[3],iupper[3]; 2519 HYPRE_Int hlower[3],hupper[3]; 2520 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2521 PetscInt ordering= mx->dofs_order; 2522 PetscInt nvars = mx->nvars; 2523 PetscInt part = 0; 2524 PetscInt size; 2525 PetscInt i; 2526 2527 PetscFunctionBegin; 2528 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2529 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2530 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */ 2531 iupper[0] += ilower[0] - 1; 2532 iupper[1] += ilower[1] - 1; 2533 iupper[2] += ilower[2] - 1; 2534 hlower[0] = (HYPRE_Int)ilower[0]; 2535 hlower[1] = (HYPRE_Int)ilower[1]; 2536 hlower[2] = (HYPRE_Int)ilower[2]; 2537 hupper[0] = (HYPRE_Int)iupper[0]; 2538 hupper[1] = (HYPRE_Int)iupper[1]; 2539 hupper[2] = (HYPRE_Int)iupper[2]; 2540 2541 size = 1; 2542 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 2543 2544 /* copy x values over to hypre for variable ordering */ 2545 if (ordering) { 2546 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2547 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2548 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i)))); 2549 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2550 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2551 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 2552 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2553 2554 /* copy solution values back to PETSc */ 2555 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2556 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i)))); 2557 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2558 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2559 PetscScalar *z; 2560 PetscInt j, k; 2561 2562 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 2563 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2564 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2565 2566 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2567 for (i= 0; i< size; i++) { 2568 k= i*nvars; 2569 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 2570 } 2571 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)))); 2572 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2573 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2574 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2575 2576 /* copy solution values back to PETSc */ 2577 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2578 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)))); 2579 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2580 for (i= 0; i< size; i++) { 2581 k= i*nvars; 2582 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 2583 } 2584 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2585 ierr = PetscFree(z);CHKERRQ(ierr); 2586 } 2587 PetscFunctionReturn(0); 2588 } 2589 2590 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) 2591 { 2592 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 2593 PetscErrorCode ierr; 2594 HYPRE_Int oits; 2595 2596 PetscFunctionBegin; 2597 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2598 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 2599 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 2600 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 2601 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits)); 2602 *outits = oits; 2603 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2604 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2605 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 2606 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 2607 PetscFunctionReturn(0); 2608 } 2609 2610 PetscErrorCode PCSetUp_SysPFMG(PC pc) 2611 { 2612 PetscErrorCode ierr; 2613 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2614 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2615 PetscBool flg; 2616 2617 PetscFunctionBegin; 2618 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 2619 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 2620 2621 /* create the hypre sstruct solver object and set its information */ 2622 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2623 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2624 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 2625 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2626 PetscFunctionReturn(0); 2627 } 2628 2629 /*MC 2630 PCSysPFMG - the hypre SysPFMG multigrid solver 2631 2632 Level: advanced 2633 2634 Options Database: 2635 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 2636 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2637 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2638 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 2639 - -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 2640 2641 Notes: 2642 This is for CELL-centered descretizations 2643 2644 This must be used with the MATHYPRESSTRUCT matrix type. 2645 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 2646 Also, only cell-centered variables. 2647 2648 .seealso: PCMG, MATHYPRESSTRUCT 2649 M*/ 2650 2651 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2652 { 2653 PetscErrorCode ierr; 2654 PC_SysPFMG *ex; 2655 2656 PetscFunctionBegin; 2657 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2658 pc->data = ex; 2659 2660 ex->its = 1; 2661 ex->tol = 1.e-8; 2662 ex->relax_type = 1; 2663 ex->num_pre_relax = 1; 2664 ex->num_post_relax = 1; 2665 2666 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2667 pc->ops->view = PCView_SysPFMG; 2668 pc->ops->destroy = PCDestroy_SysPFMG; 2669 pc->ops->apply = PCApply_SysPFMG; 2670 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2671 pc->ops->setup = PCSetUp_SysPFMG; 2672 2673 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRMPI(ierr); 2674 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2675 PetscFunctionReturn(0); 2676 } 2677