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=%D\n", da->dummy);
15: } else if (isstring) {
16: PetscViewerStringSPrintf(viewer,"dummy %D", 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; indices[0] = 0; oldtypes[0] = MPIU_INT;
44: blockLen[1] = numValues; indices[1] = (1+c->numIds)*sizeof(PetscInt); oldtypes[1] = MPIU_SCALAR;
45: MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);
46: MPI_Type_commit(&c->itemType);
48: /* Initialize the local queue for char foot values */
49: VecGetLocalSize(c->velocity, &c->queueMax);
50: PetscMalloc1(c->queueMax, &c->queue);
51: c->queueSize = 0;
53: /* Allocate communication structures */
55: PetscMalloc1(c->numNeighbors, &c->needCount);
56: PetscMalloc1(c->numNeighbors, &c->localOffsets);
57: PetscMalloc1(c->numNeighbors, &c->fillCount);
58: PetscMalloc1(c->numNeighbors, &c->remoteOffsets);
59: PetscMalloc1(c->numNeighbors-1, &c->request);
60: PetscMalloc1(c->numNeighbors-1, &c->status);
61: return 0;
62: }
64: PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
65: {
66: Characteristic_DA *da;
68: PetscNew(&da);
69: PetscMemzero(da, sizeof(Characteristic_DA));
70: PetscLogObjectMemory((PetscObject)c, sizeof(Characteristic_DA));
71: c->data = (void*) da;
73: c->ops->setup = CharacteristicSetUp_DA;
74: c->ops->destroy = CharacteristicDestroy_DA;
75: c->ops->view = CharacteristicView_DA;
77: da->dummy = 0;
78: return 0;
79: }
81: /* -----------------------------------------------------------------------------
82: Checks for periodicity of a DM and Maps points outside of a domain back onto the domain
83: using appropriate periodicity. At the moment assumes only a 2-D DMDA.
84: ----------------------------------------------------------------------------------------*/
85: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
86: {
87: DMBoundaryType bx, by;
88: PetscInt dim, gx, gy;
90: DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);
92: if (bx == DM_BOUNDARY_PERIODIC) {
93: while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
94: while (*x < 0.0) *x += (PetscScalar)gx;
95: }
96: if (by == DM_BOUNDARY_PERIODIC) {
97: while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
98: while (*y < 0.0) *y += (PetscScalar)gy;
99: }
100: return 0;
101: }