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