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