1 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 2 #include <petscdraw.h> 3 4 static PetscErrorCode DMView_DA_2d(DM da, PetscViewer viewer) 5 { 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 { 156 PetscScalar h, *aa, *ww, v; 157 PetscReal epsilon = PETSC_SQRT_MACHINE_EPSILON, umin = 100.0 * PETSC_SQRT_MACHINE_EPSILON; 158 PetscInt gI, nI; 159 MatStencil stencil; 160 DMDALocalInfo info; 161 162 PetscFunctionBegin; 163 PetscCall((*ctx->func)(0, U, a, ctx->funcctx)); 164 PetscCall((*ctx->funcisetbase)(U, ctx->funcctx)); 165 166 PetscCall(VecGetArray(U, &ww)); 167 PetscCall(VecGetArray(a, &aa)); 168 169 nI = 0; 170 h = ww[gI]; 171 if (h == 0.0) h = 1.0; 172 if (PetscAbsScalar(h) < umin && PetscRealPart(h) >= 0.0) h = umin; 173 else if (PetscRealPart(h) < 0.0 && PetscAbsScalar(h) < umin) h = -umin; 174 h *= epsilon; 175 176 ww[gI] += h; 177 PetscCall((*ctx->funci)(i, w, &v, ctx->funcctx)); 178 aa[nI] = (v - aa[nI]) / h; 179 ww[gI] -= h; 180 nI++; 181 182 PetscCall(VecRestoreArray(U, &ww)); 183 PetscCall(VecRestoreArray(a, &aa)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 #endif 187 188 PetscErrorCode DMSetUp_DA_2D(DM da) 189 { 190 DM_DA *dd = (DM_DA *)da->data; 191 const PetscInt M = dd->M; 192 const PetscInt N = dd->N; 193 PetscInt m = dd->m; 194 PetscInt n = dd->n; 195 const PetscInt dof = dd->w; 196 const PetscInt s = dd->s; 197 DMBoundaryType bx = dd->bx; 198 DMBoundaryType by = dd->by; 199 DMDAStencilType stencil_type = dd->stencil_type; 200 PetscInt *lx = dd->lx; 201 PetscInt *ly = dd->ly; 202 MPI_Comm comm; 203 PetscMPIInt rank, size; 204 PetscInt xs, xe, ys, ye, x, y, Xs, Xe, Ys, Ye, IXs, IXe, IYs, IYe; 205 PetscInt up, down, left, right, i, n0, n1, n2, n3, n5, n6, n7, n8, *idx, nn; 206 PetscInt xbase, *bases, *ldims, j, x_t, y_t, s_t, base, count; 207 PetscInt s_x, s_y; /* s proportionalized to w */ 208 PetscInt sn0 = 0, sn2 = 0, sn6 = 0, sn8 = 0; 209 Vec local, global; 210 VecScatter gtol; 211 IS to, from; 212 213 PetscFunctionBegin; 214 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"); 215 PetscCall(PetscObjectGetComm((PetscObject)da, &comm)); 216 #if !defined(PETSC_USE_64BIT_INDICES) 217 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); 218 #endif 219 220 PetscCallMPI(MPI_Comm_size(comm, &size)); 221 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 222 223 dd->p = 1; 224 if (m != PETSC_DECIDE) { 225 PetscCheck(m >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in X direction: %" PetscInt_FMT, m); 226 PetscCheck(m <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in X direction: %" PetscInt_FMT " %d", m, size); 227 } 228 if (n != PETSC_DECIDE) { 229 PetscCheck(n >= 1, comm, PETSC_ERR_ARG_OUTOFRANGE, "Non-positive number of processors in Y direction: %" PetscInt_FMT, n); 230 PetscCheck(n <= size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Too many processors in Y direction: %" PetscInt_FMT " %d", n, size); 231 } 232 233 if (m == PETSC_DECIDE || n == PETSC_DECIDE) { 234 if (n != PETSC_DECIDE) { 235 m = size / n; 236 } else if (m != PETSC_DECIDE) { 237 n = size / m; 238 } else { 239 /* try for squarish distribution */ 240 m = (PetscInt)(0.5 + PetscSqrtReal(((PetscReal)M) * ((PetscReal)size) / ((PetscReal)N))); 241 if (!m) m = 1; 242 while (m > 0) { 243 n = size / m; 244 if (m * n == size) break; 245 m--; 246 } 247 if (M > N && m < n) { 248 PetscInt _m = m; 249 m = n; 250 n = _m; 251 } 252 } 253 PetscCheck(m * n == size, comm, PETSC_ERR_PLIB, "Unable to create partition, check the size of the communicator and input m and n "); 254 } else PetscCheck(m * n == size, comm, PETSC_ERR_ARG_OUTOFRANGE, "Given Bad partition"); 255 256 PetscCheck(M >= m, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in x direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, M, m); 257 PetscCheck(N >= n, comm, PETSC_ERR_ARG_OUTOFRANGE, "Partition in y direction is too fine! %" PetscInt_FMT " %" PetscInt_FMT, N, n); 258 259 /* 260 Determine locally owned region 261 xs is the first local node number, x is the number of local nodes 262 */ 263 if (!lx) { 264 PetscCall(PetscMalloc1(m, &dd->lx)); 265 lx = dd->lx; 266 for (i = 0; i < m; i++) lx[i] = M / m + ((M % m) > i); 267 } 268 x = lx[rank % m]; 269 xs = 0; 270 for (i = 0; i < (rank % m); i++) xs += lx[i]; 271 if (PetscDefined(USE_DEBUG)) { 272 left = xs; 273 for (i = (rank % m); i < m; i++) left += lx[i]; 274 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); 275 } 276 277 /* 278 Determine locally owned region 279 ys is the first local node number, y is the number of local nodes 280 */ 281 if (!ly) { 282 PetscCall(PetscMalloc1(n, &dd->ly)); 283 ly = dd->ly; 284 for (i = 0; i < n; i++) ly[i] = N / n + ((N % n) > i); 285 } 286 y = ly[rank / m]; 287 ys = 0; 288 for (i = 0; i < (rank / m); i++) ys += ly[i]; 289 if (PetscDefined(USE_DEBUG)) { 290 left = ys; 291 for (i = (rank / m); i < n; i++) left += ly[i]; 292 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); 293 } 294 295 /* 296 check if the scatter requires more than one process neighbor or wraps around 297 the domain more than once 298 */ 299 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); 300 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); 301 xe = xs + x; 302 ye = ys + y; 303 304 /* determine ghost region (Xs) and region scattered into (IXs) */ 305 if (xs - s > 0) { 306 Xs = xs - s; 307 IXs = xs - s; 308 } else { 309 if (bx) { 310 Xs = xs - s; 311 } else { 312 Xs = 0; 313 } 314 IXs = 0; 315 } 316 if (xe + s <= M) { 317 Xe = xe + s; 318 IXe = xe + s; 319 } else { 320 if (bx) { 321 Xs = xs - s; 322 Xe = xe + s; 323 } else { 324 Xe = M; 325 } 326 IXe = M; 327 } 328 329 if (bx == DM_BOUNDARY_PERIODIC || bx == DM_BOUNDARY_MIRROR) { 330 IXs = xs - s; 331 IXe = xe + s; 332 Xs = xs - s; 333 Xe = xe + s; 334 } 335 336 if (ys - s > 0) { 337 Ys = ys - s; 338 IYs = ys - s; 339 } else { 340 if (by) { 341 Ys = ys - s; 342 } else { 343 Ys = 0; 344 } 345 IYs = 0; 346 } 347 if (ye + s <= N) { 348 Ye = ye + s; 349 IYe = ye + s; 350 } else { 351 if (by) { 352 Ye = ye + s; 353 } else { 354 Ye = N; 355 } 356 IYe = N; 357 } 358 359 if (by == DM_BOUNDARY_PERIODIC || by == DM_BOUNDARY_MIRROR) { 360 IYs = ys - s; 361 IYe = ye + s; 362 Ys = ys - s; 363 Ye = ye + s; 364 } 365 366 /* stencil length in each direction */ 367 s_x = s; 368 s_y = s; 369 370 /* determine starting point of each processor */ 371 nn = x * y; 372 PetscCall(PetscMalloc2(size + 1, &bases, size, &ldims)); 373 PetscCallMPI(MPI_Allgather(&nn, 1, MPIU_INT, ldims, 1, MPIU_INT, comm)); 374 bases[0] = 0; 375 for (i = 1; i <= size; i++) bases[i] = ldims[i - 1]; 376 for (i = 1; i <= size; i++) bases[i] += bases[i - 1]; 377 base = bases[rank] * dof; 378 379 /* allocate the base parallel and sequential vectors */ 380 dd->Nlocal = x * y * dof; 381 PetscCall(VecCreateMPIWithArray(comm, dof, dd->Nlocal, PETSC_DECIDE, NULL, &global)); 382 dd->nlocal = (Xe - Xs) * (Ye - Ys) * dof; 383 PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, dof, dd->nlocal, NULL, &local)); 384 385 /* generate global to local vector scatter and local to global mapping*/ 386 387 /* global to local must include ghost points within the domain, 388 but not ghost points outside the domain that aren't periodic */ 389 PetscCall(PetscMalloc1((IXe - IXs) * (IYe - IYs), &idx)); 390 if (stencil_type == DMDA_STENCIL_BOX) { 391 left = IXs - Xs; 392 right = left + (IXe - IXs); 393 down = IYs - Ys; 394 up = down + (IYe - IYs); 395 count = 0; 396 for (i = down; i < up; i++) { 397 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 398 } 399 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 400 401 } else { 402 /* must drop into cross shape region */ 403 /* ---------| 404 | top | 405 |--- ---| up 406 | middle | 407 | | 408 ---- ---- down 409 | bottom | 410 ----------- 411 Xs xs xe Xe */ 412 left = xs - Xs; 413 right = left + x; 414 down = ys - Ys; 415 up = down + y; 416 count = 0; 417 /* bottom */ 418 for (i = (IYs - Ys); i < down; i++) { 419 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 420 } 421 /* middle */ 422 for (i = down; i < up; i++) { 423 for (j = (IXs - Xs); j < (IXe - Xs); j++) idx[count++] = j + i * (Xe - Xs); 424 } 425 /* top */ 426 for (i = up; i < up + IYe - ye; i++) { 427 for (j = left; j < right; j++) idx[count++] = j + i * (Xe - Xs); 428 } 429 PetscCall(ISCreateBlock(comm, dof, count, idx, PETSC_OWN_POINTER, &to)); 430 } 431 432 /* determine who lies on each side of us stored in n6 n7 n8 433 n3 n5 434 n0 n1 n2 435 */ 436 437 /* Assume the Non-Periodic Case */ 438 n1 = rank - m; 439 if (rank % m) { 440 n0 = n1 - 1; 441 } else { 442 n0 = -1; 443 } 444 if ((rank + 1) % m) { 445 n2 = n1 + 1; 446 n5 = rank + 1; 447 n8 = rank + m + 1; 448 if (n8 >= m * n) n8 = -1; 449 } else { 450 n2 = -1; 451 n5 = -1; 452 n8 = -1; 453 } 454 if (rank % m) { 455 n3 = rank - 1; 456 n6 = n3 + m; 457 if (n6 >= m * n) n6 = -1; 458 } else { 459 n3 = -1; 460 n6 = -1; 461 } 462 n7 = rank + m; 463 if (n7 >= m * n) n7 = -1; 464 465 if (bx == DM_BOUNDARY_PERIODIC && by == DM_BOUNDARY_PERIODIC) { 466 /* Modify for Periodic Cases */ 467 /* Handle all four corners */ 468 if ((n6 < 0) && (n7 < 0) && (n3 < 0)) n6 = m - 1; 469 if ((n8 < 0) && (n7 < 0) && (n5 < 0)) n8 = 0; 470 if ((n2 < 0) && (n5 < 0) && (n1 < 0)) n2 = size - m; 471 if ((n0 < 0) && (n3 < 0) && (n1 < 0)) n0 = size - 1; 472 473 /* Handle Top and Bottom Sides */ 474 if (n1 < 0) n1 = rank + m * (n - 1); 475 if (n7 < 0) n7 = rank - m * (n - 1); 476 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 477 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1; 478 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 479 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1; 480 481 /* Handle Left and Right Sides */ 482 if (n3 < 0) n3 = rank + (m - 1); 483 if (n5 < 0) n5 = rank - (m - 1); 484 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1; 485 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1; 486 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1; 487 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1; 488 } else if (by == DM_BOUNDARY_PERIODIC) { /* Handle Top and Bottom Sides */ 489 if (n1 < 0) n1 = rank + m * (n - 1); 490 if (n7 < 0) n7 = rank - m * (n - 1); 491 if ((n3 >= 0) && (n0 < 0)) n0 = size - m + rank - 1; 492 if ((n3 >= 0) && (n6 < 0)) n6 = (rank % m) - 1; 493 if ((n5 >= 0) && (n2 < 0)) n2 = size - m + rank + 1; 494 if ((n5 >= 0) && (n8 < 0)) n8 = (rank % m) + 1; 495 } else if (bx == DM_BOUNDARY_PERIODIC) { /* Handle Left and Right Sides */ 496 if (n3 < 0) n3 = rank + (m - 1); 497 if (n5 < 0) n5 = rank - (m - 1); 498 if ((n1 >= 0) && (n0 < 0)) n0 = rank - 1; 499 if ((n1 >= 0) && (n2 < 0)) n2 = rank - 2 * m + 1; 500 if ((n7 >= 0) && (n6 < 0)) n6 = rank + 2 * m - 1; 501 if ((n7 >= 0) && (n8 < 0)) n8 = rank + 1; 502 } 503 504 PetscCall(PetscMalloc1(9, &dd->neighbors)); 505 506 dd->neighbors[0] = n0; 507 dd->neighbors[1] = n1; 508 dd->neighbors[2] = n2; 509 dd->neighbors[3] = n3; 510 dd->neighbors[4] = rank; 511 dd->neighbors[5] = n5; 512 dd->neighbors[6] = n6; 513 dd->neighbors[7] = n7; 514 dd->neighbors[8] = n8; 515 516 if (stencil_type == DMDA_STENCIL_STAR) { 517 /* save corner processor numbers */ 518 sn0 = n0; 519 sn2 = n2; 520 sn6 = n6; 521 sn8 = n8; 522 n0 = n2 = n6 = n8 = -1; 523 } 524 525 PetscCall(PetscMalloc1((Xe - Xs) * (Ye - Ys), &idx)); 526 527 nn = 0; 528 xbase = bases[rank]; 529 for (i = 1; i <= s_y; i++) { 530 if (n0 >= 0) { /* left below */ 531 x_t = lx[n0 % m]; 532 y_t = ly[(n0 / m)]; 533 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x; 534 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 535 } 536 537 if (n1 >= 0) { /* directly below */ 538 x_t = x; 539 y_t = ly[(n1 / m)]; 540 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t; 541 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 542 } else if (by == DM_BOUNDARY_MIRROR) { 543 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j; 544 } 545 546 if (n2 >= 0) { /* right below */ 547 x_t = lx[n2 % m]; 548 y_t = ly[(n2 / m)]; 549 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t; 550 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 551 } 552 } 553 554 for (i = 0; i < y; i++) { 555 if (n3 >= 0) { /* directly left */ 556 x_t = lx[n3 % m]; 557 /* y_t = y; */ 558 s_t = bases[n3] + (i + 1) * x_t - s_x; 559 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 560 } else if (bx == DM_BOUNDARY_MIRROR) { 561 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j; 562 } 563 564 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */ 565 566 if (n5 >= 0) { /* directly right */ 567 x_t = lx[n5 % m]; 568 /* y_t = y; */ 569 s_t = bases[n5] + (i)*x_t; 570 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 571 } else if (bx == DM_BOUNDARY_MIRROR) { 572 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j; 573 } 574 } 575 576 for (i = 1; i <= s_y; i++) { 577 if (n6 >= 0) { /* left above */ 578 x_t = lx[n6 % m]; 579 /* y_t = ly[(n6/m)]; */ 580 s_t = bases[n6] + (i)*x_t - s_x; 581 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 582 } 583 584 if (n7 >= 0) { /* directly above */ 585 x_t = x; 586 /* y_t = ly[(n7/m)]; */ 587 s_t = bases[n7] + (i - 1) * x_t; 588 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 589 } else if (by == DM_BOUNDARY_MIRROR) { 590 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j; 591 } 592 593 if (n8 >= 0) { /* right above */ 594 x_t = lx[n8 % m]; 595 /* y_t = ly[(n8/m)]; */ 596 s_t = bases[n8] + (i - 1) * x_t; 597 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 598 } 599 } 600 601 PetscCall(ISCreateBlock(comm, dof, nn, idx, PETSC_USE_POINTER, &from)); 602 PetscCall(VecScatterCreate(global, from, local, to, >ol)); 603 PetscCall(ISDestroy(&to)); 604 PetscCall(ISDestroy(&from)); 605 606 if (stencil_type == DMDA_STENCIL_STAR) { 607 n0 = sn0; 608 n2 = sn2; 609 n6 = sn6; 610 n8 = sn8; 611 } 612 613 if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) { 614 /* 615 Recompute the local to global mappings, this time keeping the 616 information about the cross corner processor numbers and any ghosted 617 but not periodic indices. 618 */ 619 nn = 0; 620 xbase = bases[rank]; 621 for (i = 1; i <= s_y; i++) { 622 if (n0 >= 0) { /* left below */ 623 x_t = lx[n0 % m]; 624 y_t = ly[(n0 / m)]; 625 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x; 626 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 627 } else if (xs - Xs > 0 && ys - Ys > 0) { 628 for (j = 0; j < s_x; j++) idx[nn++] = -1; 629 } 630 if (n1 >= 0) { /* directly below */ 631 x_t = x; 632 y_t = ly[(n1 / m)]; 633 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t; 634 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 635 } else if (ys - Ys > 0) { 636 if (by == DM_BOUNDARY_MIRROR) { 637 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j; 638 } else { 639 for (j = 0; j < x; j++) idx[nn++] = -1; 640 } 641 } 642 if (n2 >= 0) { /* right below */ 643 x_t = lx[n2 % m]; 644 y_t = ly[(n2 / m)]; 645 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t; 646 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 647 } else if (Xe - xe > 0 && ys - Ys > 0) { 648 for (j = 0; j < s_x; j++) idx[nn++] = -1; 649 } 650 } 651 652 for (i = 0; i < y; i++) { 653 if (n3 >= 0) { /* directly left */ 654 x_t = lx[n3 % m]; 655 /* y_t = y; */ 656 s_t = bases[n3] + (i + 1) * x_t - s_x; 657 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 658 } else if (xs - Xs > 0) { 659 if (bx == DM_BOUNDARY_MIRROR) { 660 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j; 661 } else { 662 for (j = 0; j < s_x; j++) idx[nn++] = -1; 663 } 664 } 665 666 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */ 667 668 if (n5 >= 0) { /* directly right */ 669 x_t = lx[n5 % m]; 670 /* y_t = y; */ 671 s_t = bases[n5] + (i)*x_t; 672 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 673 } else if (Xe - xe > 0) { 674 if (bx == DM_BOUNDARY_MIRROR) { 675 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j; 676 } else { 677 for (j = 0; j < s_x; j++) idx[nn++] = -1; 678 } 679 } 680 } 681 682 for (i = 1; i <= s_y; i++) { 683 if (n6 >= 0) { /* left above */ 684 x_t = lx[n6 % m]; 685 /* y_t = ly[(n6/m)]; */ 686 s_t = bases[n6] + (i)*x_t - s_x; 687 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 688 } else if (xs - Xs > 0 && Ye - ye > 0) { 689 for (j = 0; j < s_x; j++) idx[nn++] = -1; 690 } 691 if (n7 >= 0) { /* directly above */ 692 x_t = x; 693 /* y_t = ly[(n7/m)]; */ 694 s_t = bases[n7] + (i - 1) * x_t; 695 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 696 } else if (Ye - ye > 0) { 697 if (by == DM_BOUNDARY_MIRROR) { 698 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j; 699 } else { 700 for (j = 0; j < x; j++) idx[nn++] = -1; 701 } 702 } 703 if (n8 >= 0) { /* right above */ 704 x_t = lx[n8 % m]; 705 /* y_t = ly[(n8/m)]; */ 706 s_t = bases[n8] + (i - 1) * x_t; 707 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 708 } else if (Xe - xe > 0 && Ye - ye > 0) { 709 for (j = 0; j < s_x; j++) idx[nn++] = -1; 710 } 711 } 712 } 713 /* 714 Set the local to global ordering in the global vector, this allows use 715 of VecSetValuesLocal(). 716 */ 717 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap)); 718 719 PetscCall(PetscFree2(bases, ldims)); 720 dd->m = m; 721 dd->n = n; 722 /* note petsc expects xs/xe/Xs/Xe to be multiplied by #dofs in many places */ 723 dd->xs = xs * dof; 724 dd->xe = xe * dof; 725 dd->ys = ys; 726 dd->ye = ye; 727 dd->zs = 0; 728 dd->ze = 1; 729 dd->Xs = Xs * dof; 730 dd->Xe = Xe * dof; 731 dd->Ys = Ys; 732 dd->Ye = Ye; 733 dd->Zs = 0; 734 dd->Ze = 1; 735 736 PetscCall(VecDestroy(&local)); 737 PetscCall(VecDestroy(&global)); 738 739 dd->gtol = gtol; 740 dd->base = base; 741 da->ops->view = DMView_DA_2d; 742 dd->ltol = NULL; 743 dd->ao = NULL; 744 PetscFunctionReturn(PETSC_SUCCESS); 745 } 746 747 /*@C 748 DMDACreate2d - Creates an object that will manage the communication of two-dimensional 749 regular array data that is distributed across some processors. 750 751 Collective 752 753 Input Parameters: 754 + comm - MPI communicator 755 . bx - type of ghost nodes the x array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 756 . by - type of ghost nodes the y array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 757 . stencil_type - stencil type. Use either `DMDA_STENCIL_BOX` or `DMDA_STENCIL_STAR`. 758 . M - global dimension in x direction of the array 759 . N - global dimension in y direction of the array 760 . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated) 761 . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated) 762 . dof - number of degrees of freedom per node 763 . s - stencil width 764 . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`. 765 - ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`. 766 767 Output Parameter: 768 . da - the resulting distributed array object 769 770 Options Database Keys: 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 If `lx` or `ly` are non-null, these must be of length as `m` and `n`, and the corresponding 784 `m` and `n` cannot be `PETSC_DECIDE`. The sum of the `lx` entries must be `M`, and the sum of 785 the `ly` entries must be `N`. 786 787 The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the 788 standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes 789 the standard 9-pt stencil. 790 791 The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects; 792 The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()` 793 and DMCreateLocalVector() and calls to `VecDuplicate()` if more are needed. 794 795 You must call `DMSetUp()` after this call before using this `DM`. 796 797 If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call 798 but before `DMSetUp()`. 799 800 .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, 801 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, 802 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 803 `DMStagCreate2d()` 804 @*/ 805 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) 806 { 807 PetscFunctionBegin; 808 PetscCall(DMDACreate(comm, da)); 809 PetscCall(DMSetDimension(*da, 2)); 810 PetscCall(DMDASetSizes(*da, M, N, 1)); 811 PetscCall(DMDASetNumProcs(*da, m, n, PETSC_DECIDE)); 812 PetscCall(DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE)); 813 PetscCall(DMDASetDof(*da, dof)); 814 PetscCall(DMDASetStencilType(*da, stencil_type)); 815 PetscCall(DMDASetStencilWidth(*da, s)); 816 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, NULL)); 817 PetscFunctionReturn(PETSC_SUCCESS); 818 } 819