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