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