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