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