1 #include <../src/ts/characteristic/impls/da/slda.h> /*I "petsccharacteristic.h" I*/
2 #include <petscdmda.h>
3 #include <petscviewer.h>
4
CharacteristicView_DA(Characteristic c,PetscViewer viewer)5 static PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer)
6 {
7 Characteristic_DA *da = (Characteristic_DA *)c->data;
8 PetscBool isascii, isstring;
9
10 PetscFunctionBegin;
11 /* Pull out field names from DM */
12 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
13 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
14 if (isascii) {
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 }
21
CharacteristicDestroy_DA(Characteristic c)22 static PetscErrorCode CharacteristicDestroy_DA(Characteristic c)
23 {
24 Characteristic_DA *da = (Characteristic_DA *)c->data;
25
26 PetscFunctionBegin;
27 PetscCall(PetscFree(da));
28 PetscFunctionReturn(PETSC_SUCCESS);
29 }
30
CharacteristicSetUp_DA(Characteristic c)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;
37
38 PetscCall(DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
39 if (c->structured) PetscCall(PetscMPIIntCast(dim, &c->numIds));
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;
43
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 PetscCall(PetscMPIIntCast(numValues, &blockLen[1]));
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));
53
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;
58
59 /* Allocate communication structures */
60 PetscCheck(c->numNeighbors > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. 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 }
69
CharacteristicCreate_DA(Characteristic c)70 PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
71 {
72 Characteristic_DA *da;
73
74 PetscFunctionBegin;
75 PetscCall(PetscNew(&da));
76 PetscCall(PetscMemzero(da, sizeof(Characteristic_DA)));
77 c->data = (void *)da;
78
79 c->ops->setup = CharacteristicSetUp_DA;
80 c->ops->destroy = CharacteristicDestroy_DA;
81 c->ops->view = CharacteristicView_DA;
82
83 da->dummy = 0;
84 PetscFunctionReturn(PETSC_SUCCESS);
85 }
86
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 */
DMDAMapCoordsToPeriodicDomain(DM da,PetscScalar * x,PetscScalar * y)91 PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
92 {
93 DMBoundaryType bx, by;
94 PetscInt dim, gx, gy;
95
96 PetscFunctionBegin;
97 PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL));
98
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 }
109