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