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