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