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