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