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