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 DMCreateDomainDecompositionScatters_DA(DM,PetscInt,DM*,VecScatter **,VecScatter**,VecScatter**); 128 129 #undef __FUNCT__ 130 #define __FUNCT__ "DMLoad_DA" 131 PetscErrorCode DMLoad_DA(DM da,PetscViewer viewer) 132 { 133 PetscErrorCode ierr; 134 PetscInt dim,m,n,p,dof,swidth; 135 DMDAStencilType stencil; 136 DMDABoundaryType bx,by,bz; 137 PetscBool coors; 138 DM dac; 139 Vec c; 140 141 PetscFunctionBegin; 142 ierr = PetscViewerBinaryRead(viewer,&dim,1,PETSC_INT);CHKERRQ(ierr); 143 ierr = PetscViewerBinaryRead(viewer,&m,1,PETSC_INT);CHKERRQ(ierr); 144 ierr = PetscViewerBinaryRead(viewer,&n,1,PETSC_INT);CHKERRQ(ierr); 145 ierr = PetscViewerBinaryRead(viewer,&p,1,PETSC_INT);CHKERRQ(ierr); 146 ierr = PetscViewerBinaryRead(viewer,&dof,1,PETSC_INT);CHKERRQ(ierr); 147 ierr = PetscViewerBinaryRead(viewer,&swidth,1,PETSC_INT);CHKERRQ(ierr); 148 ierr = PetscViewerBinaryRead(viewer,&bx,1,PETSC_ENUM);CHKERRQ(ierr); 149 ierr = PetscViewerBinaryRead(viewer,&by,1,PETSC_ENUM);CHKERRQ(ierr); 150 ierr = PetscViewerBinaryRead(viewer,&bz,1,PETSC_ENUM);CHKERRQ(ierr); 151 ierr = PetscViewerBinaryRead(viewer,&stencil,1,PETSC_ENUM);CHKERRQ(ierr); 152 153 ierr = DMDASetDim(da, dim);CHKERRQ(ierr); 154 ierr = DMDASetSizes(da, m,n,p);CHKERRQ(ierr); 155 ierr = DMDASetBoundaryType(da, bx, by, bz);CHKERRQ(ierr); 156 ierr = DMDASetDof(da, dof);CHKERRQ(ierr); 157 ierr = DMDASetStencilType(da, stencil);CHKERRQ(ierr); 158 ierr = DMDASetStencilWidth(da, swidth);CHKERRQ(ierr); 159 ierr = DMSetUp(da);CHKERRQ(ierr); 160 ierr = PetscViewerBinaryRead(viewer,&coors,1,PETSC_ENUM);CHKERRQ(ierr); 161 if (coors) { 162 ierr = DMGetCoordinateDM(da,&dac);CHKERRQ(ierr); 163 ierr = DMCreateGlobalVector(dac,&c);CHKERRQ(ierr); 164 ierr = VecLoad(c,viewer);CHKERRQ(ierr); 165 ierr = DMSetCoordinates(da,c);CHKERRQ(ierr); 166 ierr = VecDestroy(&c);CHKERRQ(ierr); 167 } 168 PetscFunctionReturn(0); 169 } 170 171 #undef __FUNCT__ 172 #define __FUNCT__ "DMCreateSubDM_DA" 173 PetscErrorCode DMCreateSubDM_DA(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm) 174 { 175 PetscErrorCode ierr; 176 DM_DA *da = (DM_DA*)dm->data; 177 178 PetscFunctionBegin; 179 if (da->dim != 2) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"Support only implemented for 2d"); 180 if (subdm) { 181 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); 182 } 183 if (is) { 184 PetscInt *indices,cnt = 0, dof = da->w,i,j; 185 ierr = PetscMalloc(sizeof(PetscInt)*da->Nlocal*numFields/dof,&indices);CHKERRQ(ierr); 186 for (i=da->base/dof; i<(da->base+da->Nlocal)/dof; i++) { 187 for (j=0; j<numFields; j++) { 188 indices[cnt++] = dof*i + fields[j]; 189 } 190 } 191 if (cnt != da->Nlocal*numFields/dof) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Count does not equal expected value"); 192 ierr = ISCreateGeneral(((PetscObject)dm)->comm,da->Nlocal*numFields/dof,indices,PETSC_OWN_POINTER,is);CHKERRQ(ierr); 193 } 194 PetscFunctionReturn(0); 195 } 196 197 #undef __FUNCT__ 198 #define __FUNCT__ "DMCreateFieldDecomposition_DA" 199 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist) 200 { 201 PetscInt i; 202 PetscErrorCode ierr; 203 DM_DA *dd = (DM_DA*)dm->data; 204 PetscInt dof = dd->w; 205 206 PetscFunctionBegin; 207 if (len) *len = dof; 208 if (islist) { 209 Vec v; 210 PetscInt rstart,n; 211 212 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 213 ierr = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr); 214 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 215 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 216 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 217 for (i=0; i<dof; i++) { 218 ierr = ISCreateStride(((PetscObject)dm)->comm,n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 219 } 220 } 221 if (namelist) { 222 ierr = PetscMalloc(dof*sizeof(const char *), namelist);CHKERRQ(ierr); 223 if (dd->fieldname) { 224 for (i=0; i<dof; i++) { 225 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 226 } 227 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 228 } 229 if (dmlist) { 230 DM da; 231 232 ierr = DMDACreate(((PetscObject)dm)->comm, &da);CHKERRQ(ierr); 233 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 234 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 235 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 236 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 237 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 238 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 239 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 240 ierr = DMSetUp(da);CHKERRQ(ierr); 241 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 242 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 243 for (i=0; i<dof; i++) (*dmlist)[i] = da; 244 } 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 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 297 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 298 299 dd->overlap = 0; 300 dd->xo = 0; 301 dd->yo = 0; 302 dd->zo = 0; 303 dd->Mo = -1; 304 dd->No = -1; 305 dd->Po = -1; 306 307 dd->gtol = PETSC_NULL; 308 dd->ltog = PETSC_NULL; 309 dd->ltol = PETSC_NULL; 310 dd->ao = PETSC_NULL; 311 dd->base = -1; 312 dd->bx = DMDA_BOUNDARY_NONE; 313 dd->by = DMDA_BOUNDARY_NONE; 314 dd->bz = DMDA_BOUNDARY_NONE; 315 dd->stencil_type = DMDA_STENCIL_BOX; 316 dd->interptype = DMDA_Q1; 317 dd->idx = PETSC_NULL; 318 dd->Nl = -1; 319 dd->lx = PETSC_NULL; 320 dd->ly = PETSC_NULL; 321 dd->lz = PETSC_NULL; 322 323 dd->elementtype = DMDA_ELEMENT_Q1; 324 325 ierr = PetscStrallocpy(VECSTANDARD,(char**)&da->vectype);CHKERRQ(ierr); 326 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 327 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 328 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 329 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 330 da->ops->createglobalvector = DMCreateGlobalVector_DA; 331 da->ops->createlocalvector = DMCreateLocalVector_DA; 332 da->ops->createinterpolation = DMCreateInterpolation_DA; 333 da->ops->getcoloring = DMCreateColoring_DA; 334 da->ops->creatematrix = DMCreateMatrix_DA; 335 da->ops->refine = DMRefine_DA; 336 da->ops->coarsen = DMCoarsen_DA; 337 da->ops->refinehierarchy = DMRefineHierarchy_DA; 338 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 339 da->ops->getinjection = DMCreateInjection_DA; 340 da->ops->getaggregates = DMCreateAggregates_DA; 341 da->ops->destroy = DMDestroy_DA; 342 da->ops->view = 0; 343 da->ops->setfromoptions = DMSetFromOptions_DA; 344 da->ops->setup = DMSetUp_DA; 345 da->ops->load = DMLoad_DA; 346 da->ops->createcoordinatedm = DMCreateCoordinateDM_DA; 347 da->ops->createsubdm = DMCreateSubDM_DA; 348 da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA; 349 da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA; 350 da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA; 351 PetscFunctionReturn(0); 352 } 353 EXTERN_C_END 354 355 #undef __FUNCT__ 356 #define __FUNCT__ "DMDACreate" 357 /*@ 358 DMDACreate - Creates a DMDA object. 359 360 Collective on MPI_Comm 361 362 Input Parameter: 363 . comm - The communicator for the DMDA object 364 365 Output Parameter: 366 . da - The DMDA object 367 368 Level: advanced 369 370 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 371 372 .keywords: DMDA, create 373 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 374 @*/ 375 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 376 { 377 PetscErrorCode ierr; 378 379 PetscFunctionBegin; 380 PetscValidPointer(da,2); 381 ierr = DMCreate(comm,da);CHKERRQ(ierr); 382 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 383 PetscFunctionReturn(0); 384 } 385 386 387