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