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