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