Actual source code: slda.c
1: #include <../src/ts/characteristic/impls/da/slda.h>
2: #include <petscdmda.h>
3: #include <petscviewer.h>
5: PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer)
6: {
7: Characteristic_DA *da = (Characteristic_DA *)c->data;
8: PetscBool iascii, isstring;
10: /* Pull out field names from DM */
11: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
12: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
13: if (iascii) {
14: PetscViewerASCIIPrintf(viewer, " DMDA: dummy=%" PetscInt_FMT "\n", da->dummy);
15: } else if (isstring) {
16: PetscViewerStringSPrintf(viewer, "dummy %" PetscInt_FMT, da->dummy);
17: }
18: return 0;
19: }
21: PetscErrorCode CharacteristicDestroy_DA(Characteristic c)
22: {
23: Characteristic_DA *da = (Characteristic_DA *)c->data;
25: PetscFree(da);
26: return 0;
27: }
29: PetscErrorCode CharacteristicSetUp_DA(Characteristic c)
30: {
31: PetscMPIInt blockLen[2];
32: MPI_Aint indices[2];
33: MPI_Datatype oldtypes[2];
34: PetscInt dim, numValues;
36: DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);
37: if (c->structured) c->numIds = dim;
38: else c->numIds = 3;
40: numValues = 4 + MAX_COMPONENTS;
42: /* Create new MPI datatype for communication of characteristic point structs */
43: blockLen[0] = 1 + c->numIds;
44: indices[0] = 0;
45: oldtypes[0] = MPIU_INT;
46: blockLen[1] = numValues;
47: indices[1] = (1 + c->numIds) * sizeof(PetscInt);
48: oldtypes[1] = MPIU_SCALAR;
49: MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);
50: MPI_Type_commit(&c->itemType);
52: /* Initialize the local queue for char foot values */
53: VecGetLocalSize(c->velocity, &c->queueMax);
54: PetscMalloc1(c->queueMax, &c->queue);
55: c->queueSize = 0;
57: /* Allocate communication structures */
59: PetscMalloc1(c->numNeighbors, &c->needCount);
60: PetscMalloc1(c->numNeighbors, &c->localOffsets);
61: PetscMalloc1(c->numNeighbors, &c->fillCount);
62: PetscMalloc1(c->numNeighbors, &c->remoteOffsets);
63: PetscMalloc1(c->numNeighbors - 1, &c->request);
64: PetscMalloc1(c->numNeighbors - 1, &c->status);
65: return 0;
66: }
68: PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
69: {
70: Characteristic_DA *da;
72: PetscNew(&da);
73: PetscMemzero(da, sizeof(Characteristic_DA));
74: c->data = (void *)da;
76: c->ops->setup = CharacteristicSetUp_DA;
77: c->ops->destroy = CharacteristicDestroy_DA;
78: c->ops->view = CharacteristicView_DA;
80: da->dummy = 0;
81: return 0;
82: }
84: /* -----------------------------------------------------------------------------
85: Checks for periodicity of a DM and Maps points outside of a domain back onto the domain
86: using appropriate periodicity. At the moment assumes only a 2-D DMDA.
87: ----------------------------------------------------------------------------------------*/
88: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
89: {
90: DMBoundaryType bx, by;
91: PetscInt dim, gx, gy;
93: DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);
95: if (bx == DM_BOUNDARY_PERIODIC) {
96: while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
97: while (*x < 0.0) *x += (PetscScalar)gx;
98: }
99: if (by == DM_BOUNDARY_PERIODIC) {
100: while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
101: while (*y < 0.0) *y += (PetscScalar)gy;
102: }
103: return 0;
104: }