1 2 #include <petsc-private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "DMSetFromOptions_DA" 6 PetscErrorCode DMSetFromOptions_DA(DM da) 7 { 8 PetscErrorCode ierr; 9 DM_DA *dd = (DM_DA*)da->data; 10 PetscInt refine = 0,maxnlevels = 100,refx[100],refy[100],refz[100],n,i; 11 PetscBool negativeMNP = PETSC_FALSE,bM = PETSC_FALSE,bN = PETSC_FALSE, bP = PETSC_FALSE,flg; 12 13 PetscFunctionBegin; 14 PetscValidHeaderSpecific(da,DM_CLASSID,1); 15 16 if (dd->M < 0) { 17 dd->M = -dd->M; 18 bM = PETSC_TRUE; 19 negativeMNP = PETSC_TRUE; 20 } 21 if (dd->dim > 1 && dd->N < 0) { 22 dd->N = -dd->N; 23 bN = PETSC_TRUE; 24 negativeMNP = PETSC_TRUE; 25 } 26 if (dd->dim > 2 && dd->P < 0) { 27 dd->P = -dd->P; 28 bP = PETSC_TRUE; 29 negativeMNP = PETSC_TRUE; 30 } 31 32 ierr = PetscOptionsHead("DMDA Options");CHKERRQ(ierr); 33 if (bM) {ierr = PetscOptionsInt("-da_grid_x","Number of grid points in x direction","DMDASetSizes",dd->M,&dd->M,NULL);CHKERRQ(ierr);} 34 if (bN) {ierr = PetscOptionsInt("-da_grid_y","Number of grid points in y direction","DMDASetSizes",dd->N,&dd->N,NULL);CHKERRQ(ierr);} 35 if (bP) {ierr = PetscOptionsInt("-da_grid_z","Number of grid points in z direction","DMDASetSizes",dd->P,&dd->P,NULL);CHKERRQ(ierr);} 36 37 ierr = PetscOptionsInt("-da_overlap","Decomposition overlap in all directions","DMDASetOverlap",dd->xol,&dd->xol,&flg);CHKERRQ(ierr); 38 if (flg) {ierr = DMDASetOverlap(da,dd->xol,dd->xol,dd->xol);CHKERRQ(ierr);} 39 ierr = PetscOptionsInt("-da_overlap_x","Decomposition overlap in x direction","DMDASetOverlap",dd->xol,&dd->xol,NULL);CHKERRQ(ierr); 40 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_overlap_y","Decomposition overlap in y direction","DMDASetOverlap",dd->yol,&dd->yol,NULL);CHKERRQ(ierr);} 41 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_overlap_z","Decomposition overlap in z direction","DMDASetOverlap",dd->zol,&dd->zol,NULL);CHKERRQ(ierr);} 42 43 ierr = PetscOptionsInt("-da_local_subdomains","","DMDASetNumLocalSubdomains",dd->Nsub,&dd->Nsub,&flg);CHKERRQ(ierr); 44 if (flg) {ierr = DMDASetNumLocalSubDomains(da,dd->Nsub);CHKERRQ(ierr);} 45 46 /* Handle DMDA parallel distibution */ 47 ierr = PetscOptionsInt("-da_processors_x","Number of processors in x direction","DMDASetNumProcs",dd->m,&dd->m,NULL);CHKERRQ(ierr); 48 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_processors_y","Number of processors in y direction","DMDASetNumProcs",dd->n,&dd->n,NULL);CHKERRQ(ierr);} 49 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_processors_z","Number of processors in z direction","DMDASetNumProcs",dd->p,&dd->p,NULL);CHKERRQ(ierr);} 50 /* Handle DMDA refinement */ 51 ierr = PetscOptionsInt("-da_refine_x","Refinement ratio in x direction","DMDASetRefinementFactor",dd->refine_x,&dd->refine_x,NULL);CHKERRQ(ierr); 52 if (dd->dim > 1) {ierr = PetscOptionsInt("-da_refine_y","Refinement ratio in y direction","DMDASetRefinementFactor",dd->refine_y,&dd->refine_y,NULL);CHKERRQ(ierr);} 53 if (dd->dim > 2) {ierr = PetscOptionsInt("-da_refine_z","Refinement ratio in z direction","DMDASetRefinementFactor",dd->refine_z,&dd->refine_z,NULL);CHKERRQ(ierr);} 54 dd->coarsen_x = dd->refine_x; dd->coarsen_y = dd->refine_y; dd->coarsen_z = dd->refine_z; 55 56 /* Get refinement factors, defaults taken from the coarse DMDA */ 57 ierr = DMDAGetRefinementFactor(da,&refx[0],&refy[0],&refz[0]);CHKERRQ(ierr); 58 for (i=1; i<maxnlevels; i++) { 59 refx[i] = refx[0]; 60 refy[i] = refy[0]; 61 refz[i] = refz[0]; 62 } 63 n = maxnlevels; 64 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_x",refx,&n,NULL);CHKERRQ(ierr); 65 if (da->levelup - da->leveldown >= 0) dd->refine_x = refx[da->levelup - da->leveldown]; 66 if (da->levelup - da->leveldown >= 1) dd->coarsen_x = refx[da->levelup - da->leveldown - 1]; 67 if (dd->dim > 1) { 68 n = maxnlevels; 69 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_y",refy,&n,NULL);CHKERRQ(ierr); 70 if (da->levelup - da->leveldown >= 0) dd->refine_y = refy[da->levelup - da->leveldown]; 71 if (da->levelup - da->leveldown >= 1) dd->coarsen_y = refy[da->levelup - da->leveldown - 1]; 72 } 73 if (dd->dim > 2) { 74 n = maxnlevels; 75 ierr = PetscOptionsGetIntArray(((PetscObject)da)->prefix,"-da_refine_hierarchy_z",refz,&n,NULL);CHKERRQ(ierr); 76 if (da->levelup - da->leveldown >= 0) dd->refine_z = refz[da->levelup - da->leveldown]; 77 if (da->levelup - da->leveldown >= 1) dd->coarsen_z = refz[da->levelup - da->leveldown - 1]; 78 } 79 80 if (negativeMNP) {ierr = PetscOptionsInt("-da_refine","Uniformly refine DA one or more times","None",refine,&refine,NULL);CHKERRQ(ierr);} 81 ierr = PetscOptionsTail();CHKERRQ(ierr); 82 83 while (refine--) { 84 if (dd->bx == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) { 85 dd->M = dd->refine_x*dd->M; 86 } else { 87 dd->M = 1 + dd->refine_x*(dd->M - 1); 88 } 89 if (dd->by == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) { 90 dd->N = dd->refine_y*dd->N; 91 } else { 92 dd->N = 1 + dd->refine_y*(dd->N - 1); 93 } 94 if (dd->bz == DMDA_BOUNDARY_PERIODIC || dd->interptype == DMDA_Q0) { 95 dd->P = dd->refine_z*dd->P; 96 } else { 97 dd->P = 1 + dd->refine_z*(dd->P - 1); 98 } 99 da->levelup++; 100 if (da->levelup - da->leveldown >= 0) { 101 dd->refine_x = refx[da->levelup - da->leveldown]; 102 dd->refine_y = refy[da->levelup - da->leveldown]; 103 dd->refine_z = refz[da->levelup - da->leveldown]; 104 } 105 if (da->levelup - da->leveldown >= 1) { 106 dd->coarsen_x = refx[da->levelup - da->leveldown - 1]; 107 dd->coarsen_y = refy[da->levelup - da->leveldown - 1]; 108 dd->coarsen_z = refz[da->levelup - da->leveldown - 1]; 109 } 110 } 111 PetscFunctionReturn(0); 112 } 113 114 extern PetscErrorCode DMCreateGlobalVector_DA(DM,Vec*); 115 extern PetscErrorCode DMCreateLocalVector_DA(DM,Vec*); 116 extern PetscErrorCode DMGlobalToLocalBegin_DA(DM,Vec,InsertMode,Vec); 117 extern PetscErrorCode DMGlobalToLocalEnd_DA(DM,Vec,InsertMode,Vec); 118 extern PetscErrorCode DMLocalToGlobalBegin_DA(DM,Vec,InsertMode,Vec); 119 extern PetscErrorCode DMLocalToGlobalEnd_DA(DM,Vec,InsertMode,Vec); 120 extern PetscErrorCode DMCreateInterpolation_DA(DM,DM,Mat*,Vec*); 121 extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,MatType,ISColoring*); 122 extern PetscErrorCode DMCreateMatrix_DA(DM,MatType,Mat*); 123 extern PetscErrorCode DMCreateCoordinateDM_DA(DM,DM*); 124 extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*); 125 extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*); 126 extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]); 127 extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]); 128 extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*); 129 extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*); 130 extern PetscErrorCode DMView_DA(DM,PetscViewer); 131 extern PetscErrorCode DMSetUp_DA(DM); 132 extern PetscErrorCode DMDestroy_DA(DM); 133 extern PetscErrorCode DMCreateDomainDecomposition_DA(DM,PetscInt*,char***,IS**,IS**,DM**); 134 extern PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**); 135 136 #undef __FUNCT__ 137 #define __FUNCT__ "DMLoad_DA" 138 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer) 139 { 140 PetscErrorCode ierr; 141 PetscInt dim,m,n,p,dof,swidth; 142 DMDAStencilType stencil; 143 DMDABoundaryType bx,by,bz; 144 PetscBool coors; 145 DM dac; 146 Vec c; 147 148 PetscFunctionBegin; 149 ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 150 ierr = PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 151 ierr = PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 152 ierr = PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 153 ierr = PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 154 ierr = PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 155 ierr = PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 156 ierr = PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 157 ierr = PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 158 ierr = PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 159 160 ierr = DMDASetDim(da, dim);CHKERRQ(ierr); 161 ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr); 162 ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr); 163 ierr = DMDASetDof(da, dof);CHKERRQ(ierr); 164 ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr); 165 ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr); 166 ierr = DMSetUp(da);CHKERRQ(ierr); 167 ierr = PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);CHKERRQ(ierr); 168 if (coors) { 169 ierr = DMGetCoordinateDM(da,&dac);CHKERRQ(ierr); 170 ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr); 171 ierr = VecLoad(c,viewer);CHKERRQ(ierr); 172 ierr = DMSetCoordinates(da,c);CHKERRQ(ierr); 173 ierr = VecDestroy(&c);CHKERRQ(ierr); 174 } 175 PetscFunctionReturn(0); 176 } 177 178 #undef __FUNCT__ 179 #define __FUNCT__ "DMCreateSubDM_DA" 180 PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 181 { 182 PetscErrorCode ierr; 183 DM_DA *da = (DM_DA*)dm->data; 184 185 PetscFunctionBegin; 186 if (da->dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only implemented for 2d"); 187 if (subdm) { 188 ierr = DMDACreate2d(PetscObjectComm((PetscObject)dm),da->bx,da->by,da->stencil_type,da->M,da->N,da->m,da->n,numFields,da->s,da->lx,da->ly,subdm);CHKERRQ(ierr); 189 } 190 if (is) { 191 PetscInt *indices,cnt = 0, dof = da->w,i,j; 192 ierr = PetscMalloc(sizeof(PetscInt)*da->Nlocal*numFields/dof,&indices);CHKERRQ(ierr); 193 for (i=da->base/dof; i<(da->base+da->Nlocal)/dof; i++) { 194 for (j=0; j<numFields; j++) { 195 indices[cnt++] = dof*i + fields[j]; 196 } 197 } 198 if (cnt != da->Nlocal*numFields/dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not equal expected value"); 199 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),da->Nlocal*numFields/dof,indices,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 200 } 201 PetscFunctionReturn(0); 202 } 203 204 #undef __FUNCT__ 205 #define __FUNCT__ "DMCreateFieldDecomposition_DA" 206 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist) 207 { 208 PetscInt i; 209 PetscErrorCode ierr; 210 DM_DA *dd = (DM_DA*)dm->data; 211 PetscInt dof = dd->w; 212 213 PetscFunctionBegin; 214 if (len) *len = dof; 215 if (islist) { 216 Vec v; 217 PetscInt rstart,n; 218 219 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 220 ierr = VecGetOwnershipRange(v,&rstart,NULL);CHKERRQ(ierr); 221 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 222 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 223 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 224 for (i=0; i<dof; i++) { 225 ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 226 } 227 } 228 if (namelist) { 229 ierr = PetscMalloc(dof*sizeof(const char*), namelist);CHKERRQ(ierr); 230 if (dd->fieldname) { 231 for (i=0; i<dof; i++) { 232 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 233 } 234 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 235 } 236 if (dmlist) { 237 DM da; 238 239 ierr = DMDACreate(PetscObjectComm((PetscObject)dm), &da);CHKERRQ(ierr); 240 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 241 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 242 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 243 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 244 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 245 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 246 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 247 ierr = DMSetUp(da);CHKERRQ(ierr); 248 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 249 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 250 for (i=0; i<dof; i++) (*dmlist)[i] = da; 251 } 252 PetscFunctionReturn(0); 253 } 254 255 /*MC 256 DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions. 257 In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points. 258 In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width. 259 260 The vectors can be thought of as either cell centered or vertex centered on the mesh. But some variables cannot be cell centered and others 261 vertex centered. 262 263 264 Level: intermediate 265 266 .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType() 267 M*/ 268 269 270 #undef __FUNCT__ 271 #define __FUNCT__ "DMCreate_DA" 272 PETSC_EXTERN PetscErrorCode DMCreate_DA(DM da) 273 { 274 PetscErrorCode ierr; 275 DM_DA *dd; 276 277 PetscFunctionBegin; 278 PetscValidPointer(da,1); 279 ierr = PetscNewLog(da,DM_DA,&dd);CHKERRQ(ierr); 280 da->data = dd; 281 282 dd->dim = -1; 283 dd->interptype = DMDA_Q1; 284 dd->refine_x = 2; 285 dd->refine_y = 2; 286 dd->refine_z = 2; 287 dd->coarsen_x = 2; 288 dd->coarsen_y = 2; 289 dd->coarsen_z = 2; 290 dd->fieldname = NULL; 291 dd->nlocal = -1; 292 dd->Nlocal = -1; 293 dd->M = -1; 294 dd->N = -1; 295 dd->P = -1; 296 dd->m = -1; 297 dd->n = -1; 298 dd->p = -1; 299 dd->w = -1; 300 dd->s = -1; 301 302 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 303 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 304 305 dd->Nsub = 1; 306 dd->xol = 0; 307 dd->yol = 0; 308 dd->zol = 0; 309 dd->xo = 0; 310 dd->yo = 0; 311 dd->zo = 0; 312 dd->Mo = -1; 313 dd->No = -1; 314 dd->Po = -1; 315 316 dd->gtol = NULL; 317 dd->ltog = NULL; 318 dd->ltol = NULL; 319 dd->ao = NULL; 320 dd->base = -1; 321 dd->bx = DMDA_BOUNDARY_NONE; 322 dd->by = DMDA_BOUNDARY_NONE; 323 dd->bz = DMDA_BOUNDARY_NONE; 324 dd->stencil_type = DMDA_STENCIL_BOX; 325 dd->interptype = DMDA_Q1; 326 dd->idx = NULL; 327 dd->Nl = -1; 328 dd->lx = NULL; 329 dd->ly = NULL; 330 dd->lz = NULL; 331 332 dd->elementtype = DMDA_ELEMENT_Q1; 333 334 ierr = PetscStrallocpy(VECSTANDARD,(char**)&da->vectype);CHKERRQ(ierr); 335 336 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 337 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 338 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 339 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 340 da->ops->createglobalvector = DMCreateGlobalVector_DA; 341 da->ops->createlocalvector = DMCreateLocalVector_DA; 342 da->ops->createinterpolation = DMCreateInterpolation_DA; 343 da->ops->getcoloring = DMCreateColoring_DA; 344 da->ops->creatematrix = DMCreateMatrix_DA; 345 da->ops->refine = DMRefine_DA; 346 da->ops->coarsen = DMCoarsen_DA; 347 da->ops->refinehierarchy = DMRefineHierarchy_DA; 348 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 349 da->ops->getinjection = DMCreateInjection_DA; 350 da->ops->getaggregates = DMCreateAggregates_DA; 351 da->ops->destroy = DMDestroy_DA; 352 da->ops->view = 0; 353 da->ops->setfromoptions = DMSetFromOptions_DA; 354 da->ops->setup = DMSetUp_DA; 355 da->ops->load = DMLoad_DA; 356 da->ops->createcoordinatedm = DMCreateCoordinateDM_DA; 357 da->ops->createsubdm = DMCreateSubDM_DA; 358 da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA; 359 da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA; 360 da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA; 361 PetscFunctionReturn(0); 362 } 363 364 #undef __FUNCT__ 365 #define __FUNCT__ "DMDACreate" 366 /*@ 367 DMDACreate - Creates a DMDA object. 368 369 Collective on MPI_Comm 370 371 Input Parameter: 372 . comm - The communicator for the DMDA object 373 374 Output Parameter: 375 . da - The DMDA object 376 377 Level: advanced 378 379 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 380 381 .keywords: DMDA, create 382 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 383 @*/ 384 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 385 { 386 PetscErrorCode ierr; 387 388 PetscFunctionBegin; 389 PetscValidPointer(da,2); 390 ierr = DMCreate(comm,da);CHKERRQ(ierr); 391 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 392 PetscFunctionReturn(0); 393 } 394 395 396