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