1 2 #include <petscdm.h> 3 #include <petscctable.h> 4 #include <petsc/private/matimpl.h> 5 /* Need to access the hypre private data */ 6 #include <_hypre_parcsr_ls.h> 7 #include <petsc/private/pcmgimpl.h> 8 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 9 10 typedef struct { 11 PC hypre; 12 PetscBool reuseinterp; 13 PetscInt blocksize; 14 } PC_HMG; 15 16 PetscErrorCode MatConvert_ParCSRMatrix_AIJ(MPI_Comm, hypre_ParCSRMatrix*, MatType, MatReuse, Mat*); 17 PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems*,PC); 18 PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems*,PC); 19 PetscErrorCode PCReset_MG(PC); 20 PetscErrorCode PCHYPREGetSolver(PC,HYPRE_Solver*); 21 22 static PetscErrorCode PCHMGExtractSubMatrix_HMG(Mat pmat,Mat *submat,MatReuse reuse,PetscInt blocksize) 23 { 24 IS isrow; 25 PetscErrorCode ierr; 26 PetscInt rstart,rend, row,subsize, *rowsindices; 27 MPI_Comm comm; 28 29 PetscFunctionBegin; 30 ierr = PetscObjectGetComm((PetscObject)pmat,&comm);CHKERRQ(ierr); 31 ierr = MatGetOwnershipRange(pmat,&rstart,&rend);CHKERRQ(ierr); 32 if ((rend-rstart)%blocksize != 0) SETERRQ3(comm,PETSC_ERR_ARG_INCOMP,"Block size %d is inconsisent for [%d, %d) \n",blocksize,rstart,rend); 33 subsize = (rend-rstart)/blocksize; 34 ierr = PetscCalloc1(subsize,&rowsindices);CHKERRQ(ierr); 35 subsize = 0; 36 for (row=rstart; row<rend; row+=blocksize){ 37 rowsindices[subsize++]=row; 38 } 39 ierr = ISCreateGeneral(comm,subsize,rowsindices,PETSC_OWN_POINTER,&isrow);CHKERRQ(ierr); 40 ierr = MatCreateSubMatrix(pmat,isrow,isrow,reuse,submat);CHKERRQ(ierr); 41 ierr = ISDestroy(&isrow);CHKERRQ(ierr); 42 PetscFunctionReturn(0); 43 } 44 45 static PetscErrorCode PCHMGExpandInterpolation_HMG(MPI_Comm comm,hypre_ParCSRMatrix *subinterp, Mat *interp, PetscInt blocksize) 46 { 47 PetscInt subrstart,subrend,subrowsize,subcolsize,subcstart,subcend,rowsize,colsize; 48 PetscInt subrow,row,*idx,nz,*d_nnz,*o_nnz,i,j,dnz,onz,max_nz,*indices; 49 PetscScalar *values; 50 PetscErrorCode ierr; 51 52 PetscFunctionBegin; 53 subrstart = hypre_ParCSRMatrixFirstRowIndex(subinterp); 54 subrend = hypre_ParCSRMatrixLastRowIndex(subinterp); 55 subrowsize = subrend-subrstart+1; 56 rowsize = subrowsize*blocksize; 57 subcstart = hypre_ParCSRMatrixFirstColDiag(subinterp); 58 subcend = hypre_ParCSRMatrixLastColDiag(subinterp); 59 subcolsize = subcend-subcstart+1; 60 colsize = subcolsize*blocksize; 61 ierr = PetscCalloc2(rowsize,&d_nnz,rowsize,&o_nnz);CHKERRQ(ierr); 62 max_nz = 0; 63 for(subrow=subrstart; subrow<=subrend; subrow++){ 64 HYPRE_ParCSRMatrixGetRow(subinterp,subrow,(HYPRE_Int*)&nz,(HYPRE_Int**)&idx,NULL); 65 if (max_nz<nz) max_nz = nz; 66 dnz = 0; onz = 0; 67 for(i=0;i<nz;i++){ 68 if(idx[i]<subcstart || idx[i]>subcend) onz++; 69 else dnz++; 70 } 71 for(i=0;i<blocksize;i++){ 72 d_nnz[(subrow-subrstart)*blocksize+i] = dnz; 73 o_nnz[(subrow-subrstart)*blocksize+i] = onz; 74 } 75 HYPRE_ParCSRMatrixRestoreRow(subinterp,subrow,(HYPRE_Int*)&nz,(HYPRE_Int**)&idx,NULL); 76 } 77 ierr = MatCreateAIJ(comm,rowsize,colsize,PETSC_DETERMINE,PETSC_DETERMINE,0,d_nnz,0,o_nnz,interp);CHKERRQ(ierr); 78 ierr = MatSetOption(*interp,MAT_IGNORE_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 79 ierr = MatSetOption(*interp,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 80 ierr = MatSetOption(*interp,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 81 ierr = MatSetFromOptions(*interp);CHKERRQ(ierr); 82 83 ierr = MatSetUp(*interp);CHKERRQ(ierr); 84 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 85 ierr = PetscCalloc1(max_nz,&indices);CHKERRQ(ierr); 86 for(subrow=subrstart; subrow<=subrend; subrow++){ 87 HYPRE_ParCSRMatrixGetRow(subinterp,subrow,(HYPRE_Int*)&nz,(HYPRE_Int**)&idx,&values); 88 for(i=0;i<blocksize;i++){ 89 row = subrow*blocksize+i; 90 for (j=0;j<nz;j++){ 91 indices[j] = idx[j]*blocksize+i; 92 } 93 ierr = MatSetValues(*interp,1,&row,nz,indices,values,INSERT_VALUES);CHKERRQ(ierr); 94 } 95 HYPRE_ParCSRMatrixRestoreRow(subinterp,subrow,(HYPRE_Int*)&nz,(HYPRE_Int**)&idx,&values); 96 } 97 ierr = PetscFree(indices);CHKERRQ(ierr); 98 ierr = MatAssemblyBegin(*interp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 99 ierr = MatAssemblyEnd(*interp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 PetscErrorCode PCSetUp_HMG(PC pc) 104 { 105 PetscErrorCode ierr; 106 Mat PA, submat; 107 PC_MG *mg = (PC_MG*)pc->data; 108 PC_HMG *hmg = (PC_HMG*) mg->innerctx; 109 MPI_Comm comm; 110 PetscInt level; 111 hypre_ParCSRMatrix **P_array, **A_array; 112 HYPRE_Solver hsolver; 113 hypre_ParAMGData *amg_data; 114 PetscInt num_levels; 115 PetscReal global_nonzeros, num_rows; 116 PetscReal *sparse; 117 118 PetscFunctionBegin; 119 120 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 121 122 if(pc->setupcalled){ 123 if (pc->flag == SAME_NONZERO_PATTERN && hmg->reuseinterp) { 124 ierr = PCMGSetGalerkin(pc,PC_MG_GALERKIN_PMAT);CHKERRQ(ierr); 125 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 }else { 128 ierr = PCReset_MG(pc);CHKERRQ(ierr); 129 pc->setupcalled = PETSC_FALSE; 130 } 131 } 132 133 if (!hmg->hypre){ 134 ierr = PCCreate(comm,&hmg->hypre);CHKERRQ(ierr); 135 ierr = PCSetType(hmg->hypre,PCHYPRE);CHKERRQ(ierr); 136 ierr = PCHYPRESetType(hmg->hypre,"boomeramg");CHKERRQ(ierr); 137 } 138 ierr = PCGetOperators(pc,NULL,&PA);CHKERRQ(ierr); 139 if(hmg->blocksize>1) { 140 ierr = PCHMGExtractSubMatrix_HMG(PA,&submat,MAT_INITIAL_MATRIX,hmg->blocksize);CHKERRQ(ierr); 141 PA = submat; 142 } 143 ierr = PCSetOperators(hmg->hypre,PA,PA);CHKERRQ(ierr); 144 if (hmg->blocksize>1){ 145 ierr = MatDestroy(&PA);CHKERRQ(ierr); 146 } 147 ierr = PCSetUseAmat(hmg->hypre,PETSC_FALSE);CHKERRQ(ierr); 148 ierr = PetscObjectOptionsBegin((PetscObject)hmg->hypre);CHKERRQ(ierr); 149 ierr = PCSetFromOptions_HYPRE(PetscOptionsObject,hmg->hypre);CHKERRQ(ierr); 150 ierr = PetscOptionsEnd();CHKERRQ(ierr); 151 ierr = PCSetUp(hmg->hypre); 152 153 ierr = PCHYPREGetSolver(hmg->hypre,&hsolver);CHKERRQ(ierr); 154 amg_data = (hypre_ParAMGData*) (hsolver); 155 156 if (!amg_data) SETERRQ(comm,PETSC_ERR_ARG_BADPTR,"Fails to setup Hypre \n"); 157 158 num_levels = hypre_ParAMGDataNumLevels(amg_data); 159 ierr = PetscCalloc1(num_levels,&sparse);CHKERRQ(ierr); 160 161 A_array= hypre_ParAMGDataAArray(amg_data); 162 P_array = hypre_ParAMGDataPArray(amg_data); 163 for(level=0;level<num_levels;level++){ 164 global_nonzeros = (PetscReal) hypre_ParCSRMatrixDNumNonzeros(A_array[level]); 165 num_rows = (PetscReal) hypre_ParCSRMatrixGlobalNumRows(A_array[level]); 166 sparse[num_levels-1-level] =global_nonzeros/(num_rows*num_rows); 167 } 168 ierr = PCMGSetLevels_MG(pc,num_levels,NULL);CHKERRQ(ierr); 169 ierr = PetscFree(sparse);CHKERRQ(ierr); 170 171 for(level=num_levels-1;level>0;level--){ 172 Mat P=0, pmat=0; 173 Vec b, x,r; 174 if (hmg->blocksize>1){ 175 ierr = PCHMGExpandInterpolation_HMG(comm,P_array[num_levels-1-level],&P,hmg->blocksize);CHKERRQ(ierr); 176 }else{ 177 ierr = MatConvert_ParCSRMatrix_AIJ(comm, P_array[num_levels-1-level],MATAIJ, MAT_INPLACE_MATRIX,&P);CHKERRQ(ierr); 178 } 179 /*ierr = MatView(P,NULL);CHKERRQ(ierr);*/ 180 ierr = MatCreateVecs(P,&b,&r);CHKERRQ(ierr); 181 ierr = PCMGSetInterpolation(pc,level,P);CHKERRQ(ierr); 182 ierr = PCMGSetRestriction(pc,level,P);CHKERRQ(ierr); 183 ierr = MatDestroy(&P);CHKERRQ(ierr); 184 if ((level-1)>=0 && hmg->blocksize<=1) { 185 ierr = MatConvert_ParCSRMatrix_AIJ(comm, A_array[num_levels-level],MATAIJ, MAT_INPLACE_MATRIX,&pmat);CHKERRQ(ierr); 186 ierr = PCMGSetOperators(pc,level-1,pmat);CHKERRQ(ierr); 187 ierr = MatDestroy(&pmat);CHKERRQ(ierr); 188 } 189 ierr = PCMGSetRhs(pc,level-1,b);CHKERRQ(ierr); 190 191 ierr = PCMGSetR(pc,level,r);CHKERRQ(ierr); 192 ierr = VecDestroy(&r);CHKERRQ(ierr); 193 194 ierr = VecDuplicate(b,&x); 195 ierr = PCMGSetX(pc,level-1,x);CHKERRQ(ierr); 196 ierr = VecDestroy(&x);CHKERRQ(ierr); 197 ierr = VecDestroy(&b);CHKERRQ(ierr); 198 } 199 ierr = PCDestroy(&hmg->hypre);CHKERRQ(ierr); 200 hmg->hypre = 0; 201 ierr = PCMGSetGalerkin(pc,(hmg->blocksize>1) ? PC_MG_GALERKIN_PMAT:PC_MG_GALERKIN_NONE);CHKERRQ(ierr); 202 ierr = PCSetDM(pc,NULL);CHKERRQ(ierr); 203 ierr = PCSetUseAmat(pc,PETSC_FALSE);CHKERRQ(ierr); 204 ierr = PetscObjectOptionsBegin((PetscObject)pc);CHKERRQ(ierr); 205 ierr = PCSetFromOptions_MG(PetscOptionsObject,pc);CHKERRQ(ierr); /* should be called in PCSetFromOptions_HMG(), but cannot be called prior to PCMGSetLevels() */ 206 ierr = PetscOptionsEnd();CHKERRQ(ierr); 207 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 211 PetscErrorCode PCDestroy_HMG(PC pc) 212 { 213 PetscErrorCode ierr; 214 PC_MG *mg = (PC_MG*)pc->data; 215 PC_HMG *ctx = (PC_HMG*) mg->innerctx; 216 217 PetscFunctionBegin; 218 ierr = PCDestroy(&ctx->hypre);CHKERRQ(ierr); 219 ierr = PetscFree(ctx);CHKERRQ(ierr); 220 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 221 PetscFunctionReturn(0); 222 } 223 224 PetscErrorCode PCView_HMG(PC pc,PetscViewer viewer) 225 { 226 PC_MG *mg = (PC_MG*)pc->data; 227 PC_HMG *ctx = (PC_HMG*) mg->innerctx; 228 PetscErrorCode ierr; 229 PetscBool iascii; 230 231 PetscFunctionBegin; 232 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 233 if (iascii) { 234 ierr = PetscViewerASCIIPrintf(viewer," Reuse interpolation %d\n",ctx->reuseinterp);CHKERRQ(ierr); 235 ierr = PetscViewerASCIIPrintf(viewer," Matrix block size %D \n",ctx->blocksize);CHKERRQ(ierr); 236 } 237 ierr = PCView_MG(pc,viewer);CHKERRQ(ierr); 238 PetscFunctionReturn(0); 239 } 240 241 PetscErrorCode PCSetFromOptions_HMG(PetscOptionItems *PetscOptionsObject,PC pc) 242 { 243 PetscErrorCode ierr; 244 PC_MG *mg = (PC_MG*)pc->data; 245 PC_HMG *ctx = (PC_HMG*) mg->innerctx; 246 247 PetscFunctionBegin; 248 ierr = PetscOptionsHead(PetscOptionsObject,"HMG");CHKERRQ(ierr); 249 ierr = PetscOptionsBool("-pc_hmg_reuse_interpolation","Reuse the interpolation operators when possible (cheaper, weaker when matrix entries change a lot)","None",ctx->reuseinterp,&ctx->reuseinterp,NULL);CHKERRQ(ierr); 250 ierr = PetscOptionsInt("-pc_hmg_pmat_blocksize","Block size for each grid point","hmg",ctx->blocksize,&ctx->blocksize,NULL);CHKERRQ(ierr); 251 ierr = PetscOptionsTail();CHKERRQ(ierr); 252 PetscFunctionReturn(0); 253 } 254 255 /*MC 256 PCHMG - Hybrid of PETSc preconditioners (such as ASM, BJacobi, SOR, etc.) and Hypre BoomerAMG. BoomerAMG is used to coarsen matrix and generate 257 a sequence of coarse matrices and interpolations. The matrices and interpolations are employed to construct PCMG, and then any available 258 PETSc preconditioners can be chosen as smoothers and the coarse solver. 259 260 Options Database Keys: 261 + -pc_hmg_reuse_interpolation <true | false> - Whether or not or not to reuse the interpolations. If true, it potentially save the compute time. 262 . -pc_hmg_pmat_blocksize - Block size of the underlying matrix. 263 264 265 Notes: 266 For multicomponent problems, we can just coarsen one submatrix associated with one particular component. In this way, the preconditioner setup 267 time is significantly reduced. One typical use case is neutron transport equations. There are many variables on each mesh vertex due to the 268 of angle and energy. Each variable, in fact, corresponds to the same PDEs but with different material properties. 269 270 Level: beginner 271 272 Concepts: additive Schwarz method 273 274 References: 275 + 1. - Fande Kong, Yaqi Wang, Derek R Gaston, Cody J Permann, Andrew E Slaughter, Alexander D Lindsay, Richard C Martineau, A highly parallel multilevel 276 Newton-Krylov-Schwarz method with subspace-based coarsening and partition-based balancing for the multigroup neutron transport equations on 277 3D unstructured meshes, arXiv preprint arXiv:1903.03659, 2019 278 279 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMG, PCHYPRE 280 281 M*/ 282 PETSC_EXTERN PetscErrorCode PCCreate_HMG(PC pc) 283 { 284 PetscErrorCode ierr; 285 PC_HMG *hmg; 286 PC_MG *mg; 287 288 PetscFunctionBegin; 289 /* if type was previously mg; must manually destroy it because call to PCSetType(pc,PCMG) will not destroy it */ 290 if (pc->ops->destroy) { 291 ierr = (*pc->ops->destroy)(pc);CHKERRQ(ierr); 292 pc->data = 0; 293 } 294 ierr = PetscFree(((PetscObject)pc)->type_name);CHKERRQ(ierr); 295 ((PetscObject)pc)->type_name = 0; 296 297 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); 298 ierr = PetscNew(&hmg);CHKERRQ(ierr); \ 299 300 mg = (PC_MG*) pc->data; 301 mg->innerctx = hmg; 302 hmg->reuseinterp = PETSC_FALSE; 303 hmg->blocksize = 1; 304 305 pc->ops->setfromoptions = PCSetFromOptions_HMG; 306 pc->ops->view = PCView_HMG; 307 pc->ops->destroy = PCDestroy_HMG; 308 pc->ops->setup = PCSetUp_HMG; 309 PetscFunctionReturn(0); 310 } 311