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