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 /* duplicate communicator for hypre */ 1941 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr); 1942 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr); 1943 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr); 1944 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr); 1945 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr); 1946 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr); 1947 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);CHKERRQ(ierr); 1948 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);CHKERRQ(ierr); 1949 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);CHKERRQ(ierr); 1950 PetscFunctionReturn(0); 1951 } 1952 1953 /* ---------------------------------------------------------------------------------------------------------------------------------*/ 1954 1955 typedef struct { 1956 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */ 1957 HYPRE_StructSolver hsolver; 1958 1959 /* keep copy of PFMG options used so may view them */ 1960 PetscInt its; 1961 double tol; 1962 PetscInt relax_type; 1963 PetscInt rap_type; 1964 PetscInt num_pre_relax,num_post_relax; 1965 PetscInt max_levels; 1966 } PC_PFMG; 1967 1968 PetscErrorCode PCDestroy_PFMG(PC pc) 1969 { 1970 PetscErrorCode ierr; 1971 PC_PFMG *ex = (PC_PFMG*) pc->data; 1972 1973 PetscFunctionBegin; 1974 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 1975 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 1976 ierr = PetscFree(pc->data);CHKERRQ(ierr); 1977 PetscFunctionReturn(0); 1978 } 1979 1980 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"}; 1981 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"}; 1982 1983 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer) 1984 { 1985 PetscErrorCode ierr; 1986 PetscBool iascii; 1987 PC_PFMG *ex = (PC_PFMG*) pc->data; 1988 1989 PetscFunctionBegin; 1990 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1991 if (iascii) { 1992 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr); 1993 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 1994 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 1995 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 1996 ierr = PetscViewerASCIIPrintf(viewer," RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr); 1997 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 1998 ierr = PetscViewerASCIIPrintf(viewer," max levels %d\n",ex->max_levels);CHKERRQ(ierr); 1999 } 2000 PetscFunctionReturn(0); 2001 } 2002 2003 PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2004 { 2005 PetscErrorCode ierr; 2006 PC_PFMG *ex = (PC_PFMG*) pc->data; 2007 PetscBool flg = PETSC_FALSE; 2008 2009 PetscFunctionBegin; 2010 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr); 2011 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2012 if (flg) { 2013 int level=3; 2014 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level)); 2015 } 2016 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2017 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its)); 2018 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); 2019 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax)); 2020 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); 2021 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax)); 2022 2023 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr); 2024 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels)); 2025 2026 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2027 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol)); 2028 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); 2029 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type)); 2030 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr); 2031 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type)); 2032 ierr = PetscOptionsTail();CHKERRQ(ierr); 2033 PetscFunctionReturn(0); 2034 } 2035 2036 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y) 2037 { 2038 PetscErrorCode ierr; 2039 PC_PFMG *ex = (PC_PFMG*) pc->data; 2040 PetscScalar *yy; 2041 const PetscScalar *xx; 2042 PetscInt ilower[3],iupper[3]; 2043 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2044 2045 PetscFunctionBegin; 2046 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2047 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2048 iupper[0] += ilower[0] - 1; 2049 iupper[1] += ilower[1] - 1; 2050 iupper[2] += ilower[2] - 1; 2051 2052 /* copy x values over to hypre */ 2053 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0)); 2054 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2055 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx)); 2056 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2057 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb)); 2058 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2059 2060 /* copy solution values back to PETSc */ 2061 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2062 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy)); 2063 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 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) 2068 { 2069 PC_PFMG *jac = (PC_PFMG*)pc->data; 2070 PetscErrorCode ierr; 2071 PetscInt oits; 2072 2073 PetscFunctionBegin; 2074 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2075 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its)); 2076 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol)); 2077 2078 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr); 2079 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits)); 2080 *outits = oits; 2081 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2082 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2083 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol)); 2084 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its)); 2085 PetscFunctionReturn(0); 2086 } 2087 2088 2089 PetscErrorCode PCSetUp_PFMG(PC pc) 2090 { 2091 PetscErrorCode ierr; 2092 PC_PFMG *ex = (PC_PFMG*) pc->data; 2093 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data); 2094 PetscBool flg; 2095 2096 PetscFunctionBegin; 2097 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr); 2098 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner"); 2099 2100 /* create the hypre solver object and set its information */ 2101 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver)); 2102 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2103 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx)); 2104 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver)); 2105 PetscFunctionReturn(0); 2106 } 2107 2108 /*MC 2109 PCPFMG - the hypre PFMG multigrid solver 2110 2111 Level: advanced 2112 2113 Options Database: 2114 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner 2115 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2116 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2117 . -pc_pfmg_tol <tol> tolerance of PFMG 2118 . -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 2119 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin 2120 2121 Notes: This is for CELL-centered descretizations 2122 2123 This must be used with the MATHYPRESTRUCT matrix type. 2124 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA. 2125 2126 .seealso: PCMG, MATHYPRESTRUCT 2127 M*/ 2128 2129 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) 2130 { 2131 PetscErrorCode ierr; 2132 PC_PFMG *ex; 2133 2134 PetscFunctionBegin; 2135 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2136 pc->data = ex; 2137 2138 ex->its = 1; 2139 ex->tol = 1.e-8; 2140 ex->relax_type = 1; 2141 ex->rap_type = 0; 2142 ex->num_pre_relax = 1; 2143 ex->num_post_relax = 1; 2144 ex->max_levels = 0; 2145 2146 pc->ops->setfromoptions = PCSetFromOptions_PFMG; 2147 pc->ops->view = PCView_PFMG; 2148 pc->ops->destroy = PCDestroy_PFMG; 2149 pc->ops->apply = PCApply_PFMG; 2150 pc->ops->applyrichardson = PCApplyRichardson_PFMG; 2151 pc->ops->setup = PCSetUp_PFMG; 2152 2153 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 2154 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver)); 2155 PetscFunctionReturn(0); 2156 } 2157 2158 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/ 2159 2160 /* we know we are working with a HYPRE_SStructMatrix */ 2161 typedef struct { 2162 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */ 2163 HYPRE_SStructSolver ss_solver; 2164 2165 /* keep copy of SYSPFMG options used so may view them */ 2166 PetscInt its; 2167 double tol; 2168 PetscInt relax_type; 2169 PetscInt num_pre_relax,num_post_relax; 2170 } PC_SysPFMG; 2171 2172 PetscErrorCode PCDestroy_SysPFMG(PC pc) 2173 { 2174 PetscErrorCode ierr; 2175 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2176 2177 PetscFunctionBegin; 2178 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2179 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr); 2180 ierr = PetscFree(pc->data);CHKERRQ(ierr); 2181 PetscFunctionReturn(0); 2182 } 2183 2184 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"}; 2185 2186 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer) 2187 { 2188 PetscErrorCode ierr; 2189 PetscBool iascii; 2190 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2191 2192 PetscFunctionBegin; 2193 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 2194 if (iascii) { 2195 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr); 2196 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr); 2197 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr); 2198 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr); 2199 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr); 2200 } 2201 PetscFunctionReturn(0); 2202 } 2203 2204 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc) 2205 { 2206 PetscErrorCode ierr; 2207 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2208 PetscBool flg = PETSC_FALSE; 2209 2210 PetscFunctionBegin; 2211 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr); 2212 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr); 2213 if (flg) { 2214 int level=3; 2215 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level)); 2216 } 2217 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr); 2218 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its)); 2219 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); 2220 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax)); 2221 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); 2222 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax)); 2223 2224 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr); 2225 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol)); 2226 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); 2227 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type)); 2228 ierr = PetscOptionsTail();CHKERRQ(ierr); 2229 PetscFunctionReturn(0); 2230 } 2231 2232 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y) 2233 { 2234 PetscErrorCode ierr; 2235 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2236 PetscScalar *yy; 2237 const PetscScalar *xx; 2238 PetscInt ilower[3],iupper[3]; 2239 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2240 PetscInt ordering= mx->dofs_order; 2241 PetscInt nvars = mx->nvars; 2242 PetscInt part = 0; 2243 PetscInt size; 2244 PetscInt i; 2245 2246 PetscFunctionBegin; 2247 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2248 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr); 2249 iupper[0] += ilower[0] - 1; 2250 iupper[1] += ilower[1] - 1; 2251 iupper[2] += ilower[2] - 1; 2252 2253 size = 1; 2254 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1); 2255 2256 /* copy x values over to hypre for variable ordering */ 2257 if (ordering) { 2258 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2259 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2260 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i))); 2261 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2262 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2263 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x)); 2264 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2265 2266 /* copy solution values back to PETSc */ 2267 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2268 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i))); 2269 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2270 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */ 2271 PetscScalar *z; 2272 PetscInt j, k; 2273 2274 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr); 2275 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0)); 2276 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 2277 2278 /* transform nodal to hypre's variable ordering for sys_pfmg */ 2279 for (i= 0; i< size; i++) { 2280 k= i*nvars; 2281 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j]; 2282 } 2283 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2284 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 2285 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b)); 2286 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2287 2288 /* copy solution values back to PETSc */ 2289 ierr = VecGetArray(y,&yy);CHKERRQ(ierr); 2290 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i))); 2291 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */ 2292 for (i= 0; i< size; i++) { 2293 k= i*nvars; 2294 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i]; 2295 } 2296 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr); 2297 ierr = PetscFree(z);CHKERRQ(ierr); 2298 } 2299 PetscFunctionReturn(0); 2300 } 2301 2302 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) 2303 { 2304 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data; 2305 PetscErrorCode ierr; 2306 PetscInt oits; 2307 2308 PetscFunctionBegin; 2309 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr); 2310 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its)); 2311 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol)); 2312 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr); 2313 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits)); 2314 *outits = oits; 2315 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS; 2316 else *reason = PCRICHARDSON_CONVERGED_RTOL; 2317 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol)); 2318 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its)); 2319 PetscFunctionReturn(0); 2320 } 2321 2322 PetscErrorCode PCSetUp_SysPFMG(PC pc) 2323 { 2324 PetscErrorCode ierr; 2325 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data; 2326 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data); 2327 PetscBool flg; 2328 2329 PetscFunctionBegin; 2330 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr); 2331 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner"); 2332 2333 /* create the hypre sstruct solver object and set its information */ 2334 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver)); 2335 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2336 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver)); 2337 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x)); 2338 PetscFunctionReturn(0); 2339 } 2340 2341 /*MC 2342 PCSysPFMG - the hypre SysPFMG multigrid solver 2343 2344 Level: advanced 2345 2346 Options Database: 2347 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner 2348 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid 2349 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid 2350 . -pc_syspfmg_tol <tol> tolerance of SysPFMG 2351 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel 2352 2353 Notes: This is for CELL-centered descretizations 2354 2355 This must be used with the MATHYPRESSTRUCT matrix type. 2356 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA. 2357 Also, only cell-centered variables. 2358 2359 .seealso: PCMG, MATHYPRESSTRUCT 2360 M*/ 2361 2362 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) 2363 { 2364 PetscErrorCode ierr; 2365 PC_SysPFMG *ex; 2366 2367 PetscFunctionBegin; 2368 ierr = PetscNew(&ex);CHKERRQ(ierr); \ 2369 pc->data = ex; 2370 2371 ex->its = 1; 2372 ex->tol = 1.e-8; 2373 ex->relax_type = 1; 2374 ex->num_pre_relax = 1; 2375 ex->num_post_relax = 1; 2376 2377 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG; 2378 pc->ops->view = PCView_SysPFMG; 2379 pc->ops->destroy = PCDestroy_SysPFMG; 2380 pc->ops->apply = PCApply_SysPFMG; 2381 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG; 2382 pc->ops->setup = PCSetUp_SysPFMG; 2383 2384 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr); 2385 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver)); 2386 PetscFunctionReturn(0); 2387 } 2388