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 DMCreateInterpolation_DA(DM,DM,Mat*,Vec*); 123 extern PetscErrorCode DMCreateColoring_DA(DM,ISColoringType,MatType,ISColoring*); 124 extern PetscErrorCode DMCreateMatrix_DA(DM,MatType,Mat*); 125 extern PetscErrorCode DMCreateCoordinateDM_DA(DM,DM*); 126 extern PetscErrorCode DMRefine_DA(DM,MPI_Comm,DM*); 127 extern PetscErrorCode DMCoarsen_DA(DM,MPI_Comm,DM*); 128 extern PetscErrorCode DMRefineHierarchy_DA(DM,PetscInt,DM[]); 129 extern PetscErrorCode DMCoarsenHierarchy_DA(DM,PetscInt,DM[]); 130 extern PetscErrorCode DMCreateInjection_DA(DM,DM,VecScatter*); 131 extern PetscErrorCode DMCreateAggregates_DA(DM,DM,Mat*); 132 extern PetscErrorCode DMView_DA(DM,PetscViewer); 133 extern PetscErrorCode DMSetUp_DA(DM); 134 extern PetscErrorCode DMDestroy_DA(DM); 135 extern PetscErrorCode DMCreateDomainDecomposition_DA(DM,PetscInt*,char***,IS**,IS**,DM**); 136 extern PetscErrorCode DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter**,VecScatter**,VecScatter**); 137 138 #undef __FUNCT__ 139 #define __FUNCT__ "DMLoad_DA" 140 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer) 141 { 142 PetscErrorCode ierr; 143 PetscInt dim,m,n,p,dof,swidth; 144 DMDAStencilType stencil; 145 DMDABoundaryType bx,by,bz; 146 PetscBool coors; 147 DM dac; 148 Vec c; 149 150 PetscFunctionBegin; 151 ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 152 ierr = PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 153 ierr = PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 154 ierr = PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 155 ierr = PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 156 ierr = PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 157 ierr = PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 158 ierr = PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 159 ierr = PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 160 ierr = PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 161 162 ierr = DMDASetDim(da, dim);CHKERRQ(ierr); 163 ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr); 164 ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr); 165 ierr = DMDASetDof(da, dof);CHKERRQ(ierr); 166 ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr); 167 ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr); 168 ierr = DMSetUp(da);CHKERRQ(ierr); 169 ierr = PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);CHKERRQ(ierr); 170 if (coors) { 171 ierr = DMGetCoordinateDM(da,&dac);CHKERRQ(ierr); 172 ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr); 173 ierr = VecLoad(c,viewer);CHKERRQ(ierr); 174 ierr = DMSetCoordinates(da,c);CHKERRQ(ierr); 175 ierr = VecDestroy(&c);CHKERRQ(ierr); 176 } 177 PetscFunctionReturn(0); 178 } 179 180 #undef __FUNCT__ 181 #define __FUNCT__ "DMCreateSubDM_DA" 182 PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 183 { 184 PetscErrorCode ierr; 185 DM_DA *da = (DM_DA*)dm->data; 186 187 PetscFunctionBegin; 188 if (da->dim != 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Support only implemented for 2d"); 189 if (subdm) { 190 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); 191 } 192 if (is) { 193 PetscInt *indices,cnt = 0, dof = da->w,i,j; 194 ierr = PetscMalloc(sizeof(PetscInt)*da->Nlocal*numFields/dof,&indices);CHKERRQ(ierr); 195 for (i=da->base/dof; i<(da->base+da->Nlocal)/dof; i++) { 196 for (j=0; j<numFields; j++) { 197 indices[cnt++] = dof*i + fields[j]; 198 } 199 } 200 if (cnt != da->Nlocal*numFields/dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not equal expected value"); 201 ierr = ISCreateGeneral(PetscObjectComm((PetscObject)dm),da->Nlocal*numFields/dof,indices,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 202 } 203 PetscFunctionReturn(0); 204 } 205 206 #undef __FUNCT__ 207 #define __FUNCT__ "DMCreateFieldDecomposition_DA" 208 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM **dmlist) 209 { 210 PetscInt i; 211 PetscErrorCode ierr; 212 DM_DA *dd = (DM_DA*)dm->data; 213 PetscInt dof = dd->w; 214 215 PetscFunctionBegin; 216 if (len) *len = dof; 217 if (islist) { 218 Vec v; 219 PetscInt rstart,n; 220 221 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 222 ierr = VecGetOwnershipRange(v,&rstart,NULL);CHKERRQ(ierr); 223 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 224 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 225 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 226 for (i=0; i<dof; i++) { 227 ierr = ISCreateStride(PetscObjectComm((PetscObject)dm),n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 228 } 229 } 230 if (namelist) { 231 ierr = PetscMalloc(dof*sizeof(const char*), namelist);CHKERRQ(ierr); 232 if (dd->fieldname) { 233 for (i=0; i<dof; i++) { 234 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 235 } 236 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 237 } 238 if (dmlist) { 239 DM da; 240 241 ierr = DMDACreate(PetscObjectComm((PetscObject)dm), &da);CHKERRQ(ierr); 242 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 243 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 244 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 245 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 246 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 247 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 248 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 249 ierr = DMSetUp(da);CHKERRQ(ierr); 250 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 251 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 252 for (i=0; i<dof; i++) (*dmlist)[i] = da; 253 } 254 PetscFunctionReturn(0); 255 } 256 257 /*MC 258 DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions. 259 In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points. 260 In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width. 261 262 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 263 vertex centered. 264 265 266 Level: intermediate 267 268 .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType() 269 M*/ 270 271 272 #undef __FUNCT__ 273 #define __FUNCT__ "DMCreate_DA" 274 PETSC_EXTERN 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->Nsub = 1; 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 366 #undef __FUNCT__ 367 #define __FUNCT__ "DMDACreate" 368 /*@ 369 DMDACreate - Creates a DMDA object. 370 371 Collective on MPI_Comm 372 373 Input Parameter: 374 . comm - The communicator for the DMDA object 375 376 Output Parameter: 377 . da - The DMDA object 378 379 Level: advanced 380 381 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 382 383 .keywords: DMDA, create 384 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 385 @*/ 386 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 387 { 388 PetscErrorCode ierr; 389 390 PetscFunctionBegin; 391 PetscValidPointer(da,2); 392 ierr = DMCreate(comm,da);CHKERRQ(ierr); 393 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 394 PetscFunctionReturn(0); 395 } 396 397 398