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->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */ 1370 jac->relaxtype[2] = 9; /*G.E. */ 1371 jac->relaxweight = 1.0; 1372 jac->outerrelaxweight = 1.0; 1373 jac->relaxorder = 1; 1374 jac->interptype = 0; 1375 jac->agg_nl = 0; 1376 jac->pmax = 0; 1377 jac->truncfactor = 0.0; 1378 jac->agg_num_paths = 1; 1379 1380 jac->nodal_coarsen = 0; 1381 jac->nodal_relax = PETSC_FALSE; 1382 jac->nodal_relax_levels = 1; 1383 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 1384 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 1385 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 1386 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 1387 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 1388 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 1389 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 1390 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 1391 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 1392 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 1393 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 1394 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 1395 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 1396 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 1397 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/ 1398 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */ 1399 PetscFunctionReturn(0); 1400 } 1401 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr); 1402 if (flag) { 1403 ierr = HYPRE_AMSCreate(&jac->hsolver); 1404 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 1405 pc->ops->view = PCView_HYPRE_AMS; 1406 jac->destroy = HYPRE_AMSDestroy; 1407 jac->setup = HYPRE_AMSSetup; 1408 jac->solve = HYPRE_AMSSolve; 1409 jac->setdgrad = HYPRE_AMSSetDiscreteGradient; 1410 jac->setcoord = HYPRE_AMSSetCoordinateVectors; 1411 jac->setdim = HYPRE_AMSSetDimension; 1412 jac->coords[0] = NULL; 1413 jac->coords[1] = NULL; 1414 jac->coords[2] = NULL; 1415 jac->G = NULL; 1416 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 1417 jac->as_print = 0; 1418 jac->as_max_iter = 1; /* used as a preconditioner */ 1419 jac->as_tol = 0.; /* used as a preconditioner */ 1420 jac->ams_cycle_type = 13; 1421 /* Smoothing options */ 1422 jac->as_relax_type = 2; 1423 jac->as_relax_times = 1; 1424 jac->as_relax_weight = 1.0; 1425 jac->as_omega = 1.0; 1426 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1427 jac->as_amg_alpha_opts[0] = 10; 1428 jac->as_amg_alpha_opts[1] = 1; 1429 jac->as_amg_alpha_opts[2] = 8; 1430 jac->as_amg_alpha_opts[3] = 6; 1431 jac->as_amg_alpha_opts[4] = 4; 1432 jac->as_amg_alpha_theta = 0.25; 1433 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1434 jac->ams_beta_is_zero = PETSC_FALSE; 1435 jac->as_amg_beta_opts[0] = 10; 1436 jac->as_amg_beta_opts[1] = 1; 1437 jac->as_amg_beta_opts[2] = 8; 1438 jac->as_amg_beta_opts[3] = 6; 1439 jac->as_amg_beta_opts[4] = 4; 1440 jac->as_amg_beta_theta = 0.25; 1441 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print)); 1442 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1443 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1444 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol)); 1445 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1446 jac->as_relax_times, 1447 jac->as_relax_weight, 1448 jac->as_omega)); 1449 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1450 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1451 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1452 jac->as_amg_alpha_theta, 1453 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1454 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1455 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1456 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1457 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1458 jac->as_amg_beta_theta, 1459 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1460 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1461 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr); 1462 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr); 1463 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);CHKERRQ(ierr); 1464 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1465 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1466 PetscFunctionReturn(0); 1467 } 1468 ierr = PetscStrcmp("ads",jac->hypre_type,&flag);CHKERRQ(ierr); 1469 if (flag) { 1470 ierr = HYPRE_ADSCreate(&jac->hsolver); 1471 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS; 1472 pc->ops->view = PCView_HYPRE_ADS; 1473 jac->destroy = HYPRE_ADSDestroy; 1474 jac->setup = HYPRE_ADSSetup; 1475 jac->solve = HYPRE_ADSSolve; 1476 jac->setdgrad = HYPRE_ADSSetDiscreteGradient; 1477 jac->setdcurl = HYPRE_ADSSetDiscreteCurl; 1478 jac->setcoord = HYPRE_ADSSetCoordinateVectors; 1479 jac->coords[0] = NULL; 1480 jac->coords[1] = NULL; 1481 jac->coords[2] = NULL; 1482 jac->G = NULL; 1483 jac->C = NULL; 1484 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */ 1485 jac->as_print = 0; 1486 jac->as_max_iter = 1; /* used as a preconditioner */ 1487 jac->as_tol = 0.; /* used as a preconditioner */ 1488 jac->ads_cycle_type = 13; 1489 /* Smoothing options */ 1490 jac->as_relax_type = 2; 1491 jac->as_relax_times = 1; 1492 jac->as_relax_weight = 1.0; 1493 jac->as_omega = 1.0; 1494 /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1495 jac->ams_cycle_type = 14; 1496 jac->as_amg_alpha_opts[0] = 10; 1497 jac->as_amg_alpha_opts[1] = 1; 1498 jac->as_amg_alpha_opts[2] = 6; 1499 jac->as_amg_alpha_opts[3] = 6; 1500 jac->as_amg_alpha_opts[4] = 4; 1501 jac->as_amg_alpha_theta = 0.25; 1502 /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1503 jac->as_amg_beta_opts[0] = 10; 1504 jac->as_amg_beta_opts[1] = 1; 1505 jac->as_amg_beta_opts[2] = 6; 1506 jac->as_amg_beta_opts[3] = 6; 1507 jac->as_amg_beta_opts[4] = 4; 1508 jac->as_amg_beta_theta = 0.25; 1509 PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print)); 1510 PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1511 PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1512 PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol)); 1513 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1514 jac->as_relax_times, 1515 jac->as_relax_weight, 1516 jac->as_omega)); 1517 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */ 1518 jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1519 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1520 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1521 jac->as_amg_alpha_theta, 1522 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1523 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1524 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1525 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1526 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1527 jac->as_amg_beta_theta, 1528 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1529 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1530 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr); 1531 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr); 1532 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr); 1533 PetscFunctionReturn(0); 1534 } 1535 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 1536 1537 jac->hypre_type = NULL; 1538 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 /* 1543 It only gets here if the HYPRE type has not been set before the call to 1544 ...SetFromOptions() which actually is most of the time 1545 */ 1546 #undef __FUNCT__ 1547 #define __FUNCT__ "PCSetFromOptions_HYPRE" 1548 static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc) 1549 { 1550 PetscErrorCode ierr; 1551 PetscInt indx; 1552 const char *type[] = {"pilut","parasails","boomeramg","ams","ads"}; 1553 PetscBool flg; 1554 1555 PetscFunctionBegin; 1556 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr); 1557 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);CHKERRQ(ierr); 1558 if (flg) { 1559 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr); 1560 } else { 1561 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr); 1562 } 1563 if (pc->ops->setfromoptions) { 1564 ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr); 1565 } 1566 ierr = PetscOptionsTail();CHKERRQ(ierr); 1567 PetscFunctionReturn(0); 1568 } 1569 1570 #undef __FUNCT__ 1571 #define __FUNCT__ "PCHYPRESetType" 1572 /*@C 1573 PCHYPRESetType - Sets which hypre preconditioner you wish to use 1574 1575 Input Parameters: 1576 + pc - the preconditioner context 1577 - name - either pilut, parasails, boomeramg, ams, ads 1578 1579 Options Database Keys: 1580 -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads 1581 1582 Level: intermediate 1583 1584 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1585 PCHYPRE 1586 1587 @*/ 1588 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 1589 { 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1594 PetscValidCharPointer(name,2); 1595 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr); 1596 PetscFunctionReturn(0); 1597 } 1598 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "PCHYPREGetType" 1601 /*@C 1602 PCHYPREGetType - Gets which hypre preconditioner you are using 1603 1604 Input Parameter: 1605 . pc - the preconditioner context 1606 1607 Output Parameter: 1608 . name - either pilut, parasails, boomeramg, ams, ads 1609 1610 Level: intermediate 1611 1612 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC, 1613 PCHYPRE 1614 1615 @*/ 1616 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 1617 { 1618 PetscErrorCode ierr; 1619 1620 PetscFunctionBegin; 1621 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1622 PetscValidPointer(name,2); 1623 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr); 1624 PetscFunctionReturn(0); 1625 } 1626 1627 /*MC 1628 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 1629 1630 Options Database Keys: 1631 + -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads 1632 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX 1633 preconditioner 1634 1635 Level: intermediate 1636 1637 Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()), 1638 the many hypre options can ONLY be set via the options database (e.g. the command line 1639 or with PetscOptionsSetValue(), there are no functions to set them) 1640 1641 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations 1642 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if 1643 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner 1644 (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of 1645 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations 1646 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 1647 then AT MOST twenty V-cycles of boomeramg will be called. 1648 1649 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation 1650 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 1651 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi. 1652 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly 1653 and use -ksp_max_it to control the number of V-cycles. 1654 (see the PETSc FAQ.html at the PETSc website under the Documentation tab). 1655 1656 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option 1657 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L. 1658 1659 See PCPFMG for access to the hypre Struct PFMG solver 1660 1661 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1662 PCHYPRESetType(), PCPFMG 1663 1664 M*/ 1665 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "PCCreate_HYPRE" 1668 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 1669 { 1670 PC_HYPRE *jac; 1671 PetscErrorCode ierr; 1672 1673 PetscFunctionBegin; 1674 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 1675 1676 pc->data = jac; 1677 pc->ops->destroy = PCDestroy_HYPRE; 1678 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 1679 pc->ops->setup = PCSetUp_HYPRE; 1680 pc->ops->apply = PCApply_HYPRE; 1681 jac->comm_hypre = MPI_COMM_NULL; 1682 jac->hypre_type = NULL; 1683 jac->coords[0] = NULL; 1684 jac->coords[1] = NULL; 1685 jac->coords[2] = NULL; 1686 jac->constants[0] = NULL; 1687 jac->constants[1] = NULL; 1688 jac->constants[2] = NULL; 1689 jac->G = NULL; 1690 jac->C = NULL; 1691 jac->alpha_Poisson = NULL; 1692 jac->beta_Poisson = NULL; 1693 jac->setdim = NULL; 1694 /* duplicate communicator for hypre */ 1695 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr); 1696 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 1697 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 1698 PetscFunctionReturn(0); 1699 } 1700 1701 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 1702 1703 /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */ 1704 #include <petsc/private/matimpl.h> 1705 1706 typedef struct { 1707 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 1708 HYPRE_StructSolver hsolver; 1709 1710 /* keep copy of PFMG options used so may view them */ 1711 PetscInt its; 1712 double tol; 1713 PetscInt relax_type; 1714 PetscInt rap_type; 1715 PetscInt num_pre_relax,num_post_relax; 1716 PetscInt max_levels; 1717 } PC_PFMG; 1718 1719 #undef __FUNCT__ 1720 #define __FUNCT__ "PCDestroy_PFMG" 1721 PetscErrorCode PCDestroy_PFMG(PC pc) 1722 { 1723 PetscErrorCode ierr; 1724 PC_PFMG *ex = (PC_PFMG*) pc->data; 1725 1726 PetscFunctionBegin; 1727 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1728 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1729 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1730 PetscFunctionReturn(0); 1731 } 1732 1733 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 1734 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 1735 1736 #undef __FUNCT__ 1737 #define __FUNCT__ "PCView_PFMG" 1738 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 1739 { 1740 PetscErrorCode ierr; 1741 PetscBool iascii; 1742 PC_PFMG *ex = (PC_PFMG*) pc->data; 1743 1744 PetscFunctionBegin; 1745 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1746 if (iascii) { 1747 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 1748 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1749 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1750 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1751 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 1752 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1753 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);CHKERRQ(ierr); 1754 } 1755 PetscFunctionReturn(0); 1756 } 1757 1758 1759 #undef __FUNCT__ 1760 #define __FUNCT__ "PCSetFromOptions_PFMG" 1761 PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc) 1762 { 1763 PetscErrorCode ierr; 1764 PC_PFMG *ex = (PC_PFMG*) pc->data; 1765 PetscBool flg = PETSC_FALSE; 1766 1767 PetscFunctionBegin; 1768 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr); 1769 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1770 if (flg) { 1771 int level=3; 1772 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level)); 1773 } 1774 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1775 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 1776 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); 1777 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 1778 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); 1779 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 1780 1781 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 1782 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 1783 1784 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1785 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 1786 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); 1787 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 1788 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 1789 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 1790 ierr = PetscOptionsTail();CHKERRQ(ierr); 1791 PetscFunctionReturn(0); 1792 } 1793 1794 #undef __FUNCT__ 1795 #define __FUNCT__ "PCApply_PFMG" 1796 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 1797 { 1798 PetscErrorCode ierr; 1799 PC_PFMG *ex = (PC_PFMG*) pc->data; 1800 PetscScalar *yy; 1801 const PetscScalar *xx; 1802 PetscInt ilower[3],iupper[3]; 1803 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1804 1805 PetscFunctionBegin; 1806 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1807 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1808 iupper[0] += ilower[0] - 1; 1809 iupper[1] += ilower[1] - 1; 1810 iupper[2] += ilower[2] - 1; 1811 1812 /* copy x values over to hypre */ 1813 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 1814 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1815 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx)); 1816 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1817 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 1818 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1819 1820 /* copy solution values back to PETSc */ 1821 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1822 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 1823 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1824 PetscFunctionReturn(0); 1825 } 1826 1827 #undef __FUNCT__ 1828 #define __FUNCT__ "PCApplyRichardson_PFMG" 1829 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) 1830 { 1831 PC_PFMG *jac = (PC_PFMG*)pc->data; 1832 PetscErrorCode ierr; 1833 PetscInt oits; 1834 1835 PetscFunctionBegin; 1836 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1837 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 1838 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 1839 1840 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 1841 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 1842 *outits = oits; 1843 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1844 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1845 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 1846 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 1847 PetscFunctionReturn(0); 1848 } 1849 1850 1851 #undef __FUNCT__ 1852 #define __FUNCT__ "PCSetUp_PFMG" 1853 PetscErrorCode PCSetUp_PFMG(PC pc) 1854 { 1855 PetscErrorCode ierr; 1856 PC_PFMG *ex = (PC_PFMG*) pc->data; 1857 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1858 PetscBool flg; 1859 1860 PetscFunctionBegin; 1861 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 1862 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 1863 1864 /* create the hypre solver object and set its information */ 1865 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1866 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1867 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1868 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 1869 PetscFunctionReturn(0); 1870 } 1871 1872 1873 /*MC 1874 PCPFMG - the hypre PFMG multigrid solver 1875 1876 Level: advanced 1877 1878 Options Database: 1879 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 1880 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 1881 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 1882 . -pc_pfmg_tol <tol> tolerance of PFMG 1883 . -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 1884 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 1885 1886 Notes: This is for CELL-centered descretizations 1887 1888 This must be used with the MATHYPRESTRUCT matrix type. 1889 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 1890 1891 .seealso: PCMG, MATHYPRESTRUCT 1892 M*/ 1893 1894 #undef __FUNCT__ 1895 #define __FUNCT__ "PCCreate_PFMG" 1896 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 1897 { 1898 PetscErrorCode ierr; 1899 PC_PFMG *ex; 1900 1901 PetscFunctionBegin; 1902 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 1903 pc->data = ex; 1904 1905 ex->its = 1; 1906 ex->tol = 1.e-8; 1907 ex->relax_type = 1; 1908 ex->rap_type = 0; 1909 ex->num_pre_relax = 1; 1910 ex->num_post_relax = 1; 1911 ex->max_levels = 0; 1912 1913 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 1914 pc->ops->view = PCView_PFMG; 1915 pc->ops->destroy = PCDestroy_PFMG; 1916 pc->ops->apply = PCApply_PFMG; 1917 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 1918 pc->ops->setup = PCSetUp_PFMG; 1919 1920 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 1921 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1922 PetscFunctionReturn(0); 1923 } 1924 1925 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 1926 1927 /* we know we are working with a HYPRE_SStructMatrix */ 1928 typedef struct { 1929 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 1930 HYPRE_SStructSolver ss_solver; 1931 1932 /* keep copy of SYSPFMG options used so may view them */ 1933 PetscInt its; 1934 double tol; 1935 PetscInt relax_type; 1936 PetscInt num_pre_relax,num_post_relax; 1937 } PC_SysPFMG; 1938 1939 #undef __FUNCT__ 1940 #define __FUNCT__ "PCDestroy_SysPFMG" 1941 PetscErrorCode PCDestroy_SysPFMG(PC pc) 1942 { 1943 PetscErrorCode ierr; 1944 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1945 1946 PetscFunctionBegin; 1947 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 1948 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1949 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 1954 1955 #undef __FUNCT__ 1956 #define __FUNCT__ "PCView_SysPFMG" 1957 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 1958 { 1959 PetscErrorCode ierr; 1960 PetscBool iascii; 1961 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1962 1963 PetscFunctionBegin; 1964 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1965 if (iascii) { 1966 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 1967 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1968 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1969 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1970 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1971 } 1972 PetscFunctionReturn(0); 1973 } 1974 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "PCSetFromOptions_SysPFMG" 1978 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc) 1979 { 1980 PetscErrorCode ierr; 1981 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1982 PetscBool flg = PETSC_FALSE; 1983 1984 PetscFunctionBegin; 1985 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr); 1986 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1987 if (flg) { 1988 int level=3; 1989 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level)); 1990 } 1991 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1992 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 1993 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); 1994 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 1995 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); 1996 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 1997 1998 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1999 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 2000 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); 2001 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 2002 ierr = PetscOptionsTail();CHKERRQ(ierr); 2003 PetscFunctionReturn(0); 2004 } 2005 2006 #undef __FUNCT__ 2007 #define __FUNCT__ "PCApply_SysPFMG" 2008 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 2009 { 2010 PetscErrorCode ierr; 2011 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2012 PetscScalar *yy; 2013 const PetscScalar *xx; 2014 PetscInt ilower[3],iupper[3]; 2015 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2016 PetscInt ordering= mx->dofs_order; 2017 PetscInt nvars = mx->nvars; 2018 PetscInt part = 0; 2019 PetscInt size; 2020 PetscInt i; 2021 2022 PetscFunctionBegin; 2023 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2024 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2025 iupper[0] += ilower[0] - 1; 2026 iupper[1] += ilower[1] - 1; 2027 iupper[2] += ilower[2] - 1; 2028 2029 size = 1; 2030 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 2031 2032 /* copy x values over to hypre for variable ordering */ 2033 if (ordering) { 2034 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2035 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2036 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i))); 2037 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2038 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2039 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 2040 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2041 2042 /* copy solution values back to PETSc */ 2043 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2044 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 2045 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2046 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2047 PetscScalar *z; 2048 PetscInt j, k; 2049 2050 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 2051 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2052 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2053 2054 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2055 for (i= 0; i< size; i++) { 2056 k= i*nvars; 2057 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 2058 } 2059 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2060 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2061 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2062 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2063 2064 /* copy solution values back to PETSc */ 2065 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2066 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2067 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2068 for (i= 0; i< size; i++) { 2069 k= i*nvars; 2070 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 2071 } 2072 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2073 ierr = PetscFree(z);CHKERRQ(ierr); 2074 } 2075 PetscFunctionReturn(0); 2076 } 2077 2078 #undef __FUNCT__ 2079 #define __FUNCT__ "PCApplyRichardson_SysPFMG" 2080 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) 2081 { 2082 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 2083 PetscErrorCode ierr; 2084 PetscInt oits; 2085 2086 PetscFunctionBegin; 2087 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2088 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 2089 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 2090 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 2091 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits)); 2092 *outits = oits; 2093 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2094 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2095 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 2096 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 2097 PetscFunctionReturn(0); 2098 } 2099 2100 2101 #undef __FUNCT__ 2102 #define __FUNCT__ "PCSetUp_SysPFMG" 2103 PetscErrorCode PCSetUp_SysPFMG(PC pc) 2104 { 2105 PetscErrorCode ierr; 2106 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2107 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2108 PetscBool flg; 2109 2110 PetscFunctionBegin; 2111 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 2112 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 2113 2114 /* create the hypre sstruct solver object and set its information */ 2115 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2116 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2117 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 2118 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2119 PetscFunctionReturn(0); 2120 } 2121 2122 2123 /*MC 2124 PCSysPFMG - the hypre SysPFMG multigrid solver 2125 2126 Level: advanced 2127 2128 Options Database: 2129 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 2130 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2131 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2132 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 2133 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 2134 2135 Notes: This is for CELL-centered descretizations 2136 2137 This must be used with the MATHYPRESSTRUCT matrix type. 2138 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 2139 Also, only cell-centered variables. 2140 2141 .seealso: PCMG, MATHYPRESSTRUCT 2142 M*/ 2143 2144 #undef __FUNCT__ 2145 #define __FUNCT__ "PCCreate_SysPFMG" 2146 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2147 { 2148 PetscErrorCode ierr; 2149 PC_SysPFMG *ex; 2150 2151 PetscFunctionBegin; 2152 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2153 pc->data = ex; 2154 2155 ex->its = 1; 2156 ex->tol = 1.e-8; 2157 ex->relax_type = 1; 2158 ex->num_pre_relax = 1; 2159 ex->num_post_relax = 1; 2160 2161 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2162 pc->ops->view = PCView_SysPFMG; 2163 pc->ops->destroy = PCDestroy_SysPFMG; 2164 pc->ops->apply = PCApply_SysPFMG; 2165 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2166 pc->ops->setup = PCSetUp_SysPFMG; 2167 2168 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 2169 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2170 PetscFunctionReturn(0); 2171 } 2172