1 2 /* 3 Code for manipulating distributed regular 1d arrays in parallel. 4 This file was created by Peter Mell 6/30/95 5 */ 6 7 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 8 9 #include <petscdraw.h> 10 static PetscErrorCode DMView_DA_1d(DM da, PetscViewer viewer) { 11 PetscMPIInt rank; 12 PetscBool iascii, isdraw, isglvis, isbinary; 13 DM_DA *dd = (DM_DA *)da->data; 14 #if defined(PETSC_HAVE_MATLAB_ENGINE) 15 PetscBool ismatlab; 16 #endif 17 18 PetscFunctionBegin; 19 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank)); 20 21 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 22 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 23 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis)); 24 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 25 #if defined(PETSC_HAVE_MATLAB_ENGINE) 26 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab)); 27 #endif 28 if (iascii) { 29 PetscViewerFormat format; 30 31 PetscCall(PetscViewerGetFormat(viewer, &format)); 32 if (format == PETSC_VIEWER_LOAD_BALANCE) { 33 PetscInt i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal; 34 DMDALocalInfo info; 35 PetscMPIInt size; 36 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size)); 37 PetscCall(DMDAGetLocalInfo(da, &info)); 38 nzlocal = info.xm; 39 PetscCall(PetscMalloc1(size, &nz)); 40 PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da))); 41 for (i = 0; i < (PetscInt)size; i++) { 42 nmax = PetscMax(nmax, nz[i]); 43 nmin = PetscMin(nmin, nz[i]); 44 navg += nz[i]; 45 } 46 PetscCall(PetscFree(nz)); 47 navg = navg / size; 48 PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax)); 49 PetscFunctionReturn(0); 50 } 51 if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) { 52 DMDALocalInfo info; 53 PetscCall(DMDAGetLocalInfo(da, &info)); 54 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 55 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " m %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->m, dd->w, dd->s)); 56 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs, info.xs + info.xm)); 57 PetscCall(PetscViewerFlush(viewer)); 58 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 59 } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer)); 60 else PetscCall(DMView_DA_VTK(da, viewer)); 61 } else if (isdraw) { 62 PetscDraw draw; 63 double ymin = -1, ymax = 1, xmin = -1, xmax = dd->M, x; 64 PetscInt base; 65 char node[10]; 66 PetscBool isnull; 67 68 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 69 PetscCall(PetscDrawIsNull(draw, &isnull)); 70 if (isnull) PetscFunctionReturn(0); 71 72 PetscCall(PetscDrawCheckResizedWindow(draw)); 73 PetscCall(PetscDrawClear(draw)); 74 PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax)); 75 76 PetscDrawCollectiveBegin(draw); 77 /* first processor draws all node lines */ 78 if (rank == 0) { 79 PetscInt xmin_tmp; 80 ymin = 0.0; 81 ymax = 0.3; 82 for (xmin_tmp = 0; xmin_tmp < dd->M; xmin_tmp++) PetscCall(PetscDrawLine(draw, (double)xmin_tmp, ymin, (double)xmin_tmp, ymax, PETSC_DRAW_BLACK)); 83 xmin = 0.0; 84 xmax = dd->M - 1; 85 PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK)); 86 PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_BLACK)); 87 } 88 PetscDrawCollectiveEnd(draw); 89 PetscCall(PetscDrawFlush(draw)); 90 PetscCall(PetscDrawPause(draw)); 91 92 PetscDrawCollectiveBegin(draw); 93 /* draw my box */ 94 ymin = 0; 95 ymax = 0.3; 96 xmin = dd->xs / dd->w; 97 xmax = (dd->xe / dd->w) - 1; 98 PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED)); 99 PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED)); 100 PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED)); 101 PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED)); 102 /* Put in index numbers */ 103 base = dd->base / dd->w; 104 for (x = xmin; x <= xmax; x++) { 105 PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++)); 106 PetscCall(PetscDrawString(draw, x, ymin, PETSC_DRAW_RED, node)); 107 } 108 PetscDrawCollectiveEnd(draw); 109 PetscCall(PetscDrawFlush(draw)); 110 PetscCall(PetscDrawPause(draw)); 111 PetscCall(PetscDrawSave(draw)); 112 } else if (isglvis) { 113 PetscCall(DMView_DA_GLVis(da, viewer)); 114 } else if (isbinary) { 115 PetscCall(DMView_DA_Binary(da, viewer)); 116 #if defined(PETSC_HAVE_MATLAB_ENGINE) 117 } else if (ismatlab) { 118 PetscCall(DMView_DA_Matlab(da, viewer)); 119 #endif 120 } 121 PetscFunctionReturn(0); 122 } 123 124 PetscErrorCode DMSetUp_DA_1D(DM da) { 125 DM_DA *dd = (DM_DA *)da->data; 126 const PetscInt M = dd->M; 127 const PetscInt dof = dd->w; 128 const PetscInt s = dd->s; 129 const PetscInt sDist = s; /* stencil distance in points */ 130 const PetscInt *lx = dd->lx; 131 DMBoundaryType bx = dd->bx; 132 MPI_Comm comm; 133 Vec local, global; 134 VecScatter gtol; 135 IS to, from; 136 PetscBool flg1 = PETSC_FALSE, flg2 = PETSC_FALSE; 137 PetscMPIInt rank, size; 138 PetscInt i, *idx, nn, left, xs, xe, x, Xs, Xe, start, m, IXs, IXe; 139 140 PetscFunctionBegin; 141 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 142 PetscCallMPI(MPI_Comm_size(comm, &size)); 143 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 144 145 dd->p = 1; 146 dd->n = 1; 147 dd->m = size; 148 m = dd->m; 149 150 if (s > 0) { 151 /* if not communicating data then should be ok to have nothing on some processes */ 152 PetscCheck(M >= m, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "More processes than data points! %" PetscInt_FMT " %" PetscInt_FMT, m, M); 153 PetscCheck((M - 1) >= s || size <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Array is too small for stencil! %" PetscInt_FMT " %" PetscInt_FMT, M - 1, s); 154 } 155 156 /* 157 Determine locally owned region 158 xs is the first local node number, x is the number of local nodes 159 */ 160 if (!lx) { 161 PetscCall(PetscMalloc1(m, &dd->lx)); 162 PetscCall(PetscOptionsGetBool(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_partition_blockcomm", &flg1, NULL)); 163 PetscCall(PetscOptionsGetBool(((PetscObject)da)->options, ((PetscObject)da)->prefix, "-da_partition_nodes_at_end", &flg2, NULL)); 164 if (flg1) { /* Block Comm type Distribution */ 165 xs = rank * M / m; 166 x = (rank + 1) * M / m - xs; 167 } else if (flg2) { /* The odd nodes are evenly distributed across last nodes */ 168 x = (M + rank) / m; 169 if (M / m == x) xs = rank * x; 170 else xs = rank * (x - 1) + (M + rank) % (x * m); 171 } else { /* The odd nodes are evenly distributed across the first k nodes */ 172 /* Regular PETSc Distribution */ 173 x = M / m + ((M % m) > rank); 174 if (rank >= (M % m)) xs = (rank * (PetscInt)(M / m) + M % m); 175 else xs = rank * (PetscInt)(M / m) + rank; 176 } 177 PetscCallMPI(MPI_Allgather(&xs, 1, MPIU_INT, dd->lx, 1, MPIU_INT, comm)); 178 for (i = 0; i < m - 1; i++) dd->lx[i] = dd->lx[i + 1] - dd->lx[i]; 179 dd->lx[m - 1] = M - dd->lx[m - 1]; 180 } else { 181 x = lx[rank]; 182 xs = 0; 183 for (i = 0; i < rank; i++) xs += lx[i]; 184 /* verify that data user provided is consistent */ 185 left = xs; 186 for (i = rank; i < size; i++) left += lx[i]; 187 PetscCheck(left == M, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of lx across processors not equal to M %" PetscInt_FMT " %" PetscInt_FMT, left, M); 188 } 189 190 /* 191 check if the scatter requires more than one process neighbor or wraps around 192 the domain more than once 193 */ 194 PetscCheck((x >= s) || ((M <= 1) && (bx != DM_BOUNDARY_PERIODIC)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local x-width of domain x %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, x, s); 195 196 xe = xs + x; 197 198 /* determine ghost region (Xs) and region scattered into (IXs) */ 199 if (xs - sDist > 0) { 200 Xs = xs - sDist; 201 IXs = xs - sDist; 202 } else { 203 if (bx) Xs = xs - sDist; 204 else Xs = 0; 205 IXs = 0; 206 } 207 if (xe + sDist <= M) { 208 Xe = xe + sDist; 209 IXe = xe + sDist; 210 } else { 211 if (bx) Xe = xe + sDist; 212 else Xe = M; 213 IXe = M; 214 } 215 216 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { 217 Xs = xs - sDist; 218 Xe = xe + sDist; 219 IXs = xs - sDist; 220 IXe = xe + sDist; 221 } 222 223 /* allocate the base parallel and sequential vectors */ 224 dd->Nlocal = dof * x; 225 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global)); 226 dd->nlocal = dof * (Xe - Xs); 227 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local)); 228 229 PetscCall(VecGetOwnershipRange(global, &start, NULL)); 230 231 /* Create Global to Local Vector Scatter Context */ 232 /* global to local must retrieve ghost points */ 233 PetscCall(ISCreateStride(comm, dof * (IXe - IXs), dof * (IXs - Xs), 1, &to)); 234 235 PetscCall(PetscMalloc1(x + 2 * sDist, &idx)); 236 237 for (i = 0; i < IXs - Xs; i++) idx[i] = -1; /* prepend with -1s if needed for ghosted case*/ 238 239 nn = IXs - Xs; 240 if (bx == DM_BOUNDARY_PERIODIC) { /* Handle all cases with periodic first */ 241 for (i = 0; i < sDist; i++) { /* Left ghost points */ 242 if ((xs - sDist + i) >= 0) idx[nn++] = xs - sDist + i; 243 else idx[nn++] = M + (xs - sDist + i); 244 } 245 246 for (i = 0; i < x; i++) idx[nn++] = xs + i; /* Non-ghost points */ 247 248 for (i = 0; i < sDist; i++) { /* Right ghost points */ 249 if ((xe + i) < M) idx[nn++] = xe + i; 250 else idx[nn++] = (xe + i) - M; 251 } 252 } else if (bx == DM_BOUNDARY_MIRROR) { /* Handle all cases with periodic first */ 253 for (i = 0; i < (sDist); i++) { /* Left ghost points */ 254 if ((xs - sDist + i) >= 0) idx[nn++] = xs - sDist + i; 255 else idx[nn++] = sDist - i; 256 } 257 258 for (i = 0; i < x; i++) idx[nn++] = xs + i; /* Non-ghost points */ 259 260 for (i = 0; i < (sDist); i++) { /* Right ghost points */ 261 if ((xe + i) < M) idx[nn++] = xe + i; 262 else idx[nn++] = M - (i + 2); 263 } 264 } else { /* Now do all cases with no periodicity */ 265 if (0 <= xs - sDist) { 266 for (i = 0; i < sDist; i++) idx[nn++] = xs - sDist + i; 267 } else { 268 for (i = 0; i < xs; i++) idx[nn++] = i; 269 } 270 271 for (i = 0; i < x; i++) idx[nn++] = xs + i; 272 273 if ((xe + sDist) <= M) { 274 for (i = 0; i < sDist; i++) idx[nn++] = xe + i; 275 } else { 276 for (i = xe; i < M; i++) idx[nn++] = i; 277 } 278 } 279 280 PetscCall(ISCreateBlock(comm, dof, nn - IXs + Xs, &idx[IXs - Xs], PETSC_USE_POINTER, &from)); 281 PetscCall(VecScatterCreate(global, from, local, to, >ol)); 282 PetscCall(ISDestroy(&to)); 283 PetscCall(ISDestroy(&from)); 284 PetscCall(VecDestroy(&local)); 285 PetscCall(VecDestroy(&global)); 286 287 dd->xs = dof * xs; 288 dd->xe = dof * xe; 289 dd->ys = 0; 290 dd->ye = 1; 291 dd->zs = 0; 292 dd->ze = 1; 293 dd->Xs = dof * Xs; 294 dd->Xe = dof * Xe; 295 dd->Ys = 0; 296 dd->Ye = 1; 297 dd->Zs = 0; 298 dd->Ze = 1; 299 300 dd->gtol = gtol; 301 dd->base = dof * xs; 302 da->ops->view = DMView_DA_1d; 303 304 /* 305 Set the local to global ordering in the global vector, this allows use 306 of VecSetValuesLocal(). 307 */ 308 for (i = 0; i < Xe - IXe; i++) idx[nn++] = -1; /* pad with -1s if needed for ghosted case*/ 309 310 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap)); 311 312 PetscFunctionReturn(0); 313 } 314 315 /*@C 316 DMDACreate1d - Creates an object that will manage the communication of one-dimensional 317 regular array data that is distributed across some processors. 318 319 Collective 320 321 Input Parameters: 322 + comm - MPI communicator 323 . bx - type of ghost cells at the boundary the array should have, if any. Use 324 DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, or DM_BOUNDARY_PERIODIC. 325 . M - global dimension of the array (that is the number of grid points) 326 from the command line with -da_grid_x <M>) 327 . dof - number of degrees of freedom per node 328 . s - stencil width 329 - lx - array containing number of nodes in the X direction on each processor, 330 or NULL. If non-null, must be of length as the number of processes in the MPI_Comm. 331 The sum of these entries must equal M 332 333 Output Parameter: 334 . da - the resulting distributed array object 335 336 Options Database Key: 337 + -dm_view - Calls DMView() at the conclusion of DMDACreate1d() 338 . -da_grid_x <nx> - number of grid points in x direction 339 . -da_refine_x <rx> - refinement factor 340 - -da_refine <n> - refine the DMDA n times before creating it 341 342 Level: beginner 343 344 Notes: 345 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 346 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 347 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 348 349 You must call DMSetUp() after this call before using this DM. 350 351 If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call 352 but before DMSetUp(). 353 354 .seealso: `DMDestroy()`, `DMView()`, `DMDACreate2d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDASetRefinementFactor()`, 355 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDAGetRefinementFactor()`, 356 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 357 `DMStagCreate1d()` 358 359 @*/ 360 PetscErrorCode DMDACreate1d(MPI_Comm comm, DMBoundaryType bx, PetscInt M, PetscInt dof, PetscInt s, const PetscInt lx[], DM *da) { 361 PetscMPIInt size; 362 363 PetscFunctionBegin; 364 PetscCall(DMDACreate(comm, da)); 365 PetscCall(DMSetDimension(*da, 1)); 366 PetscCall(DMDASetSizes(*da, M, 1, 1)); 367 PetscCallMPI(MPI_Comm_size(comm, &size)); 368 PetscCall(DMDASetNumProcs(*da, size, PETSC_DECIDE, PETSC_DECIDE)); 369 PetscCall(DMDASetBoundaryType(*da, bx, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE)); 370 PetscCall(DMDASetDof(*da, dof)); 371 PetscCall(DMDASetStencilWidth(*da, s)); 372 PetscCall(DMDASetOwnershipRanges(*da, lx, NULL, NULL)); 373 PetscFunctionReturn(0); 374 } 375