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