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