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