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__ "DMCreateFieldDecomposition_DA" 173 PetscErrorCode DMCreateFieldDecomposition_DA(DM dm, PetscInt *len,char ***namelist, IS **islist, DM** dmlist) 174 { 175 PetscInt i; 176 PetscErrorCode ierr; 177 DM_DA *dd = (DM_DA*)dm->data; 178 PetscInt dof = dd->w; 179 180 PetscFunctionBegin; 181 if (len) *len = dof; 182 if (islist) { 183 Vec v; 184 PetscInt rstart,n; 185 186 ierr = DMGetGlobalVector(dm,&v);CHKERRQ(ierr); 187 ierr = VecGetOwnershipRange(v,&rstart,PETSC_NULL);CHKERRQ(ierr); 188 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 189 ierr = DMRestoreGlobalVector(dm,&v);CHKERRQ(ierr); 190 ierr = PetscMalloc(dof*sizeof(IS),islist);CHKERRQ(ierr); 191 for (i=0; i<dof; i++) { 192 ierr = ISCreateStride(((PetscObject)dm)->comm,n/dof,rstart+i,dof,&(*islist)[i]);CHKERRQ(ierr); 193 } 194 } 195 if (namelist) { 196 ierr = PetscMalloc(dof*sizeof(const char *), namelist);CHKERRQ(ierr); 197 if (dd->fieldname) { 198 for (i=0; i<dof; i++) { 199 ierr = PetscStrallocpy(dd->fieldname[i],&(*namelist)[i]);CHKERRQ(ierr); 200 } 201 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently DMDA must have fieldnames"); 202 } 203 if (dmlist) { 204 DM da; 205 206 ierr = DMDACreate(((PetscObject)dm)->comm, &da);CHKERRQ(ierr); 207 ierr = DMDASetDim(da, dd->dim);CHKERRQ(ierr); 208 ierr = DMDASetSizes(da, dd->M, dd->N, dd->P);CHKERRQ(ierr); 209 ierr = DMDASetNumProcs(da, dd->m, dd->n, dd->p);CHKERRQ(ierr); 210 ierr = DMDASetBoundaryType(da, dd->bx, dd->by, dd->bz);CHKERRQ(ierr); 211 ierr = DMDASetDof(da, 1);CHKERRQ(ierr); 212 ierr = DMDASetStencilType(da, dd->stencil_type);CHKERRQ(ierr); 213 ierr = DMDASetStencilWidth(da, dd->s);CHKERRQ(ierr); 214 ierr = DMSetUp(da);CHKERRQ(ierr); 215 ierr = PetscMalloc(dof*sizeof(DM),dmlist);CHKERRQ(ierr); 216 for (i=0; i<dof-1; i++) {ierr = PetscObjectReference((PetscObject)da);CHKERRQ(ierr);} 217 for (i=0; i<dof; i++) (*dmlist)[i] = da; 218 } 219 220 PetscFunctionReturn(0); 221 } 222 223 /*MC 224 DMDA = "da" - A DM object that is used to manage data for a structured grid in 1, 2, or 3 dimensions. 225 In the global representation of the vector each process stores a non-overlapping rectangular (or slab in 3d) portion of the grid points. 226 In the local representation these rectangular regions (slabs) are extended in all directions by a stencil width. 227 228 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 229 vertex centered. 230 231 232 Level: intermediate 233 234 .seealso: DMType, DMCOMPOSITE, DMDACreate(), DMCreate(), DMSetType() 235 M*/ 236 237 238 EXTERN_C_BEGIN 239 #undef __FUNCT__ 240 #define __FUNCT__ "DMCreate_DA" 241 PetscErrorCode DMCreate_DA(DM da) 242 { 243 PetscErrorCode ierr; 244 DM_DA *dd; 245 246 PetscFunctionBegin; 247 PetscValidPointer(da,1); 248 ierr = PetscNewLog(da,DM_DA,&dd);CHKERRQ(ierr); 249 da->data = dd; 250 251 dd->dim = -1; 252 dd->interptype = DMDA_Q1; 253 dd->refine_x = 2; 254 dd->refine_y = 2; 255 dd->refine_z = 2; 256 dd->coarsen_x = 2; 257 dd->coarsen_y = 2; 258 dd->coarsen_z = 2; 259 dd->fieldname = PETSC_NULL; 260 dd->nlocal = -1; 261 dd->Nlocal = -1; 262 dd->M = -1; 263 dd->N = -1; 264 dd->P = -1; 265 dd->m = -1; 266 dd->n = -1; 267 dd->p = -1; 268 dd->w = -1; 269 dd->s = -1; 270 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 271 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 272 273 dd->overlap = 0; 274 dd->xo = 0; 275 dd->yo = 0; 276 dd->zo = 0; 277 dd->Mo = -1; 278 dd->No = -1; 279 dd->Po = -1; 280 281 dd->gtol = PETSC_NULL; 282 dd->ltog = PETSC_NULL; 283 dd->ltol = PETSC_NULL; 284 dd->ao = PETSC_NULL; 285 dd->base = -1; 286 dd->bx = DMDA_BOUNDARY_NONE; 287 dd->by = DMDA_BOUNDARY_NONE; 288 dd->bz = DMDA_BOUNDARY_NONE; 289 dd->stencil_type = DMDA_STENCIL_BOX; 290 dd->interptype = DMDA_Q1; 291 dd->idx = PETSC_NULL; 292 dd->Nl = -1; 293 dd->lx = PETSC_NULL; 294 dd->ly = PETSC_NULL; 295 dd->lz = PETSC_NULL; 296 297 dd->elementtype = DMDA_ELEMENT_Q1; 298 299 ierr = PetscStrallocpy(VECSTANDARD,(char**)&da->vectype);CHKERRQ(ierr); 300 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 301 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 302 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 303 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 304 da->ops->createglobalvector = DMCreateGlobalVector_DA; 305 da->ops->createlocalvector = DMCreateLocalVector_DA; 306 da->ops->createinterpolation = DMCreateInterpolation_DA; 307 da->ops->getcoloring = DMCreateColoring_DA; 308 da->ops->creatematrix = DMCreateMatrix_DA; 309 da->ops->refine = DMRefine_DA; 310 da->ops->coarsen = DMCoarsen_DA; 311 da->ops->refinehierarchy = DMRefineHierarchy_DA; 312 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 313 da->ops->getinjection = DMCreateInjection_DA; 314 da->ops->getaggregates = DMCreateAggregates_DA; 315 da->ops->destroy = DMDestroy_DA; 316 da->ops->view = 0; 317 da->ops->setfromoptions = DMSetFromOptions_DA; 318 da->ops->setup = DMSetUp_DA; 319 da->ops->load = DMLoad_DA; 320 da->ops->createcoordinatedm = DMCreateCoordinateDM_DA; 321 da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA; 322 da->ops->createdomaindecomposition = DMCreateDomainDecomposition_DA; 323 da->ops->createddscatters = DMCreateDomainDecompositionScatters_DA; 324 PetscFunctionReturn(0); 325 } 326 EXTERN_C_END 327 328 #undef __FUNCT__ 329 #define __FUNCT__ "DMDACreate" 330 /*@ 331 DMDACreate - Creates a DMDA object. 332 333 Collective on MPI_Comm 334 335 Input Parameter: 336 . comm - The communicator for the DMDA object 337 338 Output Parameter: 339 . da - The DMDA object 340 341 Level: advanced 342 343 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 344 345 .keywords: DMDA, create 346 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 347 @*/ 348 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 349 { 350 PetscErrorCode ierr; 351 352 PetscFunctionBegin; 353 PetscValidPointer(da,2); 354 ierr = DMCreate(comm,da);CHKERRQ(ierr); 355 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 356 PetscFunctionReturn(0); 357 } 358 359 360