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