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