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_ENGINE) 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_ENGINE) 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_ENGINE) 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(PetscLogObjectParent((PetscObject)da, (PetscObject)gtol)); 602 PetscCall(ISDestroy(&to)); 603 PetscCall(ISDestroy(&from)); 604 605 if (stencil_type == DMDA_STENCIL_STAR) { 606 n0 = sn0; 607 n2 = sn2; 608 n6 = sn6; 609 n8 = sn8; 610 } 611 612 if (((stencil_type == DMDA_STENCIL_STAR) || (bx && bx != DM_BOUNDARY_PERIODIC) || (by && by != DM_BOUNDARY_PERIODIC))) { 613 /* 614 Recompute the local to global mappings, this time keeping the 615 information about the cross corner processor numbers and any ghosted 616 but not periodic indices. 617 */ 618 nn = 0; 619 xbase = bases[rank]; 620 for (i = 1; i <= s_y; i++) { 621 if (n0 >= 0) { /* left below */ 622 x_t = lx[n0 % m]; 623 y_t = ly[(n0 / m)]; 624 s_t = bases[n0] + x_t * y_t - (s_y - i) * x_t - s_x; 625 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 626 } else if (xs - Xs > 0 && ys - Ys > 0) { 627 for (j = 0; j < s_x; j++) idx[nn++] = -1; 628 } 629 if (n1 >= 0) { /* directly below */ 630 x_t = x; 631 y_t = ly[(n1 / m)]; 632 s_t = bases[n1] + x_t * y_t - (s_y + 1 - i) * x_t; 633 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 634 } else if (ys - Ys > 0) { 635 if (by == DM_BOUNDARY_MIRROR) { 636 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (s_y - i + 1) + j; 637 } else { 638 for (j = 0; j < x; j++) idx[nn++] = -1; 639 } 640 } 641 if (n2 >= 0) { /* right below */ 642 x_t = lx[n2 % m]; 643 y_t = ly[(n2 / m)]; 644 s_t = bases[n2] + x_t * y_t - (s_y + 1 - i) * x_t; 645 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 646 } else if (Xe - xe > 0 && ys - Ys > 0) { 647 for (j = 0; j < s_x; j++) idx[nn++] = -1; 648 } 649 } 650 651 for (i = 0; i < y; i++) { 652 if (n3 >= 0) { /* directly left */ 653 x_t = lx[n3 % m]; 654 /* y_t = y; */ 655 s_t = bases[n3] + (i + 1) * x_t - s_x; 656 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 657 } else if (xs - Xs > 0) { 658 if (bx == DM_BOUNDARY_MIRROR) { 659 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * i + s_x - j; 660 } else { 661 for (j = 0; j < s_x; j++) idx[nn++] = -1; 662 } 663 } 664 665 for (j = 0; j < x; j++) idx[nn++] = xbase++; /* interior */ 666 667 if (n5 >= 0) { /* directly right */ 668 x_t = lx[n5 % m]; 669 /* y_t = y; */ 670 s_t = bases[n5] + (i)*x_t; 671 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 672 } else if (Xe - xe > 0) { 673 if (bx == DM_BOUNDARY_MIRROR) { 674 for (j = 0; j < s_x; j++) idx[nn++] = bases[rank] + x * (i + 1) - 2 - j; 675 } else { 676 for (j = 0; j < s_x; j++) idx[nn++] = -1; 677 } 678 } 679 } 680 681 for (i = 1; i <= s_y; i++) { 682 if (n6 >= 0) { /* left above */ 683 x_t = lx[n6 % m]; 684 /* y_t = ly[(n6/m)]; */ 685 s_t = bases[n6] + (i)*x_t - s_x; 686 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 687 } else if (xs - Xs > 0 && Ye - ye > 0) { 688 for (j = 0; j < s_x; j++) idx[nn++] = -1; 689 } 690 if (n7 >= 0) { /* directly above */ 691 x_t = x; 692 /* y_t = ly[(n7/m)]; */ 693 s_t = bases[n7] + (i - 1) * x_t; 694 for (j = 0; j < x_t; j++) idx[nn++] = s_t++; 695 } else if (Ye - ye > 0) { 696 if (by == DM_BOUNDARY_MIRROR) { 697 for (j = 0; j < x; j++) idx[nn++] = bases[rank] + x * (y - i - 1) + j; 698 } else { 699 for (j = 0; j < x; j++) idx[nn++] = -1; 700 } 701 } 702 if (n8 >= 0) { /* right above */ 703 x_t = lx[n8 % m]; 704 /* y_t = ly[(n8/m)]; */ 705 s_t = bases[n8] + (i - 1) * x_t; 706 for (j = 0; j < s_x; j++) idx[nn++] = s_t++; 707 } else if (Xe - xe > 0 && Ye - ye > 0) { 708 for (j = 0; j < s_x; j++) idx[nn++] = -1; 709 } 710 } 711 } 712 /* 713 Set the local to global ordering in the global vector, this allows use 714 of VecSetValuesLocal(). 715 */ 716 PetscCall(ISLocalToGlobalMappingCreate(comm, dof, nn, idx, PETSC_OWN_POINTER, &da->ltogmap)); 717 PetscCall(PetscLogObjectParent((PetscObject)da, (PetscObject)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(0); 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,by - type of ghost nodes the array have. 756 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,N - global dimension in each direction of the array 759 . m,n - corresponding number of processors in each dimension 760 (or PETSC_DECIDE to have calculated) 761 . dof - number of degrees of freedom per node 762 . s - stencil width 763 - lx, ly - arrays containing the number of nodes in each cell along 764 the x and y coordinates, or NULL. If non-null, these 765 must be of length as m and n, and the corresponding 766 m and n cannot be PETSC_DECIDE. The sum of the lx[] entries 767 must be M, and the sum of the ly[] entries must be N. 768 769 Output Parameter: 770 . da - the resulting distributed array object 771 772 Options Database Key: 773 + -dm_view - Calls DMView() at the conclusion of DMDACreate2d() 774 . -da_grid_x <nx> - number of grid points in x direction 775 . -da_grid_y <ny> - number of grid points in y direction 776 . -da_processors_x <nx> - number of processors in x direction 777 . -da_processors_y <ny> - number of processors in y direction 778 . -da_refine_x <rx> - refinement ratio in x direction 779 . -da_refine_y <ry> - refinement ratio in y direction 780 - -da_refine <n> - refine the DMDA n times before creating 781 782 Level: beginner 783 784 Notes: 785 The stencil type DMDA_STENCIL_STAR with width 1 corresponds to the 786 standard 5-pt stencil, while DMDA_STENCIL_BOX with width 1 denotes 787 the standard 9-pt stencil. 788 789 The array data itself is NOT stored in the DMDA, it is stored in Vec objects; 790 The appropriate vector objects can be obtained with calls to DMCreateGlobalVector() 791 and DMCreateLocalVector() and calls to VecDuplicate() if more are needed. 792 793 You must call DMSetUp() after this call before using this DM. 794 795 If you wish to use the options database to change values in the DMDA call DMSetFromOptions() after this call 796 but before DMSetUp(). 797 798 .seealso: `DMDestroy()`, `DMView()`, `DMDACreate1d()`, `DMDACreate3d()`, `DMGlobalToLocalBegin()`, `DMDAGetRefinementFactor()`, 799 `DMGlobalToLocalEnd()`, `DMLocalToGlobalBegin()`, `DMLocalToLocalBegin()`, `DMLocalToLocalEnd()`, `DMDASetRefinementFactor()`, 800 `DMDAGetInfo()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`, `DMDACreateNaturalVector()`, `DMLoad()`, `DMDAGetOwnershipRanges()`, 801 `DMStagCreate2d()` 802 803 @*/ 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 PetscFunctionBegin; 807 PetscCall(DMDACreate(comm, da)); 808 PetscCall(DMSetDimension(*da, 2)); 809 PetscCall(DMDASetSizes(*da, M, N, 1)); 810 PetscCall(DMDASetNumProcs(*da, m, n, PETSC_DECIDE)); 811 PetscCall(DMDASetBoundaryType(*da, bx, by, DM_BOUNDARY_NONE)); 812 PetscCall(DMDASetDof(*da, dof)); 813 PetscCall(DMDASetStencilType(*da, stencil_type)); 814 PetscCall(DMDASetStencilWidth(*da, s)); 815 PetscCall(DMDASetOwnershipRanges(*da, lx, ly, NULL)); 816 PetscFunctionReturn(0); 817 } 818