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