Actual source code: slda.c
1: #include <../src/ts/characteristic/impls/da/slda.h>
2: #include <petscdmda.h>
3: #include <petscviewer.h>
5: static PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer)
6: {
7: Characteristic_DA *da = (Characteristic_DA *)c->data;
8: PetscBool iascii, isstring;
10: PetscFunctionBegin;
11: /* Pull out field names from DM */
12: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
13: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
14: if (iascii) {
15: PetscCall(PetscViewerASCIIPrintf(viewer, " DMDA: dummy=%" PetscInt_FMT "\n", da->dummy));
16: } else if (isstring) {
17: PetscCall(PetscViewerStringSPrintf(viewer, "dummy %" PetscInt_FMT, da->dummy));
18: }
19: PetscFunctionReturn(PETSC_SUCCESS);
20: }
22: static PetscErrorCode CharacteristicDestroy_DA(Characteristic c)
23: {
24: Characteristic_DA *da = (Characteristic_DA *)c->data;
26: PetscFunctionBegin;
27: PetscCall(PetscFree(da));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: static PetscErrorCode CharacteristicSetUp_DA(Characteristic c)
32: {
33: PetscMPIInt blockLen[2];
34: MPI_Aint indices[2];
35: MPI_Datatype oldtypes[2];
36: PetscInt dim, numValues;
38: PetscCall(DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
39: if (c->structured) c->numIds = dim;
40: else c->numIds = 3;
41: PetscCheck(c->numFieldComp <= MAX_COMPONENTS, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "The maximum number of fields allowed is %d, you have %" PetscInt_FMT ". You can recompile after increasing MAX_COMPONENTS.", MAX_COMPONENTS, c->numFieldComp);
42: numValues = 4 + MAX_COMPONENTS;
44: /* Create new MPI datatype for communication of characteristic point structs */
45: blockLen[0] = 1 + c->numIds;
46: indices[0] = 0;
47: oldtypes[0] = MPIU_INT;
48: blockLen[1] = numValues;
49: indices[1] = (1 + c->numIds) * sizeof(PetscInt);
50: oldtypes[1] = MPIU_SCALAR;
51: PetscCallMPI(MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType));
52: PetscCallMPI(MPI_Type_commit(&c->itemType));
54: /* Initialize the local queue for char foot values */
55: PetscCall(VecGetLocalSize(c->velocity, &c->queueMax));
56: PetscCall(PetscMalloc1(c->queueMax, &c->queue));
57: c->queueSize = 0;
59: /* Allocate communication structures */
60: PetscCheck(c->numNeighbors > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %" PetscInt_FMT ". Call CharactersiticSetNeighbors() before setup.", c->numNeighbors);
61: PetscCall(PetscMalloc1(c->numNeighbors, &c->needCount));
62: PetscCall(PetscMalloc1(c->numNeighbors, &c->localOffsets));
63: PetscCall(PetscMalloc1(c->numNeighbors, &c->fillCount));
64: PetscCall(PetscMalloc1(c->numNeighbors, &c->remoteOffsets));
65: PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->request));
66: PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->status));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
71: {
72: Characteristic_DA *da;
74: PetscFunctionBegin;
75: PetscCall(PetscNew(&da));
76: PetscCall(PetscMemzero(da, sizeof(Characteristic_DA)));
77: c->data = (void *)da;
79: c->ops->setup = CharacteristicSetUp_DA;
80: c->ops->destroy = CharacteristicDestroy_DA;
81: c->ops->view = CharacteristicView_DA;
83: da->dummy = 0;
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: /* -----------------------------------------------------------------------------
88: Checks for periodicity of a DM and Maps points outside of a domain back onto the domain
89: using appropriate periodicity. At the moment assumes only a 2-D DMDA.
90: ----------------------------------------------------------------------------------------*/
91: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
92: {
93: DMBoundaryType bx, by;
94: PetscInt dim, gx, gy;
96: PetscFunctionBegin;
97: PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL));
99: if (bx == DM_BOUNDARY_PERIODIC) {
100: while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
101: while (*x < 0.0) *x += (PetscScalar)gx;
102: }
103: if (by == DM_BOUNDARY_PERIODIC) {
104: while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
105: while (*y < 0.0) *y += (PetscScalar)gy;
106: }
107: PetscFunctionReturn(PETSC_SUCCESS);
108: }