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 PetscErrorCode ierr; 1033 1034 PetscFunctionBegin; 1035 /* throw away any vector if already set */ 1036 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); 1037 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); 1038 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); 1039 jac->constants[0] = NULL; 1040 jac->constants[1] = NULL; 1041 jac->constants[2] = NULL; 1042 ierr = VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);CHKERRQ(ierr); 1043 ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr); 1044 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz))); 1045 ierr = VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);CHKERRQ(ierr); 1046 ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr); 1047 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz))); 1048 if (zzo) { 1049 ierr = VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);CHKERRQ(ierr); 1050 ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr); 1051 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo))); 1052 } 1053 PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo)); 1054 PetscFunctionReturn(0); 1055 } 1056 1057 #undef __FUNCT__ 1058 #define __FUNCT__ "PCHYPRESetEdgeConstantVectors" 1059 /*@ 1060 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis 1061 1062 Collective on PC 1063 1064 Input Parameters: 1065 + pc - the preconditioning context 1066 - ozz - vector representing (1,0,0) (or (1,0) in 2D) 1067 - zoz - vector representing (0,1,0) (or (0,1) in 2D) 1068 - zzo - vector representing (0,0,1) (use NULL in 2D) 1069 1070 Level: intermediate 1071 1072 Notes: 1073 1074 .seealso: 1075 @*/ 1076 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) 1077 { 1078 PetscErrorCode ierr; 1079 1080 PetscFunctionBegin; 1081 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1082 PetscValidHeaderSpecific(ozz,VEC_CLASSID,2); 1083 PetscValidHeaderSpecific(zoz,VEC_CLASSID,3); 1084 if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4); 1085 PetscCheckSameComm(pc,1,ozz,2); 1086 PetscCheckSameComm(pc,1,zoz,3); 1087 if (zzo) PetscCheckSameComm(pc,1,zzo,4); 1088 ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr); 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "PCSetCoordinates_HYPRE_AMS" 1094 static PetscErrorCode PCSetCoordinates_HYPRE_AMS(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 1095 { 1096 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1097 Vec tv; 1098 HYPRE_ParVector par_coords[3]; 1099 PetscInt i; 1100 PetscErrorCode ierr; 1101 1102 PetscFunctionBegin; 1103 /* throw away any coordinate vector if already set */ 1104 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); 1105 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); 1106 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); 1107 /* set problem's dimension */ 1108 PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim)); 1109 /* compute IJ vector for coordinates */ 1110 ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr); 1111 ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr); 1112 ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr); 1113 for (i=0;i<dim;i++) { 1114 PetscScalar *array; 1115 PetscInt j; 1116 1117 ierr = VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);CHKERRQ(ierr); 1118 ierr = VecGetArray(tv,&array);CHKERRQ(ierr); 1119 for (j=0;j<nloc;j++) { 1120 array[j] = coords[j*dim+i]; 1121 } 1122 PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array)); 1123 PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i])); 1124 ierr = VecRestoreArray(tv,&array);CHKERRQ(ierr); 1125 } 1126 ierr = VecDestroy(&tv);CHKERRQ(ierr); 1127 /* pass parCSR vectors to AMS solver */ 1128 par_coords[0] = NULL; 1129 par_coords[1] = NULL; 1130 par_coords[2] = NULL; 1131 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0]))); 1132 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1]))); 1133 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2]))); 1134 PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,par_coords[0],par_coords[1],par_coords[2])); 1135 PetscFunctionReturn(0); 1136 } 1137 1138 /* ---------------------------------------------------------------------------------*/ 1139 1140 #undef __FUNCT__ 1141 #define __FUNCT__ "PCHYPREGetType_HYPRE" 1142 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[]) 1143 { 1144 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1145 1146 PetscFunctionBegin; 1147 *name = jac->hypre_type; 1148 PetscFunctionReturn(0); 1149 } 1150 1151 #undef __FUNCT__ 1152 #define __FUNCT__ "PCHYPRESetType_HYPRE" 1153 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[]) 1154 { 1155 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1156 PetscErrorCode ierr; 1157 PetscBool flag; 1158 1159 PetscFunctionBegin; 1160 if (jac->hypre_type) { 1161 ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr); 1162 if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set"); 1163 PetscFunctionReturn(0); 1164 } else { 1165 ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr); 1166 } 1167 1168 jac->maxiter = PETSC_DEFAULT; 1169 jac->tol = PETSC_DEFAULT; 1170 jac->printstatistics = PetscLogPrintInfo; 1171 1172 ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr); 1173 if (flag) { 1174 PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver)); 1175 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; 1176 pc->ops->view = PCView_HYPRE_Pilut; 1177 jac->destroy = HYPRE_ParCSRPilutDestroy; 1178 jac->setup = HYPRE_ParCSRPilutSetup; 1179 jac->solve = HYPRE_ParCSRPilutSolve; 1180 jac->factorrowsize = PETSC_DEFAULT; 1181 PetscFunctionReturn(0); 1182 } 1183 ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr); 1184 if (flag) { 1185 PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver)); 1186 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; 1187 pc->ops->view = PCView_HYPRE_ParaSails; 1188 jac->destroy = HYPRE_ParaSailsDestroy; 1189 jac->setup = HYPRE_ParaSailsSetup; 1190 jac->solve = HYPRE_ParaSailsSolve; 1191 /* initialize */ 1192 jac->nlevels = 1; 1193 jac->threshhold = .1; 1194 jac->filter = .1; 1195 jac->loadbal = 0; 1196 if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE; 1197 else jac->logging = (int) PETSC_FALSE; 1198 1199 jac->ruse = (int) PETSC_FALSE; 1200 jac->symt = 0; 1201 PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels)); 1202 PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 1203 PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 1204 PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 1205 PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 1206 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 1207 PetscFunctionReturn(0); 1208 } 1209 ierr = PetscStrcmp("euclid",jac->hypre_type,&flag);CHKERRQ(ierr); 1210 if (flag) { 1211 ierr = HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver); 1212 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid; 1213 pc->ops->view = PCView_HYPRE_Euclid; 1214 jac->destroy = HYPRE_EuclidDestroy; 1215 jac->setup = HYPRE_EuclidSetup; 1216 jac->solve = HYPRE_EuclidSolve; 1217 /* initialization */ 1218 jac->bjilu = PETSC_FALSE; 1219 jac->levels = 1; 1220 PetscFunctionReturn(0); 1221 } 1222 ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr); 1223 if (flag) { 1224 ierr = HYPRE_BoomerAMGCreate(&jac->hsolver); 1225 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; 1226 pc->ops->view = PCView_HYPRE_BoomerAMG; 1227 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; 1228 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; 1229 jac->destroy = HYPRE_BoomerAMGDestroy; 1230 jac->setup = HYPRE_BoomerAMGSetup; 1231 jac->solve = HYPRE_BoomerAMGSolve; 1232 jac->applyrichardson = PETSC_FALSE; 1233 /* these defaults match the hypre defaults */ 1234 jac->cycletype = 1; 1235 jac->maxlevels = 25; 1236 jac->maxiter = 1; 1237 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ 1238 jac->truncfactor = 0.0; 1239 jac->strongthreshold = .25; 1240 jac->maxrowsum = .9; 1241 jac->coarsentype = 6; 1242 jac->measuretype = 0; 1243 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; 1244 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */ 1245 jac->relaxtype[2] = 9; /*G.E. */ 1246 jac->relaxweight = 1.0; 1247 jac->outerrelaxweight = 1.0; 1248 jac->relaxorder = 1; 1249 jac->interptype = 0; 1250 jac->agg_nl = 0; 1251 jac->pmax = 0; 1252 jac->truncfactor = 0.0; 1253 jac->agg_num_paths = 1; 1254 1255 jac->nodal_coarsen = 0; 1256 jac->nodal_relax = PETSC_FALSE; 1257 jac->nodal_relax_levels = 1; 1258 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 1259 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 1260 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 1261 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 1262 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 1263 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 1264 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 1265 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 1266 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 1267 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 1268 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 1269 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 1270 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 1271 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 1272 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/ 1273 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */ 1274 PetscFunctionReturn(0); 1275 } 1276 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr); 1277 if (flag) { 1278 ierr = HYPRE_AMSCreate(&jac->hsolver); 1279 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 1280 pc->ops->view = PCView_HYPRE_AMS; 1281 jac->destroy = HYPRE_AMSDestroy; 1282 jac->setup = HYPRE_AMSSetup; 1283 jac->solve = HYPRE_AMSSolve; 1284 jac->coords[0] = NULL; 1285 jac->coords[1] = NULL; 1286 jac->coords[2] = NULL; 1287 jac->G = NULL; 1288 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 1289 jac->ams_print = 0; 1290 jac->ams_max_iter = 1; /* used as a preconditioner */ 1291 jac->ams_cycle_type = 13; 1292 jac->ams_tol = 0.; /* used as a preconditioner */ 1293 /* Smoothing options */ 1294 jac->ams_relax_type = 2; 1295 jac->ams_relax_times = 1; 1296 jac->ams_relax_weight = 1.0; 1297 jac->ams_omega = 1.0; 1298 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1299 jac->ams_amg_alpha_opts[0] = 10; 1300 jac->ams_amg_alpha_opts[1] = 1; 1301 jac->ams_amg_alpha_opts[2] = 8; 1302 jac->ams_amg_alpha_opts[3] = 6; 1303 jac->ams_amg_alpha_opts[4] = 4; 1304 jac->ams_amg_alpha_theta = 0.25; 1305 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1306 jac->ams_beta_is_zero = PETSC_FALSE; 1307 jac->ams_amg_beta_opts[0] = 10; 1308 jac->ams_amg_beta_opts[1] = 1; 1309 jac->ams_amg_beta_opts[2] = 8; 1310 jac->ams_amg_beta_opts[3] = 6; 1311 jac->ams_amg_beta_opts[4] = 4; 1312 jac->ams_amg_beta_theta = 0.25; 1313 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->ams_print)); 1314 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->ams_max_iter)); 1315 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1316 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->ams_tol)); 1317 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->ams_relax_type, 1318 jac->ams_relax_times, 1319 jac->ams_relax_weight, 1320 jac->ams_omega)); 1321 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->ams_amg_alpha_opts[0], /* AMG coarsen type */ 1322 jac->ams_amg_alpha_opts[1], /* AMG agg_levels */ 1323 jac->ams_amg_alpha_opts[2], /* AMG relax_type */ 1324 jac->ams_amg_alpha_theta, 1325 jac->ams_amg_alpha_opts[3], /* AMG interp_type */ 1326 jac->ams_amg_alpha_opts[4])); /* AMG Pmax */ 1327 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->ams_amg_beta_opts[0], /* AMG coarsen type */ 1328 jac->ams_amg_beta_opts[1], /* AMG agg_levels */ 1329 jac->ams_amg_beta_opts[2], /* AMG relax_type */ 1330 jac->ams_amg_beta_theta, 1331 jac->ams_amg_beta_opts[3], /* AMG interp_type */ 1332 jac->ams_amg_beta_opts[4])); /* AMG Pmax */ 1333 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE_AMS);CHKERRQ(ierr); 1334 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE_AMS);CHKERRQ(ierr); 1335 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);CHKERRQ(ierr); 1336 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1337 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr); 1338 PetscFunctionReturn(0); 1339 } 1340 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 1341 1342 jac->hypre_type = NULL; 1343 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg, ams",name); 1344 PetscFunctionReturn(0); 1345 } 1346 1347 /* 1348 It only gets here if the HYPRE type has not been set before the call to 1349 ...SetFromOptions() which actually is most of the time 1350 */ 1351 #undef __FUNCT__ 1352 #define __FUNCT__ "PCSetFromOptions_HYPRE" 1353 static PetscErrorCode PCSetFromOptions_HYPRE(PC pc) 1354 { 1355 PetscErrorCode ierr; 1356 PetscInt indx; 1357 const char *type[] = {"pilut","parasails","boomeramg","euclid","ams"}; 1358 PetscBool flg; 1359 1360 PetscFunctionBegin; 1361 ierr = PetscOptionsHead("HYPRE preconditioner options");CHKERRQ(ierr); 1362 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,5,"boomeramg",&indx,&flg);CHKERRQ(ierr); 1363 if (flg) { 1364 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr); 1365 } else { 1366 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr); 1367 } 1368 if (pc->ops->setfromoptions) { 1369 ierr = pc->ops->setfromoptions(pc);CHKERRQ(ierr); 1370 } 1371 ierr = PetscOptionsTail();CHKERRQ(ierr); 1372 PetscFunctionReturn(0); 1373 } 1374 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "PCHYPRESetType" 1377 /*@C 1378 PCHYPRESetType - Sets which hypre preconditioner you wish to use 1379 1380 Input Parameters: 1381 + pc - the preconditioner context 1382 - name - either pilut, parasails, boomeramg, euclid, ams 1383 1384 Options Database Keys: 1385 -pc_hypre_type - One of pilut, parasails, boomeramg, euclid, ams 1386 1387 Level: intermediate 1388 1389 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1390 PCHYPRE 1391 1392 @*/ 1393 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 1394 { 1395 PetscErrorCode ierr; 1396 1397 PetscFunctionBegin; 1398 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1399 PetscValidCharPointer(name,2); 1400 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr); 1401 PetscFunctionReturn(0); 1402 } 1403 1404 #undef __FUNCT__ 1405 #define __FUNCT__ "PCHYPREGetType" 1406 /*@C 1407 PCHYPREGetType - Gets which hypre preconditioner you are using 1408 1409 Input Parameter: 1410 . pc - the preconditioner context 1411 1412 Output Parameter: 1413 . name - either pilut, parasails, boomeramg, euclid, ams 1414 1415 Level: intermediate 1416 1417 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC, 1418 PCHYPRE 1419 1420 @*/ 1421 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 1422 { 1423 PetscErrorCode ierr; 1424 1425 PetscFunctionBegin; 1426 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1427 PetscValidPointer(name,2); 1428 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr); 1429 PetscFunctionReturn(0); 1430 } 1431 1432 /*MC 1433 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 1434 1435 Options Database Keys: 1436 + -pc_hypre_type - One of pilut, parasails, boomeramg, euclid, ams 1437 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX 1438 preconditioner 1439 1440 Level: intermediate 1441 1442 Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()), 1443 the many hypre options can ONLY be set via the options database (e.g. the command line 1444 or with PetscOptionsSetValue(), there are no functions to set them) 1445 1446 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations 1447 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if 1448 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner 1449 (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of 1450 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations 1451 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 1452 then AT MOST twenty V-cycles of boomeramg will be called. 1453 1454 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation 1455 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 1456 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi. 1457 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly 1458 and use -ksp_max_it to control the number of V-cycles. 1459 (see the PETSc FAQ.html at the PETSc website under the Documentation tab). 1460 1461 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option 1462 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L. 1463 1464 See PCPFMG for access to the hypre Struct PFMG solver 1465 1466 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1467 PCHYPRESetType(), PCPFMG 1468 1469 M*/ 1470 1471 #undef __FUNCT__ 1472 #define __FUNCT__ "PCCreate_HYPRE" 1473 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 1474 { 1475 PC_HYPRE *jac; 1476 PetscErrorCode ierr; 1477 1478 PetscFunctionBegin; 1479 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 1480 1481 pc->data = jac; 1482 pc->ops->destroy = PCDestroy_HYPRE; 1483 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 1484 pc->ops->setup = PCSetUp_HYPRE; 1485 pc->ops->apply = PCApply_HYPRE; 1486 jac->comm_hypre = MPI_COMM_NULL; 1487 jac->hypre_type = NULL; 1488 jac->coords[0] = NULL; 1489 jac->coords[1] = NULL; 1490 jac->coords[2] = NULL; 1491 jac->constants[0] = NULL; 1492 jac->constants[1] = NULL; 1493 jac->constants[2] = NULL; 1494 /* duplicate communicator for hypre */ 1495 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr); 1496 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 1497 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 1498 PetscFunctionReturn(0); 1499 } 1500 1501 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 1502 1503 /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */ 1504 #include <petsc-private/matimpl.h> 1505 1506 typedef struct { 1507 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 1508 HYPRE_StructSolver hsolver; 1509 1510 /* keep copy of PFMG options used so may view them */ 1511 PetscInt its; 1512 double tol; 1513 PetscInt relax_type; 1514 PetscInt rap_type; 1515 PetscInt num_pre_relax,num_post_relax; 1516 PetscInt max_levels; 1517 } PC_PFMG; 1518 1519 #undef __FUNCT__ 1520 #define __FUNCT__ "PCDestroy_PFMG" 1521 PetscErrorCode PCDestroy_PFMG(PC pc) 1522 { 1523 PetscErrorCode ierr; 1524 PC_PFMG *ex = (PC_PFMG*) pc->data; 1525 1526 PetscFunctionBegin; 1527 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1528 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1529 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1530 PetscFunctionReturn(0); 1531 } 1532 1533 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 1534 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 1535 1536 #undef __FUNCT__ 1537 #define __FUNCT__ "PCView_PFMG" 1538 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 1539 { 1540 PetscErrorCode ierr; 1541 PetscBool iascii; 1542 PC_PFMG *ex = (PC_PFMG*) pc->data; 1543 1544 PetscFunctionBegin; 1545 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1546 if (iascii) { 1547 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 1548 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1549 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1550 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1551 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 1552 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1553 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);CHKERRQ(ierr); 1554 } 1555 PetscFunctionReturn(0); 1556 } 1557 1558 1559 #undef __FUNCT__ 1560 #define __FUNCT__ "PCSetFromOptions_PFMG" 1561 PetscErrorCode PCSetFromOptions_PFMG(PC pc) 1562 { 1563 PetscErrorCode ierr; 1564 PC_PFMG *ex = (PC_PFMG*) pc->data; 1565 PetscBool flg = PETSC_FALSE; 1566 1567 PetscFunctionBegin; 1568 ierr = PetscOptionsHead("PFMG options");CHKERRQ(ierr); 1569 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1570 if (flg) { 1571 int level=3; 1572 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level)); 1573 } 1574 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1575 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 1576 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); 1577 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 1578 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); 1579 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 1580 1581 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 1582 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 1583 1584 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1585 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 1586 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); 1587 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 1588 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 1589 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 1590 ierr = PetscOptionsTail();CHKERRQ(ierr); 1591 PetscFunctionReturn(0); 1592 } 1593 1594 #undef __FUNCT__ 1595 #define __FUNCT__ "PCApply_PFMG" 1596 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 1597 { 1598 PetscErrorCode ierr; 1599 PC_PFMG *ex = (PC_PFMG*) pc->data; 1600 PetscScalar *xx,*yy; 1601 PetscInt ilower[3],iupper[3]; 1602 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1603 1604 PetscFunctionBegin; 1605 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1606 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1607 iupper[0] += ilower[0] - 1; 1608 iupper[1] += ilower[1] - 1; 1609 iupper[2] += ilower[2] - 1; 1610 1611 /* copy x values over to hypre */ 1612 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 1613 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 1614 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,xx)); 1615 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 1616 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 1617 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1618 1619 /* copy solution values back to PETSc */ 1620 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1621 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 1622 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1623 PetscFunctionReturn(0); 1624 } 1625 1626 #undef __FUNCT__ 1627 #define __FUNCT__ "PCApplyRichardson_PFMG" 1628 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) 1629 { 1630 PC_PFMG *jac = (PC_PFMG*)pc->data; 1631 PetscErrorCode ierr; 1632 PetscInt oits; 1633 1634 PetscFunctionBegin; 1635 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1636 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 1637 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 1638 1639 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 1640 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 1641 *outits = oits; 1642 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1643 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1644 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 1645 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 1650 #undef __FUNCT__ 1651 #define __FUNCT__ "PCSetUp_PFMG" 1652 PetscErrorCode PCSetUp_PFMG(PC pc) 1653 { 1654 PetscErrorCode ierr; 1655 PC_PFMG *ex = (PC_PFMG*) pc->data; 1656 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 1657 PetscBool flg; 1658 1659 PetscFunctionBegin; 1660 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 1661 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 1662 1663 /* create the hypre solver object and set its information */ 1664 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1665 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1666 ierr = PCSetFromOptions_PFMG(pc);CHKERRQ(ierr); 1667 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 1668 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 1669 PetscFunctionReturn(0); 1670 } 1671 1672 1673 /*MC 1674 PCPFMG - the hypre PFMG multigrid solver 1675 1676 Level: advanced 1677 1678 Options Database: 1679 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 1680 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 1681 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 1682 . -pc_pfmg_tol <tol> tolerance of PFMG 1683 . -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 1684 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 1685 1686 Notes: This is for CELL-centered descretizations 1687 1688 This must be used with the MATHYPRESTRUCT matrix type. 1689 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 1690 1691 .seealso: PCMG, MATHYPRESTRUCT 1692 M*/ 1693 1694 #undef __FUNCT__ 1695 #define __FUNCT__ "PCCreate_PFMG" 1696 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 1697 { 1698 PetscErrorCode ierr; 1699 PC_PFMG *ex; 1700 1701 PetscFunctionBegin; 1702 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 1703 pc->data = ex; 1704 1705 ex->its = 1; 1706 ex->tol = 1.e-8; 1707 ex->relax_type = 1; 1708 ex->rap_type = 0; 1709 ex->num_pre_relax = 1; 1710 ex->num_post_relax = 1; 1711 ex->max_levels = 0; 1712 1713 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 1714 pc->ops->view = PCView_PFMG; 1715 pc->ops->destroy = PCDestroy_PFMG; 1716 pc->ops->apply = PCApply_PFMG; 1717 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 1718 pc->ops->setup = PCSetUp_PFMG; 1719 1720 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 1721 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 1722 PetscFunctionReturn(0); 1723 } 1724 1725 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 1726 1727 /* we know we are working with a HYPRE_SStructMatrix */ 1728 typedef struct { 1729 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 1730 HYPRE_SStructSolver ss_solver; 1731 1732 /* keep copy of SYSPFMG options used so may view them */ 1733 PetscInt its; 1734 double tol; 1735 PetscInt relax_type; 1736 PetscInt num_pre_relax,num_post_relax; 1737 } PC_SysPFMG; 1738 1739 #undef __FUNCT__ 1740 #define __FUNCT__ "PCDestroy_SysPFMG" 1741 PetscErrorCode PCDestroy_SysPFMG(PC pc) 1742 { 1743 PetscErrorCode ierr; 1744 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1745 1746 PetscFunctionBegin; 1747 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 1748 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1749 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1750 PetscFunctionReturn(0); 1751 } 1752 1753 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 1754 1755 #undef __FUNCT__ 1756 #define __FUNCT__ "PCView_SysPFMG" 1757 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 1758 { 1759 PetscErrorCode ierr; 1760 PetscBool iascii; 1761 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1762 1763 PetscFunctionBegin; 1764 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1765 if (iascii) { 1766 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 1767 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);CHKERRQ(ierr); 1768 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr); 1769 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1770 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1771 } 1772 PetscFunctionReturn(0); 1773 } 1774 1775 1776 #undef __FUNCT__ 1777 #define __FUNCT__ "PCSetFromOptions_SysPFMG" 1778 PetscErrorCode PCSetFromOptions_SysPFMG(PC pc) 1779 { 1780 PetscErrorCode ierr; 1781 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1782 PetscBool flg = PETSC_FALSE; 1783 1784 PetscFunctionBegin; 1785 ierr = PetscOptionsHead("SysPFMG options");CHKERRQ(ierr); 1786 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 1787 if (flg) { 1788 int level=3; 1789 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level)); 1790 } 1791 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 1792 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 1793 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); 1794 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 1795 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); 1796 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 1797 1798 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 1799 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 1800 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); 1801 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 1802 ierr = PetscOptionsTail();CHKERRQ(ierr); 1803 PetscFunctionReturn(0); 1804 } 1805 1806 #undef __FUNCT__ 1807 #define __FUNCT__ "PCApply_SysPFMG" 1808 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 1809 { 1810 PetscErrorCode ierr; 1811 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1812 PetscScalar *xx,*yy; 1813 PetscInt ilower[3],iupper[3]; 1814 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 1815 PetscInt ordering= mx->dofs_order; 1816 PetscInt nvars = mx->nvars; 1817 PetscInt part = 0; 1818 PetscInt size; 1819 PetscInt i; 1820 1821 PetscFunctionBegin; 1822 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1823 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 1824 iupper[0] += ilower[0] - 1; 1825 iupper[1] += ilower[1] - 1; 1826 iupper[2] += ilower[2] - 1; 1827 1828 size = 1; 1829 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 1830 1831 /* copy x values over to hypre for variable ordering */ 1832 if (ordering) { 1833 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1834 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 1835 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,xx+(size*i))); 1836 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 1837 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 1838 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 1839 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1840 1841 /* copy solution values back to PETSc */ 1842 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1843 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 1844 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1845 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 1846 PetscScalar *z; 1847 PetscInt j, k; 1848 1849 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 1850 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 1851 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 1852 1853 /* transform nodal to hypre's variable ordering for sys_pfmg */ 1854 for (i= 0; i< size; i++) { 1855 k= i*nvars; 1856 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 1857 } 1858 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 1859 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 1860 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 1861 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1862 1863 /* copy solution values back to PETSc */ 1864 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 1865 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 1866 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 1867 for (i= 0; i< size; i++) { 1868 k= i*nvars; 1869 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 1870 } 1871 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 1872 ierr = PetscFree(z);CHKERRQ(ierr); 1873 } 1874 PetscFunctionReturn(0); 1875 } 1876 1877 #undef __FUNCT__ 1878 #define __FUNCT__ "PCApplyRichardson_SysPFMG" 1879 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) 1880 { 1881 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 1882 PetscErrorCode ierr; 1883 PetscInt oits; 1884 1885 PetscFunctionBegin; 1886 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 1887 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 1888 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 1889 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 1890 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits)); 1891 *outits = oits; 1892 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 1893 else *reason = PCRICHARDSON_CONVERGED_RTOL; 1894 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 1895 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 1896 PetscFunctionReturn(0); 1897 } 1898 1899 1900 #undef __FUNCT__ 1901 #define __FUNCT__ "PCSetUp_SysPFMG" 1902 PetscErrorCode PCSetUp_SysPFMG(PC pc) 1903 { 1904 PetscErrorCode ierr; 1905 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 1906 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 1907 PetscBool flg; 1908 1909 PetscFunctionBegin; 1910 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 1911 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 1912 1913 /* create the hypre sstruct solver object and set its information */ 1914 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 1915 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 1916 ierr = PCSetFromOptions_SysPFMG(pc);CHKERRQ(ierr); 1917 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 1918 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 1919 PetscFunctionReturn(0); 1920 } 1921 1922 1923 /*MC 1924 PCSysPFMG - the hypre SysPFMG multigrid solver 1925 1926 Level: advanced 1927 1928 Options Database: 1929 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 1930 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 1931 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 1932 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 1933 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 1934 1935 Notes: This is for CELL-centered descretizations 1936 1937 This must be used with the MATHYPRESSTRUCT matrix type. 1938 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 1939 Also, only cell-centered variables. 1940 1941 .seealso: PCMG, MATHYPRESSTRUCT 1942 M*/ 1943 1944 #undef __FUNCT__ 1945 #define __FUNCT__ "PCCreate_SysPFMG" 1946 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 1947 { 1948 PetscErrorCode ierr; 1949 PC_SysPFMG *ex; 1950 1951 PetscFunctionBegin; 1952 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 1953 pc->data = ex; 1954 1955 ex->its = 1; 1956 ex->tol = 1.e-8; 1957 ex->relax_type = 1; 1958 ex->num_pre_relax = 1; 1959 ex->num_post_relax = 1; 1960 1961 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 1962 pc->ops->view = PCView_SysPFMG; 1963 pc->ops->destroy = PCDestroy_SysPFMG; 1964 pc->ops->apply = PCApply_SysPFMG; 1965 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 1966 pc->ops->setup = PCSetUp_SysPFMG; 1967 1968 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 1969 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 1970 PetscFunctionReturn(0); 1971 } 1972