1 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 #include <petscdraw.h> 4 5 static PetscErrorCode DMView_DA_2d(DM da, PetscViewer viewer) { 6 PetscMPIInt rank; 7 PetscBool iascii, isdraw, isglvis, isbinary; 8 DM_DA *dd = (DM_DA *)da->data; 9 #if defined(PETSC_HAVE_MATLAB) 10 PetscBool ismatlab; 11 #endif 12 13 PetscFunctionBegin; 14 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank)); 15 16 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 17 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw)); 18 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERGLVIS, &isglvis)); 19 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &isbinary)); 20 #if defined(PETSC_HAVE_MATLAB) 21 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERMATLAB, &ismatlab)); 22 #endif 23 if (iascii) { 24 PetscViewerFormat format; 25 26 PetscCall(PetscViewerGetFormat(viewer, &format)); 27 if (format == PETSC_VIEWER_LOAD_BALANCE) { 28 PetscInt i, nmax = 0, nmin = PETSC_MAX_INT, navg = 0, *nz, nzlocal; 29 DMDALocalInfo info; 30 PetscMPIInt size; 31 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size)); 32 PetscCall(DMDAGetLocalInfo(da, &info)); 33 nzlocal = info.xm * info.ym; 34 PetscCall(PetscMalloc1(size, &nz)); 35 PetscCallMPI(MPI_Allgather(&nzlocal, 1, MPIU_INT, nz, 1, MPIU_INT, PetscObjectComm((PetscObject)da))); 36 for (i = 0; i < (PetscInt)size; i++) { 37 nmax = PetscMax(nmax, nz[i]); 38 nmin = PetscMin(nmin, nz[i]); 39 navg += nz[i]; 40 } 41 PetscCall(PetscFree(nz)); 42 navg = navg / size; 43 PetscCall(PetscViewerASCIIPrintf(viewer, " Load Balance - Grid Points: Min %" PetscInt_FMT " avg %" PetscInt_FMT " max %" PetscInt_FMT "\n", nmin, navg, nmax)); 44 PetscFunctionReturn(0); 45 } 46 if (format != PETSC_VIEWER_ASCII_VTK_DEPRECATED && format != PETSC_VIEWER_ASCII_VTK_CELL_DEPRECATED && format != PETSC_VIEWER_ASCII_GLVIS) { 47 DMDALocalInfo info; 48 PetscCall(DMDAGetLocalInfo(da, &info)); 49 PetscCall(PetscViewerASCIIPushSynchronized(viewer)); 50 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "Processor [%d] M %" PetscInt_FMT " N %" PetscInt_FMT " m %" PetscInt_FMT " n %" PetscInt_FMT " w %" PetscInt_FMT " s %" PetscInt_FMT "\n", rank, dd->M, dd->N, dd->m, dd->n, dd->w, dd->s)); 51 PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "X range of indices: %" PetscInt_FMT " %" PetscInt_FMT ", Y range of indices: %" PetscInt_FMT " %" PetscInt_FMT "\n", info.xs, info.xs + info.xm, info.ys, info.ys + info.ym)); 52 PetscCall(PetscViewerFlush(viewer)); 53 PetscCall(PetscViewerASCIIPopSynchronized(viewer)); 54 } else if (format == PETSC_VIEWER_ASCII_GLVIS) PetscCall(DMView_DA_GLVis(da, viewer)); 55 else PetscCall(DMView_DA_VTK(da, viewer)); 56 } else if (isdraw) { 57 PetscDraw draw; 58 double ymin = -1 * dd->s - 1, ymax = dd->N + dd->s; 59 double xmin = -1 * dd->s - 1, xmax = dd->M + dd->s; 60 double x, y; 61 PetscInt base; 62 const PetscInt *idx; 63 char node[10]; 64 PetscBool isnull; 65 66 PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw)); 67 PetscCall(PetscDrawIsNull(draw, &isnull)); 68 if (isnull) PetscFunctionReturn(0); 69 70 PetscCall(PetscDrawCheckResizedWindow(draw)); 71 PetscCall(PetscDrawClear(draw)); 72 PetscCall(PetscDrawSetCoordinates(draw, xmin, ymin, xmax, ymax)); 73 74 PetscDrawCollectiveBegin(draw); 75 /* first processor draw all node lines */ 76 if (rank == 0) { 77 ymin = 0.0; 78 ymax = dd->N - 1; 79 for (xmin = 0; xmin < dd->M; xmin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_BLACK)); 80 xmin = 0.0; 81 xmax = dd->M - 1; 82 for (ymin = 0; ymin < dd->N; ymin++) PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_BLACK)); 83 } 84 PetscDrawCollectiveEnd(draw); 85 PetscCall(PetscDrawFlush(draw)); 86 PetscCall(PetscDrawPause(draw)); 87 88 PetscDrawCollectiveBegin(draw); 89 /* draw my box */ 90 xmin = dd->xs / dd->w; 91 xmax = (dd->xe - 1) / dd->w; 92 ymin = dd->ys; 93 ymax = dd->ye - 1; 94 PetscCall(PetscDrawLine(draw, xmin, ymin, xmax, ymin, PETSC_DRAW_RED)); 95 PetscCall(PetscDrawLine(draw, xmin, ymin, xmin, ymax, PETSC_DRAW_RED)); 96 PetscCall(PetscDrawLine(draw, xmin, ymax, xmax, ymax, PETSC_DRAW_RED)); 97 PetscCall(PetscDrawLine(draw, xmax, ymin, xmax, ymax, PETSC_DRAW_RED)); 98 /* put in numbers */ 99 base = (dd->base) / dd->w; 100 for (y = ymin; y <= ymax; y++) { 101 for (x = xmin; x <= xmax; x++) { 102 PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)base++)); 103 PetscCall(PetscDrawString(draw, x, y, PETSC_DRAW_BLACK, node)); 104 } 105 } 106 PetscDrawCollectiveEnd(draw); 107 PetscCall(PetscDrawFlush(draw)); 108 PetscCall(PetscDrawPause(draw)); 109 110 PetscDrawCollectiveBegin(draw); 111 /* overlay ghost numbers, useful for error checking */ 112 PetscCall(ISLocalToGlobalMappingGetBlockIndices(da->ltogmap, &idx)); 113 base = 0; 114 xmin = dd->Xs; 115 xmax = dd->Xe; 116 ymin = dd->Ys; 117 ymax = dd->Ye; 118 for (y = ymin; y < ymax; y++) { 119 for (x = xmin; x < xmax; x++) { 120 if ((base % dd->w) == 0) { 121 PetscCall(PetscSNPrintf(node, sizeof(node), "%d", (int)(idx[base / dd->w]))); 122 PetscCall(PetscDrawString(draw, x / dd->w, y, PETSC_DRAW_BLUE, node)); 123 } 124 base++; 125 } 126 } 127 PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(da->ltogmap, &idx)); 128 PetscDrawCollectiveEnd(draw); 129 PetscCall(PetscDrawFlush(draw)); 130 PetscCall(PetscDrawPause(draw)); 131 PetscCall(PetscDrawSave(draw)); 132 } else if (isglvis) { 133 PetscCall(DMView_DA_GLVis(da, viewer)); 134 } else if (isbinary) { 135 PetscCall(DMView_DA_Binary(da, viewer)); 136 #if defined(PETSC_HAVE_MATLAB) 137 } else if (ismatlab) { 138 PetscCall(DMView_DA_Matlab(da, viewer)); 139 #endif 140 } 141 PetscFunctionReturn(0); 142 } 143 144 #if defined(new) 145 /* 146 DMDAGetDiagonal_MFFD - Gets the diagonal for a matrix free matrix where local 147 function lives on a DMDA 148 149 y ~= (F(u + ha) - F(u))/h, 150 where F = nonlinear function, as set by SNESSetFunction() 151 u = current iterate 152 h = difference interval 153 */ 154 PetscErrorCode DMDAGetDiagonal_MFFD(DM da, Vec U, Vec a) { 155 PetscScalar h, *aa, *ww, v; 156 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON; 157 PetscInt gI, nI; 158 MatStencil stencil; 159 DMDALocalInfo info; 160 161 PetscFunctionBegin; 162 PetscCall((*ctx->func)(0, U, a, ctx->funcctx)); 163 PetscCall((*ctx->funcisetbase)(U, ctx->funcctx)); 164 165 PetscCall(VecGetArray(U, &ww)); 166 PetscCall(VecGetArray(a, &aa)); 167 168 nI = 0; 169 h = ww[gI]; 170 if (h == 0.0) h = 1.0; 171 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin; 172 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin; 173 h *= epsilon; 174 175 ww[gI] += h; 176 PetscCall((*ctx->funci)(i, w, &v, ctx->funcctx)); 177 aa[nI] = (v - aa[nI]) / h; 178 ww[gI] -= h; 179 nI++; 180 181 PetscCall(VecRestoreArray(U, &ww)); 182 PetscCall(VecRestoreArray(a, &aa)); 183 PetscFunctionReturn(0); 184 } 185 #endif 186 187 PetscErrorCode DMSetUp_DA_2D(DM da) { 188 DM_DA *dd = (DM_DA *)da->data; 189 const PetscInt M = dd->M; 190 const PetscInt N = dd->N; 191 PetscInt m = dd->m; 192 PetscInt n = dd->n; 193 const PetscInt dof = dd->w; 194 const PetscInt s = dd->s; 195 DMBoundaryType bx = dd->bx; 196 DMBoundaryType by = dd->by; 197 DMDAStencilType stencil_type = dd->stencil_type; 198 PetscInt *lx = dd->lx; 199 PetscInt *ly = dd->ly; 200 MPI_Comm comm; 201 PetscMPIInt rank, size; 202 PetscInt xs, xe, ys, ye, x, y, Xs, Xe, Ys, Ye, IXs, IXe, IYs, IYe; 203 PetscInt up, down, left, right, i, n0, n1, n2, n3, n5, n6, n7, n8, *idx, nn; 204 PetscInt xbase, *bases, *ldims, j, x_t, y_t, s_t, base, count; 205 PetscInt s_x, s_y; /* s proportionalized to w */ 206 PetscInt sn0 = 0, sn2 = 0, sn6 = 0, sn8 = 0; 207 Vec local, global; 208 VecScatter gtol; 209 IS to, from; 210 211 PetscFunctionBegin; 212 PetscCheck(stencil_type != DMDA_STENCIL_BOX || (bx != DM_BOUNDARY_MIRROR && by != DM_BOUNDARY_MIRROR), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Mirror boundary and box stencil"); 213 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 214 #if !defined(PETSC_USE_64BIT_INDICES) 215 PetscCheck(((PetscInt64)M) * ((PetscInt64)N) * ((PetscInt64)dof) <= (PetscInt64)PETSC_MPI_INT_MAX, comm, PETSC_ERR_INT_OVERFLOW, "Mesh of %" PetscInt_FMT " by %" PetscInt_FMT " by %" PetscInt_FMT " (dof) is too large for 32 bit indices", M, N, dof); 216 #endif 217 218 PetscCallMPI(MPI_Comm_size(comm, &size)); 219 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 220 221 dd->p = 1; 222 if (m != PETSC_DECIDE) { 223 PetscCheck(m >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m); 224 PetscCheck(m <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %d", m, size); 225 } 226 if (n != PETSC_DECIDE) { 227 PetscCheck(n >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %" PetscInt_FMT, n); 228 PetscCheck(n <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %d", n, size); 229 } 230 231 if (m == PETSC_DECIDE || n == PETSC_DECIDE) { 232 if (n != PETSC_DECIDE) { 233 m = size / n; 234 } else if (m != PETSC_DECIDE) { 235 n = size / m; 236 } else { 237 /* try for squarish distribution */ 238 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N))); 239 if (!m) m = 1; 240 while (m > 0) { 241 n = size / m; 242 if (m * n == size) break; 243 m--; 244 } 245 if (M > N && m < n) { 246 PetscInt _m = m; 247 m = n; 248 n = _m; 249 } 250 } 251 PetscCheck(m * n == size, comm, PETSC_ERR_PLIB, "Unable to create partition, check the size of the communicator and input m and n "); 252 } else PetscCheck(m * n == size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition"); 253 254 PetscCheck(M >= m, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m); 255 PetscCheck(N >= n, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n); 256 257 /* 258 Determine locally owned region 259 xs is the first local node number, x is the number of local nodes 260 */ 261 if (!lx) { 262 PetscCall(PetscMalloc1(m, &dd->lx)); 263 lx = dd->lx; 264 for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i); 265 } 266 x = lx[rank % m]; 267 xs = 0; 268 for (i = 0; i < (rank % m); i++) xs += lx[i]; 269 if (PetscDefined(USE_DEBUG)) { 270 left = xs; 271 for (i = (rank % m); i < m; i++) left += lx[i]; 272 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); 273 } 274 275 /* 276 Determine locally owned region 277 ys is the first local node number, y is the number of local nodes 278 */ 279 if (!ly) { 280 PetscCall(PetscMalloc1(n, &dd->ly)); 281 ly = dd->ly; 282 for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i); 283 } 284 y = ly[rank / m]; 285 ys = 0; 286 for (i = 0; i < (rank / m); i++) ys += ly[i]; 287 if (PetscDefined(USE_DEBUG)) { 288 left = ys; 289 for (i = (rank / m); i < n; i++) left += ly[i]; 290 PetscCheck(left == N, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Sum of ly across processors not equal to N: %" PetscInt_FMT " %" PetscInt_FMT, left, N); 291 } 292 293 /* 294 check if the scatter requires more than one process neighbor or wraps around 295 the domain more than once 296 */ 297 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); 298 PetscCheck((y >= s) || ((n <= 1) && (by != DM_BOUNDARY_PERIODIC)), PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Local y-width of domain y %" PetscInt_FMT " is smaller than stencil width s %" PetscInt_FMT, y, s); 299 xe = xs + x; 300 ye = ys + y; 301 302 /* determine ghost region (Xs) and region scattered into (IXs) */ 303 if (xs - s > 0) { 304 Xs = xs - s; 305 IXs = xs - s; 306 } else { 307 if (bx) { 308 Xs = xs - s; 309 } else { 310 Xs = 0; 311 } 312 IXs = 0; 313 } 314 if (xe + s <= M) { 315 Xe = xe + s; 316 IXe = xe + s; 317 } else { 318 if (bx) { 319 Xs = xs - s; 320 Xe = xe + s; 321 } else { 322 Xe = M; 323 } 324 IXe = M; 325 } 326 327 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { 328 IXs = xs - s; 329 IXe = xe + s; 330 Xs = xs - s; 331 Xe = xe + s; 332 } 333 334 if (ys - s > 0) { 335 Ys = ys - s; 336 IYs = ys - s; 337 } else { 338 if (by) { 339 Ys = ys - s; 340 } else { 341 Ys = 0; 342 } 343 IYs = 0; 344 } 345 if (ye + s <= N) { 346 Ye = ye + s; 347 IYe = ye + s; 348 } else { 349 if (by) { 350 Ye = ye + s; 351 } else { 352 Ye = N; 353 } 354 IYe = N; 355 } 356 357 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) { 358 IYs = ys - s; 359 IYe = ye + s; 360 Ys = ys - s; 361 Ye = ye + s; 362 } 363 364 /* stencil length in each direction */ 365 s_x = s; 366 s_y = s; 367 368 /* determine starting point of each processor */ 369 nn = x * y; 370 PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims)); 371 PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm)); 372 bases[0] = 0; 373 for (i = 1; i <= size; i++) bases[i] = ldims[i - 1]; 374 for (i = 1; i <= size; i++) bases[i] += bases[i - 1]; 375 base = bases[rank] * dof; 376 377 /* allocate the base parallel and sequential vectors */ 378 dd->Nlocal = x * y * dof; 379 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global)); 380 dd->nlocal = (Xe - Xs) * (Ye - Ys) * dof; 381 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local)); 382 383 /* generate global to local vector scatter and local to global mapping*/ 384 385 /* global to local must include ghost points within the domain, 386 but not ghost points outside the domain that aren't periodic */ 387 PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs), &idx)); 388 if (stencil_type == DMDA_STENCIL_BOX) { 389 left = IXs - Xs; 390 right = left + (IXe - IXs); 391 down = IYs - Ys; 392 up = down + (IYe - IYs); 393 count = 0; 394 for (i = down; i < up; i++) { 395 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 396 } 397 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 398 399 } else { 400 /* must drop into cross shape region */ 401 /* ---------| 402 | top | 403 |--- ---| up 404 | middle | 405 | | 406 ---- ---- down 407 | bottom | 408 ----------- 409 Xs xs xe Xe */ 410 left = xs - Xs; 411 right = left + x; 412 down = ys - Ys; 413 up = down + y; 414 count = 0; 415 /* bottom */ 416 for (i = (IYs - Ys); i < down; i++) { 417 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 418 } 419 /* middle */ 420 for (i = down; i < up; i++) { 421 for (j = (IXs - Xs); j < (IXe - Xs); j++) idx[count++] = j + i * (Xe - Xs); 422 } 423 /* top */ 424 for (i = up; i < up + IYe - ye; i++) { 425 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 426 } 427 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 428 } 429 430 /* determine who lies on each side of us stored in n6 n7 n8 431 n3 n5 432 n0 n1 n2 433 */ 434 435 /* Assume the Non-Periodic Case */ 436 n1 = rank - m; 437 if (rank % m) { 438 n0 = n1 - 1; 439 } else { 440 n0 = -1; 441 } 442 if ((rank + 1) % m) { 443 n2 = n1 + 1; 444 n5 = rank + 1; 445 n8 = rank + m + 1; 446 if (n8 >= m * n) n8 = -1; 447 } else { 448 n2 = -1; 449 n5 = -1; 450 n8 = -1; 451 } 452 if (rank % m) { 453 n3 = rank - 1; 454 n6 = n3 + m; 455 if (n6 >= m * n) n6 = -1; 456 } else { 457 n3 = -1; 458 n6 = -1; 459 } 460 n7 = rank + m; 461 if (n7 >= m * n) n7 = -1; 462 463 if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) { 464 /* Modify for Periodic Cases */ 465 /* Handle all four corners */ 466 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m - 1; 467 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0; 468 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size - m; 469 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size - 1; 470 471 /* Handle Top and Bottom Sides */ 472 if (n1 < 0) n1 = rank + m * (n - 1); 473 if (n7 < 0) n7 = rank - m * (n - 1); 474 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 475 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1; 476 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 477 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1; 478 479 /* Handle Left and Right Sides */ 480 if (n3 < 0) n3 = rank + (m - 1); 481 if (n5 < 0) n5 = rank - (m - 1); 482 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1; 483 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1; 484 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1; 485 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1; 486 } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */ 487 if (n1 < 0) n1 = rank + m * (n - 1); 488 if (n7 < 0) n7 = rank - m * (n - 1); 489 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 490 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1; 491 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 492 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1; 493 } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */ 494 if (n3 < 0) n3 = rank + (m - 1); 495 if (n5 < 0) n5 = rank - (m - 1); 496 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1; 497 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1; 498 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1; 499 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1; 500 } 501 502 PetscCall(PetscMalloc1(9, &dd->neighbors)); 503 504 dd->neighbors[0] = n0; 505 dd->neighbors[1] = n1; 506 dd->neighbors[2] = n2; 507 dd->neighbors[3] = n3; 508 dd->neighbors[4] = rank; 509 dd->neighbors[5] = n5; 510 dd->neighbors[6] = n6; 511 dd->neighbors[7] = n7; 512 dd->neighbors[8] = n8; 513 514 if (stencil_type == DMDA_STENCIL_STAR) { 515 /* save corner processor numbers */ 516 sn0 = n0; 517 sn2 = n2; 518 sn6 = n6; 519 sn8 = n8; 520 n0 = n2 = n6 = n8 = -1; 521 } 522 523 PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx)); 524 525 nn = 0; 526 xbase = bases[rank]; 527 for (i = 1; i <= s_y; i++) { 528 if (n0 >= 0) { /* left below */ 529 x_t = lx[n0 % m]; 530 y_t = ly[(n0 / m)]; 531 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x; 532 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 533 } 534 535 if (n1 >= 0) { /* directly below */ 536 x_t = x; 537 y_t = ly[(n1 / m)]; 538 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t; 539 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 540 } else if (by == DM_BOUNDARY_MIRROR) { 541 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j; 542 } 543 544 if (n2 >= 0) { /* right below */ 545 x_t = lx[n2 % m]; 546 y_t = ly[(n2 / m)]; 547 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t; 548 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 549 } 550 } 551 552 for (i = 0; i < y; i++) { 553 if (n3 >= 0) { /* directly left */ 554 x_t = lx[n3 % m]; 555 /* y_t = y; */ 556 s_t = bases[n3] + (i + 1) * x_t - s_x; 557 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 558 } else if (bx == DM_BOUNDARY_MIRROR) { 559 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j; 560 } 561 562 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */ 563 564 if (n5 >= 0) { /* directly right */ 565 x_t = lx[n5 % m]; 566 /* y_t = y; */ 567 s_t = bases[n5] + (i)*x_t; 568 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 569 } else if (bx == DM_BOUNDARY_MIRROR) { 570 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j; 571 } 572 } 573 574 for (i = 1; i <= s_y; i++) { 575 if (n6 >= 0) { /* left above */ 576 x_t = lx[n6 % m]; 577 /* y_t = ly[(n6/m)]; */ 578 s_t = bases[n6] + (i)*x_t - s_x; 579 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 580 } 581 582 if (n7 >= 0) { /* directly above */ 583 x_t = x; 584 /* y_t = ly[(n7/m)]; */ 585 s_t = bases[n7] + (i - 1) * x_t; 586 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 587 } else if (by == DM_BOUNDARY_MIRROR) { 588 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j; 589 } 590 591 if (n8 >= 0) { /* right above */ 592 x_t = lx[n8 % m]; 593 /* y_t = ly[(n8/m)]; */ 594 s_t = bases[n8] + (i - 1) * x_t; 595 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 596 } 597 } 598 599 PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from)); 600 PetscCall(VecScatterCreate(global, from, local, to, >ol)); 601 PetscCall(ISDestroy(&to)); 602 PetscCall(ISDestroy(&from)); 603 604 if (stencil_type == DMDA_STENCIL_STAR) { 605 n0 = sn0; 606 n2 = sn2; 607 n6 = sn6; 608 n8 = sn8; 609 } 610 611 if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) { 612 /* 613 Recompute the local to global mappings, this time keeping the 614 information about the cross corner processor numbers and any ghosted 615 but not periodic indices. 616 */ 617 nn = 0; 618 xbase = bases[rank]; 619 for (i = 1; i <= s_y; i++) { 620 if (n0 >= 0) { /* left below */ 621 x_t = lx[n0 % m]; 622 y_t = ly[(n0 / m)]; 623 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x; 624 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 625 } else if (xs - Xs > 0 && ys - Ys > 0) { 626 for (j = 0; j < s_x; j++) idx[nn++] = -1; 627 } 628 if (n1 >= 0) { /* directly below */ 629 x_t = x; 630 y_t = ly[(n1 / m)]; 631 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t; 632 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 633 } else if (ys - Ys > 0) { 634 if (by == DM_BOUNDARY_MIRROR) { 635 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j; 636 } else { 637 for (j = 0; j < x; j++) idx[nn++] = -1; 638 } 639 } 640 if (n2 >= 0) { /* right below */ 641 x_t = lx[n2 % m]; 642 y_t = ly[(n2 / m)]; 643 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t; 644 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 645 } else if (Xe - xe > 0 && ys - Ys > 0) { 646 for (j = 0; j < s_x; j++) idx[nn++] = -1; 647 } 648 } 649 650 for (i = 0; i < y; i++) { 651 if (n3 >= 0) { /* directly left */ 652 x_t = lx[n3 % m]; 653 /* y_t = y; */ 654 s_t = bases[n3] + (i + 1) * x_t - s_x; 655 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 656 } else if (xs - Xs > 0) { 657 if (bx == DM_BOUNDARY_MIRROR) { 658 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j; 659 } else { 660 for (j = 0; j < s_x; j++) idx[nn++] = -1; 661 } 662 } 663 664 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */ 665 666 if (n5 >= 0) { /* directly right */ 667 x_t = lx[n5 % m]; 668 /* y_t = y; */ 669 s_t = bases[n5] + (i)*x_t; 670 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 671 } else if (Xe - xe > 0) { 672 if (bx == DM_BOUNDARY_MIRROR) { 673 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j; 674 } else { 675 for (j = 0; j < s_x; j++) idx[nn++] = -1; 676 } 677 } 678 } 679 680 for (i = 1; i <= s_y; i++) { 681 if (n6 >= 0) { /* left above */ 682 x_t = lx[n6 % m]; 683 /* y_t = ly[(n6/m)]; */ 684 s_t = bases[n6] + (i)*x_t - s_x; 685 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 686 } else if (xs - Xs > 0 && Ye - ye > 0) { 687 for (j = 0; j < s_x; j++) idx[nn++] = -1; 688 } 689 if (n7 >= 0) { /* directly above */ 690 x_t = x; 691 /* y_t = ly[(n7/m)]; */ 692 s_t = bases[n7] + (i - 1) * x_t; 693 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 694 } else if (Ye - ye > 0) { 695 if (by == DM_BOUNDARY_MIRROR) { 696 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j; 697 } else { 698 for (j = 0; j < x; j++) idx[nn++] = -1; 699 } 700 } 701 if (n8 >= 0) { /* right above */ 702 x_t = lx[n8 % m]; 703 /* y_t = ly[(n8/m)]; */ 704 s_t = bases[n8] + (i - 1) * x_t; 705 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 706 } else if (Xe - xe > 0 && Ye - ye > 0) { 707 for (j = 0; j < s_x; j++) idx[nn++] = -1; 708 } 709 } 710 } 711 /* 712 Set the local to global ordering in the global vector, this allows use 713 of VecSetValuesLocal(). 714 */ 715 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap)); 716 717 PetscCall(PetscFree2(bases, ldims)); 718 dd->m = m; 719 dd->n = n; 720 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 721 dd->xs = xs * dof; 722 dd->xe = xe * dof; 723 dd->ys = ys; 724 dd->ye = ye; 725 dd->zs = 0; 726 dd->ze = 1; 727 dd->Xs = Xs * dof; 728 dd->Xe = Xe * dof; 729 dd->Ys = Ys; 730 dd->Ye = Ye; 731 dd->Zs = 0; 732 dd->Ze = 1; 733 734 PetscCall(VecDestroy(&local)); 735 PetscCall(VecDestroy(&global)); 736 737 dd->gtol = gtol; 738 dd->base = base; 739 da->ops->view = DMView_DA_2d; 740 dd->ltol = NULL; 741 dd->ao = NULL; 742 PetscFunctionReturn(0); 743 } 744 745 /*@C 746 DMDACreate2d - Creates an object that will manage the communication of two-dimensional 747 regular array data that is distributed across some processors. 748 749 Collective 750 751 Input Parameters: 752 + comm - MPI communicator 753 . bx,by - type of ghost nodes the array have. 754 Use one of DM_BOUNDARY_NONE, DM_BOUNDARY_GHOSTED, DM_BOUNDARY_PERIODIC. 755 . stencil_type - stencil type. Use either DMDA_STENCIL_BOX or DMDA_STENCIL_STAR. 756 . M,N - global dimension in each direction of the array 757 . m,n - corresponding number of processors in each dimension 758 (or PETSC_DECIDE to have calculated) 759 . dof - number of degrees of freedom per node 760 . s - stencil width 761 - lx, ly - arrays containing the number of nodes in each cell along 762 the x and y coordinates, or NULL. If non-null, these 763 must be of length as m and n, and the corresponding 764 m and n cannot be PETSC_DECIDE. The sum of the lx[] entries 765 must be M, and the sum of the ly[] entries must be N. 766 767 Output Parameter: 768 . da - the resulting distributed array object 769 770 Options Database Key: 771 + -dm_view - Calls DMView() at the conclusion of DMDACreate2d() 772 . -da_grid_x <nx> - number of grid points in x direction 773 . -da_grid_y <ny> - number of grid points in y direction 774 . -da_processors_x <nx> - number of processors in x direction 775 . -da_processors_y <ny> - number of processors in y direction 776 . -da_refine_x <rx> - refinement ratio in x direction 777 . -da_refine_y <ry> - refinement ratio in y direction 778 - -da_refine <n> - refine the DMDA n times before creating 779 780 Level: beginner 781 782 Notes: 783 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 784 standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 785 the standard 9-pt stencil. 786 787 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 788 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 789 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 790 791 You must call DMSetUp() after this call before using this DM. 792 793 If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call 794 but before DMSetUp(). 795 796 .seealso: `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, 797 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, 798 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 799 `DMStagCreate2d()` 800 801 @*/ 802 803 PetscErrorCode DMDACreate2d(MPI_Comm comm, DMBoundaryType bx, DMBoundaryType by, DMDAStencilType stencil_type, PetscInt M, PetscInt N, PetscInt m, PetscInt n, PetscInt dof, PetscInt s, const PetscInt lx[], const PetscInt ly[], DM *da) { 804 PetscFunctionBegin; 805 PetscCall(DMDACreate(comm, da)); 806 PetscCall(DMSetDimension(*da, 2)); 807 PetscCall(DMDASetSizes(*da, M, N, 1)); 808 PetscCall(DMDASetNumProcs(*da, m, n, PETSC_DECIDE)); 809 PetscCall(DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE)); 810 PetscCall(DMDASetDof(*da, dof)); 811 PetscCall(DMDASetStencilType(*da, stencil_type)); 812 PetscCall(DMDASetStencilWidth(*da, s)); 813 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, NULL)); 814 PetscFunctionReturn(0); 815 } 816