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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 - type of ghost nodes the x array have. Use one of `DM_BOUNDARY_NONE`, `DM_BOUNDARY_GHOSTED`, `DM_BOUNDARY_PERIODIC`. 757 . by - type of ghost nodes the y array have. 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 - global dimension in x direction of the array 760 . N - global dimension in y direction of the array 761 . m - corresponding number of processors in x dimension (or `PETSC_DECIDE` to have calculated) 762 . n - corresponding number of processors in y dimension (or `PETSC_DECIDE` to have calculated) 763 . dof - number of degrees of freedom per node 764 . s - stencil width 765 . lx - arrays containing the number of nodes in each cell along the x coordinates, or `NULL`. 766 - ly - arrays containing the number of nodes in each cell along the y coordinates, or `NULL`. 767 768 Output Parameter: 769 . da - the resulting distributed array object 770 771 Options Database Keys: 772 + -dm_view - Calls `DMView()` at the conclusion of `DMDACreate2d()` 773 . -da_grid_x <nx> - number of grid points in x direction 774 . -da_grid_y <ny> - number of grid points in y direction 775 . -da_processors_x <nx> - number of processors in x direction 776 . -da_processors_y <ny> - number of processors in y direction 777 . -da_refine_x <rx> - refinement ratio in x direction 778 . -da_refine_y <ry> - refinement ratio in y direction 779 - -da_refine <n> - refine the DMDA n times before creating 780 781 Level: beginner 782 783 Notes: 784 If `lx` or `ly` are non-null, these must be of length as `m` and `n`, and the corresponding 785 `m` and `n` cannot be `PETSC_DECIDE`. The sum of the `lx` entries must be `M`, and the sum of 786 the `ly` entries must be `N`. 787 788 The stencil type `DMDA_STENCIL_STAR` with width 1 corresponds to the 789 standard 5-pt stencil, while `DMDA_STENCIL_BOX` with width 1 denotes 790 the standard 9-pt stencil. 791 792 The array data itself is NOT stored in the `DMDA`, it is stored in `Vec` objects; 793 The appropriate vector objects can be obtained with calls to `DMCreateGlobalVector()` 794 and DMCreateLocalVector() and calls to `VecDuplicate()` if more are needed. 795 796 You must call `DMSetUp()` after this call before using this `DM`. 797 798 If you wish to use the options database to change values in the `DMDA` call `DMSetFromOptions()` after this call 799 but before `DMSetUp()`. 800 801 .seealso: `DM`, `DMDA`, `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, 802 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, 803 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 804 `DMStagCreate2d()` 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(PETSC_SUCCESS); 819 } 820