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