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