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