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