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