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