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 226 EXTERN_C_BEGIN 227 #undef __FUNCT__ 228 #define __FUNCT__ "DMCreate_DA" 229 PetscErrorCode DMCreate_DA(DM da) 230 { 231 PetscErrorCode ierr; 232 DM_DA *dd; 233 234 PetscFunctionBegin; 235 PetscValidPointer(da,1); 236 ierr = PetscNewLog(da,DM_DA,&dd);CHKERRQ(ierr); 237 da->data = dd; 238 239 dd->dim = -1; 240 dd->interptype = DMDA_Q1; 241 dd->refine_x = 2; 242 dd->refine_y = 2; 243 dd->refine_z = 2; 244 dd->coarsen_x = 2; 245 dd->coarsen_y = 2; 246 dd->coarsen_z = 2; 247 dd->fieldname = PETSC_NULL; 248 dd->nlocal = -1; 249 dd->Nlocal = -1; 250 dd->M = -1; 251 dd->N = -1; 252 dd->P = -1; 253 dd->m = -1; 254 dd->n = -1; 255 dd->p = -1; 256 dd->w = -1; 257 dd->s = -1; 258 dd->xs = -1; dd->xe = -1; dd->ys = -1; dd->ye = -1; dd->zs = -1; dd->ze = -1; 259 dd->Xs = -1; dd->Xe = -1; dd->Ys = -1; dd->Ye = -1; dd->Zs = -1; dd->Ze = -1; 260 261 dd->overlap = 0; 262 263 dd->gtol = PETSC_NULL; 264 dd->ltog = PETSC_NULL; 265 dd->ltol = PETSC_NULL; 266 dd->ao = PETSC_NULL; 267 dd->base = -1; 268 dd->bx = DMDA_BOUNDARY_NONE; 269 dd->by = DMDA_BOUNDARY_NONE; 270 dd->bz = DMDA_BOUNDARY_NONE; 271 dd->stencil_type = DMDA_STENCIL_BOX; 272 dd->interptype = DMDA_Q1; 273 dd->idx = PETSC_NULL; 274 dd->Nl = -1; 275 dd->lx = PETSC_NULL; 276 dd->ly = PETSC_NULL; 277 dd->lz = PETSC_NULL; 278 279 dd->elementtype = DMDA_ELEMENT_Q1; 280 281 ierr = PetscStrallocpy(VECSTANDARD,&da->vectype);CHKERRQ(ierr); 282 da->ops->globaltolocalbegin = DMGlobalToLocalBegin_DA; 283 da->ops->globaltolocalend = DMGlobalToLocalEnd_DA; 284 da->ops->localtoglobalbegin = DMLocalToGlobalBegin_DA; 285 da->ops->localtoglobalend = DMLocalToGlobalEnd_DA; 286 da->ops->createglobalvector = DMCreateGlobalVector_DA; 287 da->ops->createlocalvector = DMCreateLocalVector_DA; 288 da->ops->createinterpolation = DMCreateInterpolation_DA; 289 da->ops->getcoloring = DMCreateColoring_DA; 290 da->ops->creatematrix = DMCreateMatrix_DA; 291 da->ops->refine = DMRefine_DA; 292 da->ops->coarsen = DMCoarsen_DA; 293 da->ops->refinehierarchy = DMRefineHierarchy_DA; 294 da->ops->coarsenhierarchy = DMCoarsenHierarchy_DA; 295 da->ops->getinjection = DMCreateInjection_DA; 296 da->ops->getaggregates = DMCreateAggregates_DA; 297 da->ops->destroy = DMDestroy_DA; 298 da->ops->view = 0; 299 da->ops->setfromoptions = DMSetFromOptions_DA; 300 da->ops->setup = DMSetUp_DA; 301 da->ops->load = DMLoad_DA; 302 da->ops->createfielddecomposition = DMCreateFieldDecomposition_DA; 303 PetscFunctionReturn(0); 304 } 305 EXTERN_C_END 306 307 #undef __FUNCT__ 308 #define __FUNCT__ "DMDACreate" 309 /*@ 310 DMDACreate - Creates a DMDA object. 311 312 Collective on MPI_Comm 313 314 Input Parameter: 315 . comm - The communicator for the DMDA object 316 317 Output Parameter: 318 . da - The DMDA object 319 320 Level: advanced 321 322 Developers Note: Since there exists DMDACreate1/2/3d() should this routine even exist? 323 324 .keywords: DMDA, create 325 .seealso: DMDASetSizes(), DMDADuplicate(), DMDACreate1d(), DMDACreate2d(), DMDACreate3d() 326 @*/ 327 PetscErrorCode DMDACreate(MPI_Comm comm, DM *da) 328 { 329 PetscErrorCode ierr; 330 331 PetscFunctionBegin; 332 PetscValidPointer(da,2); 333 ierr = DMCreate(comm,da);CHKERRQ(ierr); 334 ierr = DMSetType(*da,DMDA);CHKERRQ(ierr); 335 PetscFunctionReturn(0); 336 } 337 338 339