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