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