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