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