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