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