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 */, "" /* 13 */, "" /* 14 */, 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: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh 1175 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 1176 1177 .seealso: 1178 @*/ 1179 PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G) 1180 { 1181 PetscErrorCode ierr; 1182 1183 PetscFunctionBegin; 1184 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1185 PetscValidHeaderSpecific(G,MAT_CLASSID,2); 1186 PetscCheckSameComm(pc,1,G,2); 1187 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));CHKERRQ(ierr); 1188 PetscFunctionReturn(0); 1189 } 1190 1191 static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C) 1192 { 1193 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1194 PetscBool ishypre; 1195 PetscErrorCode ierr; 1196 1197 PetscFunctionBegin; 1198 ierr = PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);CHKERRQ(ierr); 1199 if (ishypre) { 1200 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr); 1201 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1202 jac->C = C; 1203 } else { 1204 ierr = MatDestroy(&jac->C);CHKERRQ(ierr); 1205 ierr = MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr); 1206 } 1207 PetscFunctionReturn(0); 1208 } 1209 1210 /*@ 1211 PCHYPRESetDiscreteCurl - Set discrete curl matrix 1212 1213 Collective on PC 1214 1215 Input Parameters: 1216 + pc - the preconditioning context 1217 - C - the discrete curl 1218 1219 Level: intermediate 1220 1221 Notes: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh 1222 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 1223 1224 .seealso: 1225 @*/ 1226 PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C) 1227 { 1228 PetscErrorCode ierr; 1229 1230 PetscFunctionBegin; 1231 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1232 PetscValidHeaderSpecific(C,MAT_CLASSID,2); 1233 PetscCheckSameComm(pc,1,C,2); 1234 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));CHKERRQ(ierr); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) 1239 { 1240 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1241 PetscBool ishypre; 1242 PetscErrorCode ierr; 1243 PetscInt i; 1244 PetscFunctionBegin; 1245 1246 ierr = MatDestroy(&jac->RT_PiFull);CHKERRQ(ierr); 1247 ierr = MatDestroy(&jac->ND_PiFull);CHKERRQ(ierr); 1248 for (i=0;i<3;++i) { 1249 ierr = MatDestroy(&jac->RT_Pi[i]);CHKERRQ(ierr); 1250 ierr = MatDestroy(&jac->ND_Pi[i]);CHKERRQ(ierr); 1251 } 1252 1253 jac->dim = dim; 1254 if (RT_PiFull) { 1255 ierr = PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr); 1256 if (ishypre) { 1257 ierr = PetscObjectReference((PetscObject)RT_PiFull);CHKERRQ(ierr); 1258 jac->RT_PiFull = RT_PiFull; 1259 } else { 1260 ierr = MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);CHKERRQ(ierr); 1261 } 1262 } 1263 if (RT_Pi) { 1264 for (i=0;i<dim;++i) { 1265 if (RT_Pi[i]) { 1266 ierr = PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr); 1267 if (ishypre) { 1268 ierr = PetscObjectReference((PetscObject)RT_Pi[i]);CHKERRQ(ierr); 1269 jac->RT_Pi[i] = RT_Pi[i]; 1270 } else { 1271 ierr = MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);CHKERRQ(ierr); 1272 } 1273 } 1274 } 1275 } 1276 if (ND_PiFull) { 1277 ierr = PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr); 1278 if (ishypre) { 1279 ierr = PetscObjectReference((PetscObject)ND_PiFull);CHKERRQ(ierr); 1280 jac->ND_PiFull = ND_PiFull; 1281 } else { 1282 ierr = MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);CHKERRQ(ierr); 1283 } 1284 } 1285 if (ND_Pi) { 1286 for (i=0;i<dim;++i) { 1287 if (ND_Pi[i]) { 1288 ierr = PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr); 1289 if (ishypre) { 1290 ierr = PetscObjectReference((PetscObject)ND_Pi[i]);CHKERRQ(ierr); 1291 jac->ND_Pi[i] = ND_Pi[i]; 1292 } else { 1293 ierr = MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);CHKERRQ(ierr); 1294 } 1295 } 1296 } 1297 } 1298 1299 PetscFunctionReturn(0); 1300 } 1301 1302 /*@ 1303 PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner 1304 1305 Collective on PC 1306 1307 Input Parameters: 1308 + pc - the preconditioning context 1309 - dim - the dimension of the problem, only used in AMS 1310 - RT_PiFull - Raviart-Thomas interpolation matrix 1311 - RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix 1312 - ND_PiFull - Nedelec interpolation matrix 1313 - ND_Pi - x/y/z component of Nedelec interpolation matrix 1314 1315 Notes: For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL. 1316 For ADS, both type of interpolation matrices are needed. 1317 Level: intermediate 1318 1319 .seealso: 1320 @*/ 1321 PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) 1322 { 1323 PetscErrorCode ierr; 1324 PetscInt i; 1325 1326 PetscFunctionBegin; 1327 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1328 if (RT_PiFull) { 1329 PetscValidHeaderSpecific(RT_PiFull,MAT_CLASSID,3); 1330 PetscCheckSameComm(pc,1,RT_PiFull,3); 1331 } 1332 if (RT_Pi) { 1333 PetscValidPointer(RT_Pi,4); 1334 for (i=0;i<dim;++i) { 1335 if (RT_Pi[i]) { 1336 PetscValidHeaderSpecific(RT_Pi[i],MAT_CLASSID,4); 1337 PetscCheckSameComm(pc,1,RT_Pi[i],4); 1338 } 1339 } 1340 } 1341 if (ND_PiFull) { 1342 PetscValidHeaderSpecific(ND_PiFull,MAT_CLASSID,5); 1343 PetscCheckSameComm(pc,1,ND_PiFull,5); 1344 } 1345 if (ND_Pi) { 1346 PetscValidPointer(ND_Pi,6); 1347 for (i=0;i<dim;++i) { 1348 if (ND_Pi[i]) { 1349 PetscValidHeaderSpecific(ND_Pi[i],MAT_CLASSID,6); 1350 PetscCheckSameComm(pc,1,ND_Pi[i],6); 1351 } 1352 } 1353 } 1354 ierr = PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));CHKERRQ(ierr); 1355 PetscFunctionReturn(0); 1356 } 1357 1358 static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha) 1359 { 1360 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1361 PetscBool ishypre; 1362 PetscErrorCode ierr; 1363 1364 PetscFunctionBegin; 1365 ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr); 1366 if (ishypre) { 1367 if (isalpha) { 1368 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 1369 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr); 1370 jac->alpha_Poisson = A; 1371 } else { 1372 if (A) { 1373 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 1374 } else { 1375 jac->ams_beta_is_zero = PETSC_TRUE; 1376 } 1377 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1378 jac->beta_Poisson = A; 1379 } 1380 } else { 1381 if (isalpha) { 1382 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr); 1383 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);CHKERRQ(ierr); 1384 } else { 1385 if (A) { 1386 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1387 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);CHKERRQ(ierr); 1388 } else { 1389 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr); 1390 jac->ams_beta_is_zero = PETSC_TRUE; 1391 } 1392 } 1393 } 1394 PetscFunctionReturn(0); 1395 } 1396 1397 /*@ 1398 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix 1399 1400 Collective on PC 1401 1402 Input Parameters: 1403 + pc - the preconditioning context 1404 - A - the matrix 1405 1406 Level: intermediate 1407 1408 Notes: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements 1409 1410 .seealso: 1411 @*/ 1412 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) 1413 { 1414 PetscErrorCode ierr; 1415 1416 PetscFunctionBegin; 1417 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1418 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 1419 PetscCheckSameComm(pc,1,A,2); 1420 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));CHKERRQ(ierr); 1421 PetscFunctionReturn(0); 1422 } 1423 1424 /*@ 1425 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix 1426 1427 Collective on PC 1428 1429 Input Parameters: 1430 + pc - the preconditioning context 1431 - A - the matrix 1432 1433 Level: intermediate 1434 1435 Notes: A should be obtained by discretizing the Poisson problem with linear finite elements. 1436 Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL. 1437 1438 .seealso: 1439 @*/ 1440 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) 1441 { 1442 PetscErrorCode ierr; 1443 1444 PetscFunctionBegin; 1445 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1446 if (A) { 1447 PetscValidHeaderSpecific(A,MAT_CLASSID,2); 1448 PetscCheckSameComm(pc,1,A,2); 1449 } 1450 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));CHKERRQ(ierr); 1451 PetscFunctionReturn(0); 1452 } 1453 1454 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo) 1455 { 1456 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1457 PetscErrorCode ierr; 1458 1459 PetscFunctionBegin; 1460 /* throw away any vector if already set */ 1461 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); 1462 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); 1463 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); 1464 jac->constants[0] = NULL; 1465 jac->constants[1] = NULL; 1466 jac->constants[2] = NULL; 1467 ierr = VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);CHKERRQ(ierr); 1468 ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr); 1469 ierr = VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);CHKERRQ(ierr); 1470 ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr); 1471 jac->dim = 2; 1472 if (zzo) { 1473 ierr = VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);CHKERRQ(ierr); 1474 ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr); 1475 jac->dim++; 1476 } 1477 PetscFunctionReturn(0); 1478 } 1479 1480 /*@ 1481 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis 1482 1483 Collective on PC 1484 1485 Input Parameters: 1486 + pc - the preconditioning context 1487 - ozz - vector representing (1,0,0) (or (1,0) in 2D) 1488 - zoz - vector representing (0,1,0) (or (0,1) in 2D) 1489 - zzo - vector representing (0,0,1) (use NULL in 2D) 1490 1491 Level: intermediate 1492 1493 Notes: 1494 1495 .seealso: 1496 @*/ 1497 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) 1498 { 1499 PetscErrorCode ierr; 1500 1501 PetscFunctionBegin; 1502 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1503 PetscValidHeaderSpecific(ozz,VEC_CLASSID,2); 1504 PetscValidHeaderSpecific(zoz,VEC_CLASSID,3); 1505 if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4); 1506 PetscCheckSameComm(pc,1,ozz,2); 1507 PetscCheckSameComm(pc,1,zoz,3); 1508 if (zzo) PetscCheckSameComm(pc,1,zzo,4); 1509 ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr); 1510 PetscFunctionReturn(0); 1511 } 1512 1513 static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) 1514 { 1515 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1516 Vec tv; 1517 PetscInt i; 1518 PetscErrorCode ierr; 1519 1520 PetscFunctionBegin; 1521 /* throw away any coordinate vector if already set */ 1522 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); 1523 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); 1524 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); 1525 jac->dim = dim; 1526 1527 /* compute IJ vector for coordinates */ 1528 ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr); 1529 ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr); 1530 ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr); 1531 for (i=0;i<dim;i++) { 1532 PetscScalar *array; 1533 PetscInt j; 1534 1535 ierr = VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);CHKERRQ(ierr); 1536 ierr = VecGetArray(tv,&array);CHKERRQ(ierr); 1537 for (j=0;j<nloc;j++) { 1538 array[j] = coords[j*dim+i]; 1539 } 1540 PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array)); 1541 PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i])); 1542 ierr = VecRestoreArray(tv,&array);CHKERRQ(ierr); 1543 } 1544 ierr = VecDestroy(&tv);CHKERRQ(ierr); 1545 PetscFunctionReturn(0); 1546 } 1547 1548 /* ---------------------------------------------------------------------------------*/ 1549 1550 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[]) 1551 { 1552 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1553 1554 PetscFunctionBegin; 1555 *name = jac->hypre_type; 1556 PetscFunctionReturn(0); 1557 } 1558 1559 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[]) 1560 { 1561 PC_HYPRE *jac = (PC_HYPRE*)pc->data; 1562 PetscErrorCode ierr; 1563 PetscBool flag; 1564 1565 PetscFunctionBegin; 1566 if (jac->hypre_type) { 1567 ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr); 1568 if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set"); 1569 PetscFunctionReturn(0); 1570 } else { 1571 ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr); 1572 } 1573 1574 jac->maxiter = PETSC_DEFAULT; 1575 jac->tol = PETSC_DEFAULT; 1576 jac->printstatistics = PetscLogPrintInfo; 1577 1578 ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr); 1579 if (flag) { 1580 PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver)); 1581 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut; 1582 pc->ops->view = PCView_HYPRE_Pilut; 1583 jac->destroy = HYPRE_ParCSRPilutDestroy; 1584 jac->setup = HYPRE_ParCSRPilutSetup; 1585 jac->solve = HYPRE_ParCSRPilutSolve; 1586 jac->factorrowsize = PETSC_DEFAULT; 1587 PetscFunctionReturn(0); 1588 } 1589 ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr); 1590 if (flag) { 1591 PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver)); 1592 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails; 1593 pc->ops->view = PCView_HYPRE_ParaSails; 1594 jac->destroy = HYPRE_ParaSailsDestroy; 1595 jac->setup = HYPRE_ParaSailsSetup; 1596 jac->solve = HYPRE_ParaSailsSolve; 1597 /* initialize */ 1598 jac->nlevels = 1; 1599 jac->threshhold = .1; 1600 jac->filter = .1; 1601 jac->loadbal = 0; 1602 if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE; 1603 else jac->logging = (int) PETSC_FALSE; 1604 1605 jac->ruse = (int) PETSC_FALSE; 1606 jac->symt = 0; 1607 PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels)); 1608 PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter)); 1609 PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal)); 1610 PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging)); 1611 PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse)); 1612 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt)); 1613 PetscFunctionReturn(0); 1614 } 1615 ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr); 1616 if (flag) { 1617 ierr = HYPRE_BoomerAMGCreate(&jac->hsolver); 1618 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG; 1619 pc->ops->view = PCView_HYPRE_BoomerAMG; 1620 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG; 1621 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG; 1622 jac->destroy = HYPRE_BoomerAMGDestroy; 1623 jac->setup = HYPRE_BoomerAMGSetup; 1624 jac->solve = HYPRE_BoomerAMGSolve; 1625 jac->applyrichardson = PETSC_FALSE; 1626 /* these defaults match the hypre defaults */ 1627 jac->cycletype = 1; 1628 jac->maxlevels = 25; 1629 jac->maxiter = 1; 1630 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */ 1631 jac->truncfactor = 0.0; 1632 jac->strongthreshold = .25; 1633 jac->maxrowsum = .9; 1634 jac->coarsentype = 6; 1635 jac->measuretype = 0; 1636 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1; 1637 jac->smoothtype = -1; /* Not set by default */ 1638 jac->smoothnumlevels = 25; 1639 jac->eu_level = 0; 1640 jac->eu_droptolerance = 0; 1641 jac->eu_bj = 0; 1642 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */ 1643 jac->relaxtype[2] = 9; /*G.E. */ 1644 jac->relaxweight = 1.0; 1645 jac->outerrelaxweight = 1.0; 1646 jac->relaxorder = 1; 1647 jac->interptype = 0; 1648 jac->agg_nl = 0; 1649 jac->pmax = 0; 1650 jac->truncfactor = 0.0; 1651 jac->agg_num_paths = 1; 1652 1653 jac->nodal_coarsen = 0; 1654 jac->nodal_relax = PETSC_FALSE; 1655 jac->nodal_relax_levels = 1; 1656 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype)); 1657 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels)); 1658 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter)); 1659 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol)); 1660 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor)); 1661 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold)); 1662 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum)); 1663 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype)); 1664 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype)); 1665 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder)); 1666 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype)); 1667 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl)); 1668 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax)); 1669 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths)); 1670 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/ 1671 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */ 1672 PetscFunctionReturn(0); 1673 } 1674 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr); 1675 if (flag) { 1676 ierr = HYPRE_AMSCreate(&jac->hsolver); 1677 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS; 1678 pc->ops->view = PCView_HYPRE_AMS; 1679 jac->destroy = HYPRE_AMSDestroy; 1680 jac->setup = HYPRE_AMSSetup; 1681 jac->solve = HYPRE_AMSSolve; 1682 jac->coords[0] = NULL; 1683 jac->coords[1] = NULL; 1684 jac->coords[2] = NULL; 1685 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */ 1686 jac->as_print = 0; 1687 jac->as_max_iter = 1; /* used as a preconditioner */ 1688 jac->as_tol = 0.; /* used as a preconditioner */ 1689 jac->ams_cycle_type = 13; 1690 /* Smoothing options */ 1691 jac->as_relax_type = 2; 1692 jac->as_relax_times = 1; 1693 jac->as_relax_weight = 1.0; 1694 jac->as_omega = 1.0; 1695 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1696 jac->as_amg_alpha_opts[0] = 10; 1697 jac->as_amg_alpha_opts[1] = 1; 1698 jac->as_amg_alpha_opts[2] = 6; 1699 jac->as_amg_alpha_opts[3] = 6; 1700 jac->as_amg_alpha_opts[4] = 4; 1701 jac->as_amg_alpha_theta = 0.25; 1702 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1703 jac->as_amg_beta_opts[0] = 10; 1704 jac->as_amg_beta_opts[1] = 1; 1705 jac->as_amg_beta_opts[2] = 6; 1706 jac->as_amg_beta_opts[3] = 6; 1707 jac->as_amg_beta_opts[4] = 4; 1708 jac->as_amg_beta_theta = 0.25; 1709 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print)); 1710 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1711 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1712 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol)); 1713 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1714 jac->as_relax_times, 1715 jac->as_relax_weight, 1716 jac->as_omega)); 1717 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1718 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1719 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1720 jac->as_amg_alpha_theta, 1721 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1722 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1723 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1724 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1725 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1726 jac->as_amg_beta_theta, 1727 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1728 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1729 /* Zero conductivity */ 1730 jac->ams_beta_is_zero = PETSC_FALSE; 1731 jac->ams_beta_is_zero_part = PETSC_FALSE; 1732 PetscFunctionReturn(0); 1733 } 1734 ierr = PetscStrcmp("ads",jac->hypre_type,&flag);CHKERRQ(ierr); 1735 if (flag) { 1736 ierr = HYPRE_ADSCreate(&jac->hsolver); 1737 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS; 1738 pc->ops->view = PCView_HYPRE_ADS; 1739 jac->destroy = HYPRE_ADSDestroy; 1740 jac->setup = HYPRE_ADSSetup; 1741 jac->solve = HYPRE_ADSSolve; 1742 jac->coords[0] = NULL; 1743 jac->coords[1] = NULL; 1744 jac->coords[2] = NULL; 1745 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */ 1746 jac->as_print = 0; 1747 jac->as_max_iter = 1; /* used as a preconditioner */ 1748 jac->as_tol = 0.; /* used as a preconditioner */ 1749 jac->ads_cycle_type = 13; 1750 /* Smoothing options */ 1751 jac->as_relax_type = 2; 1752 jac->as_relax_times = 1; 1753 jac->as_relax_weight = 1.0; 1754 jac->as_omega = 1.0; 1755 /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1756 jac->ams_cycle_type = 14; 1757 jac->as_amg_alpha_opts[0] = 10; 1758 jac->as_amg_alpha_opts[1] = 1; 1759 jac->as_amg_alpha_opts[2] = 6; 1760 jac->as_amg_alpha_opts[3] = 6; 1761 jac->as_amg_alpha_opts[4] = 4; 1762 jac->as_amg_alpha_theta = 0.25; 1763 /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */ 1764 jac->as_amg_beta_opts[0] = 10; 1765 jac->as_amg_beta_opts[1] = 1; 1766 jac->as_amg_beta_opts[2] = 6; 1767 jac->as_amg_beta_opts[3] = 6; 1768 jac->as_amg_beta_opts[4] = 4; 1769 jac->as_amg_beta_theta = 0.25; 1770 PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print)); 1771 PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter)); 1772 PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type)); 1773 PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol)); 1774 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type, 1775 jac->as_relax_times, 1776 jac->as_relax_weight, 1777 jac->as_omega)); 1778 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */ 1779 jac->as_amg_alpha_opts[0], /* AMG coarsen type */ 1780 jac->as_amg_alpha_opts[1], /* AMG agg_levels */ 1781 jac->as_amg_alpha_opts[2], /* AMG relax_type */ 1782 jac->as_amg_alpha_theta, 1783 jac->as_amg_alpha_opts[3], /* AMG interp_type */ 1784 jac->as_amg_alpha_opts[4])); /* AMG Pmax */ 1785 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */ 1786 jac->as_amg_beta_opts[1], /* AMG agg_levels */ 1787 jac->as_amg_beta_opts[2], /* AMG relax_type */ 1788 jac->as_amg_beta_theta, 1789 jac->as_amg_beta_opts[3], /* AMG interp_type */ 1790 jac->as_amg_beta_opts[4])); /* AMG Pmax */ 1791 PetscFunctionReturn(0); 1792 } 1793 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr); 1794 1795 jac->hypre_type = NULL; 1796 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name); 1797 PetscFunctionReturn(0); 1798 } 1799 1800 /* 1801 It only gets here if the HYPRE type has not been set before the call to 1802 ...SetFromOptions() which actually is most of the time 1803 */ 1804 static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc) 1805 { 1806 PetscErrorCode ierr; 1807 PetscInt indx; 1808 const char *type[] = {"pilut","parasails","boomeramg","ams","ads"}; 1809 PetscBool flg; 1810 1811 PetscFunctionBegin; 1812 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr); 1813 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);CHKERRQ(ierr); 1814 if (flg) { 1815 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr); 1816 } else { 1817 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr); 1818 } 1819 if (pc->ops->setfromoptions) { 1820 ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr); 1821 } 1822 ierr = PetscOptionsTail();CHKERRQ(ierr); 1823 PetscFunctionReturn(0); 1824 } 1825 1826 /*@C 1827 PCHYPRESetType - Sets which hypre preconditioner you wish to use 1828 1829 Input Parameters: 1830 + pc - the preconditioner context 1831 - name - either pilut, parasails, boomeramg, ams, ads 1832 1833 Options Database Keys: 1834 -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads 1835 1836 Level: intermediate 1837 1838 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1839 PCHYPRE 1840 1841 @*/ 1842 PetscErrorCode PCHYPRESetType(PC pc,const char name[]) 1843 { 1844 PetscErrorCode ierr; 1845 1846 PetscFunctionBegin; 1847 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1848 PetscValidCharPointer(name,2); 1849 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr); 1850 PetscFunctionReturn(0); 1851 } 1852 1853 /*@C 1854 PCHYPREGetType - Gets which hypre preconditioner you are using 1855 1856 Input Parameter: 1857 . pc - the preconditioner context 1858 1859 Output Parameter: 1860 . name - either pilut, parasails, boomeramg, ams, ads 1861 1862 Level: intermediate 1863 1864 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC, 1865 PCHYPRE 1866 1867 @*/ 1868 PetscErrorCode PCHYPREGetType(PC pc,const char *name[]) 1869 { 1870 PetscErrorCode ierr; 1871 1872 PetscFunctionBegin; 1873 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 1874 PetscValidPointer(name,2); 1875 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr); 1876 PetscFunctionReturn(0); 1877 } 1878 1879 /*MC 1880 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre 1881 1882 Options Database Keys: 1883 + -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads 1884 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX 1885 preconditioner 1886 1887 Level: intermediate 1888 1889 Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()), 1890 the many hypre options can ONLY be set via the options database (e.g. the command line 1891 or with PetscOptionsSetValue(), there are no functions to set them) 1892 1893 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations 1894 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if 1895 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner 1896 (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of 1897 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations 1898 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 1899 then AT MOST twenty V-cycles of boomeramg will be called. 1900 1901 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation 1902 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry. 1903 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi. 1904 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly 1905 and use -ksp_max_it to control the number of V-cycles. 1906 (see the PETSc FAQ.html at the PETSc website under the Documentation tab). 1907 1908 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option 1909 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L. 1910 1911 MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use 1912 the two options: 1913 + -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal()) 1914 - -pc_hypre_boomeramg_vec_interp_variant <v> where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant()) 1915 1916 Depending on the linear system you may see the same or different convergence depending on the values you use. 1917 1918 See PCPFMG for access to the hypre Struct PFMG solver 1919 1920 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, 1921 PCHYPRESetType(), PCPFMG 1922 1923 M*/ 1924 1925 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) 1926 { 1927 PC_HYPRE *jac; 1928 PetscErrorCode ierr; 1929 1930 PetscFunctionBegin; 1931 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 1932 1933 pc->data = jac; 1934 pc->ops->reset = PCReset_HYPRE; 1935 pc->ops->destroy = PCDestroy_HYPRE; 1936 pc->ops->setfromoptions = PCSetFromOptions_HYPRE; 1937 pc->ops->setup = PCSetUp_HYPRE; 1938 pc->ops->apply = PCApply_HYPRE; 1939 jac->comm_hypre = MPI_COMM_NULL; 1940 jac->G = NULL; 1941 jac->C = NULL; 1942 jac->alpha_Poisson = NULL; 1943 jac->beta_Poisson = NULL; 1944 jac->coords[0] = NULL; 1945 jac->coords[1] = NULL; 1946 jac->coords[2] = NULL; 1947 jac->constants[0] = NULL; 1948 jac->constants[1] = NULL; 1949 jac->constants[2] = NULL; 1950 jac->RT_PiFull = NULL; 1951 jac->RT_Pi[0] = NULL; 1952 jac->RT_Pi[1] = NULL; 1953 jac->RT_Pi[2] = NULL; 1954 jac->ND_PiFull = NULL; 1955 jac->ND_Pi[0] = NULL; 1956 jac->ND_Pi[1] = NULL; 1957 jac->ND_Pi[2] = NULL; 1958 jac->hmnull = NULL; 1959 jac->n_hmnull = 0; 1960 /* duplicate communicator for hypre */ 1961 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr); 1962 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 1963 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 1964 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr); 1965 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr); 1966 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr); 1967 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);CHKERRQ(ierr); 1968 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);CHKERRQ(ierr); 1970 PetscFunctionReturn(0); 1971 } 1972 1973 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 1974 1975 typedef struct { 1976 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 1977 HYPRE_StructSolver hsolver; 1978 1979 /* keep copy of PFMG options used so may view them */ 1980 PetscInt its; 1981 double tol; 1982 PetscInt relax_type; 1983 PetscInt rap_type; 1984 PetscInt num_pre_relax,num_post_relax; 1985 PetscInt max_levels; 1986 } PC_PFMG; 1987 1988 PetscErrorCode PCDestroy_PFMG(PC pc) 1989 { 1990 PetscErrorCode ierr; 1991 PC_PFMG *ex = (PC_PFMG*) pc->data; 1992 1993 PetscFunctionBegin; 1994 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1995 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1996 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1997 PetscFunctionReturn(0); 1998 } 1999 2000 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 2001 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 2002 2003 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 2004 { 2005 PetscErrorCode ierr; 2006 PetscBool iascii; 2007 PC_PFMG *ex = (PC_PFMG*) pc->data; 2008 2009 PetscFunctionBegin; 2010 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2011 if (iascii) { 2012 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 2013 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 2014 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 2015 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 2016 ierr = PetscViewerASCIIPrintf(viewer," RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 2017 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 2018 ierr = PetscViewerASCIIPrintf(viewer," max levels %d\n",ex->max_levels);CHKERRQ(ierr); 2019 } 2020 PetscFunctionReturn(0); 2021 } 2022 2023 2024 PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2025 { 2026 PetscErrorCode ierr; 2027 PC_PFMG *ex = (PC_PFMG*) pc->data; 2028 PetscBool flg = PETSC_FALSE; 2029 2030 PetscFunctionBegin; 2031 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr); 2032 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2033 if (flg) { 2034 int level=3; 2035 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level)); 2036 } 2037 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2038 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 2039 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); 2040 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 2041 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); 2042 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 2043 2044 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 2045 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 2046 2047 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2048 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 2049 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); 2050 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 2051 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 2052 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 2053 ierr = PetscOptionsTail();CHKERRQ(ierr); 2054 PetscFunctionReturn(0); 2055 } 2056 2057 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 2058 { 2059 PetscErrorCode ierr; 2060 PC_PFMG *ex = (PC_PFMG*) pc->data; 2061 PetscScalar *yy; 2062 const PetscScalar *xx; 2063 PetscInt ilower[3],iupper[3]; 2064 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2065 2066 PetscFunctionBegin; 2067 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2068 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2069 iupper[0] += ilower[0] - 1; 2070 iupper[1] += ilower[1] - 1; 2071 iupper[2] += ilower[2] - 1; 2072 2073 /* copy x values over to hypre */ 2074 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 2075 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2076 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx)); 2077 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2078 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 2079 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2080 2081 /* copy solution values back to PETSc */ 2082 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2083 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 2084 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2085 PetscFunctionReturn(0); 2086 } 2087 2088 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) 2089 { 2090 PC_PFMG *jac = (PC_PFMG*)pc->data; 2091 PetscErrorCode ierr; 2092 PetscInt oits; 2093 2094 PetscFunctionBegin; 2095 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2096 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 2097 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 2098 2099 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 2100 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 2101 *outits = oits; 2102 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2103 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2104 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 2105 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 2106 PetscFunctionReturn(0); 2107 } 2108 2109 2110 PetscErrorCode PCSetUp_PFMG(PC pc) 2111 { 2112 PetscErrorCode ierr; 2113 PC_PFMG *ex = (PC_PFMG*) pc->data; 2114 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2115 PetscBool flg; 2116 2117 PetscFunctionBegin; 2118 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 2119 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 2120 2121 /* create the hypre solver object and set its information */ 2122 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 2123 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2124 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2125 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 2126 PetscFunctionReturn(0); 2127 } 2128 2129 2130 /*MC 2131 PCPFMG - the hypre PFMG multigrid solver 2132 2133 Level: advanced 2134 2135 Options Database: 2136 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 2137 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2138 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2139 . -pc_pfmg_tol <tol> tolerance of PFMG 2140 . -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 2141 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 2142 2143 Notes: This is for CELL-centered descretizations 2144 2145 This must be used with the MATHYPRESTRUCT matrix type. 2146 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 2147 2148 .seealso: PCMG, MATHYPRESTRUCT 2149 M*/ 2150 2151 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 2152 { 2153 PetscErrorCode ierr; 2154 PC_PFMG *ex; 2155 2156 PetscFunctionBegin; 2157 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2158 pc->data = ex; 2159 2160 ex->its = 1; 2161 ex->tol = 1.e-8; 2162 ex->relax_type = 1; 2163 ex->rap_type = 0; 2164 ex->num_pre_relax = 1; 2165 ex->num_post_relax = 1; 2166 ex->max_levels = 0; 2167 2168 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 2169 pc->ops->view = PCView_PFMG; 2170 pc->ops->destroy = PCDestroy_PFMG; 2171 pc->ops->apply = PCApply_PFMG; 2172 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 2173 pc->ops->setup = PCSetUp_PFMG; 2174 2175 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 2176 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2177 PetscFunctionReturn(0); 2178 } 2179 2180 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 2181 2182 /* we know we are working with a HYPRE_SStructMatrix */ 2183 typedef struct { 2184 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 2185 HYPRE_SStructSolver ss_solver; 2186 2187 /* keep copy of SYSPFMG options used so may view them */ 2188 PetscInt its; 2189 double tol; 2190 PetscInt relax_type; 2191 PetscInt num_pre_relax,num_post_relax; 2192 } PC_SysPFMG; 2193 2194 PetscErrorCode PCDestroy_SysPFMG(PC pc) 2195 { 2196 PetscErrorCode ierr; 2197 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2198 2199 PetscFunctionBegin; 2200 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2201 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 2202 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2203 PetscFunctionReturn(0); 2204 } 2205 2206 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 2207 2208 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 2209 { 2210 PetscErrorCode ierr; 2211 PetscBool iascii; 2212 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2213 2214 PetscFunctionBegin; 2215 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2216 if (iascii) { 2217 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 2218 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 2219 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 2220 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 2221 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 2222 } 2223 PetscFunctionReturn(0); 2224 } 2225 2226 2227 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2228 { 2229 PetscErrorCode ierr; 2230 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2231 PetscBool flg = PETSC_FALSE; 2232 2233 PetscFunctionBegin; 2234 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr); 2235 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2236 if (flg) { 2237 int level=3; 2238 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level)); 2239 } 2240 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2241 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 2242 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); 2243 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 2244 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); 2245 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 2246 2247 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2248 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 2249 ierr = PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr); 2250 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 2251 ierr = PetscOptionsTail();CHKERRQ(ierr); 2252 PetscFunctionReturn(0); 2253 } 2254 2255 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 2256 { 2257 PetscErrorCode ierr; 2258 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2259 PetscScalar *yy; 2260 const PetscScalar *xx; 2261 PetscInt ilower[3],iupper[3]; 2262 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2263 PetscInt ordering= mx->dofs_order; 2264 PetscInt nvars = mx->nvars; 2265 PetscInt part = 0; 2266 PetscInt size; 2267 PetscInt i; 2268 2269 PetscFunctionBegin; 2270 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2271 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2272 iupper[0] += ilower[0] - 1; 2273 iupper[1] += ilower[1] - 1; 2274 iupper[2] += ilower[2] - 1; 2275 2276 size = 1; 2277 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 2278 2279 /* copy x values over to hypre for variable ordering */ 2280 if (ordering) { 2281 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2282 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2283 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i))); 2284 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2285 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2286 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 2287 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2288 2289 /* copy solution values back to PETSc */ 2290 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2291 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 2292 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2293 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2294 PetscScalar *z; 2295 PetscInt j, k; 2296 2297 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 2298 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2299 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2300 2301 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2302 for (i= 0; i< size; i++) { 2303 k= i*nvars; 2304 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 2305 } 2306 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2307 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2308 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2309 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2310 2311 /* copy solution values back to PETSc */ 2312 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2313 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2314 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2315 for (i= 0; i< size; i++) { 2316 k= i*nvars; 2317 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 2318 } 2319 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2320 ierr = PetscFree(z);CHKERRQ(ierr); 2321 } 2322 PetscFunctionReturn(0); 2323 } 2324 2325 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) 2326 { 2327 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 2328 PetscErrorCode ierr; 2329 PetscInt oits; 2330 2331 PetscFunctionBegin; 2332 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2333 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 2334 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 2335 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 2336 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits)); 2337 *outits = oits; 2338 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2339 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2340 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 2341 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 2342 PetscFunctionReturn(0); 2343 } 2344 2345 2346 PetscErrorCode PCSetUp_SysPFMG(PC pc) 2347 { 2348 PetscErrorCode ierr; 2349 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2350 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2351 PetscBool flg; 2352 2353 PetscFunctionBegin; 2354 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 2355 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 2356 2357 /* create the hypre sstruct solver object and set its information */ 2358 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2359 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2360 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 2361 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2362 PetscFunctionReturn(0); 2363 } 2364 2365 2366 /*MC 2367 PCSysPFMG - the hypre SysPFMG multigrid solver 2368 2369 Level: advanced 2370 2371 Options Database: 2372 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 2373 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2374 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2375 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 2376 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 2377 2378 Notes: This is for CELL-centered descretizations 2379 2380 This must be used with the MATHYPRESSTRUCT matrix type. 2381 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 2382 Also, only cell-centered variables. 2383 2384 .seealso: PCMG, MATHYPRESSTRUCT 2385 M*/ 2386 2387 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2388 { 2389 PetscErrorCode ierr; 2390 PC_SysPFMG *ex; 2391 2392 PetscFunctionBegin; 2393 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2394 pc->data = ex; 2395 2396 ex->its = 1; 2397 ex->tol = 1.e-8; 2398 ex->relax_type = 1; 2399 ex->num_pre_relax = 1; 2400 ex->num_post_relax = 1; 2401 2402 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2403 pc->ops->view = PCView_SysPFMG; 2404 pc->ops->destroy = PCDestroy_SysPFMG; 2405 pc->ops->apply = PCApply_SysPFMG; 2406 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2407 pc->ops->setup = PCSetUp_SysPFMG; 2408 2409 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 2410 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2411 PetscFunctionReturn(0); 2412 } 2413