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