1: #include <../src/ts/characteristic/impls/da/slda.h>
  2: #include <petscdmda.h>

  4: static PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer)
  5: {
  6:   Characteristic_DA *da = (Characteristic_DA *)c->data;
  7:   PetscBool          isascii, isstring;

  9:   PetscFunctionBegin;
 10:   /* Pull out field names from DM */
 11:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
 12:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
 13:   if (isascii) {
 14:     PetscCall(PetscViewerASCIIPrintf(viewer, "  DMDA: dummy=%" PetscInt_FMT "\n", da->dummy));
 15:   } else if (isstring) {
 16:     PetscCall(PetscViewerStringSPrintf(viewer, "dummy %" PetscInt_FMT, da->dummy));
 17:   }
 18:   PetscFunctionReturn(PETSC_SUCCESS);
 19: }

 21: static PetscErrorCode CharacteristicDestroy_DA(Characteristic c)
 22: {
 23:   Characteristic_DA *da = (Characteristic_DA *)c->data;

 25:   PetscFunctionBegin;
 26:   PetscCall(PetscFree(da));
 27:   PetscFunctionReturn(PETSC_SUCCESS);
 28: }

 30: static PetscErrorCode CharacteristicSetUp_DA(Characteristic c)
 31: {
 32:   PetscMPIInt  blockLen[2];
 33:   MPI_Aint     indices[2];
 34:   MPI_Datatype oldtypes[2];
 35:   PetscInt     dim, numValues;

 37:   PetscCall(DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL));
 38:   if (c->structured) PetscCall(PetscMPIIntCast(dim, &c->numIds));
 39:   else c->numIds = 3;
 40:   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);
 41:   numValues = 4 + MAX_COMPONENTS;

 43:   /* Create new MPI datatype for communication of characteristic point structs */
 44:   blockLen[0] = 1 + c->numIds;
 45:   indices[0]  = 0;
 46:   oldtypes[0] = MPIU_INT;
 47:   PetscCall(PetscMPIIntCast(numValues, &blockLen[1]));
 48:   indices[1]  = (1 + c->numIds) * sizeof(PetscInt);
 49:   oldtypes[1] = MPIU_SCALAR;
 50:   PetscCallMPI(MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType));
 51:   PetscCallMPI(MPI_Type_commit(&c->itemType));

 53:   /* Initialize the local queue for char foot values */
 54:   PetscCall(VecGetLocalSize(c->velocity, &c->queueMax));
 55:   PetscCall(PetscMalloc1(c->queueMax, &c->queue));
 56:   c->queueSize = 0;

 58:   /* Allocate communication structures */
 59:   PetscCheck(c->numNeighbors > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. Call CharactersiticSetNeighbors() before setup.", c->numNeighbors);
 60:   PetscCall(PetscMalloc1(c->numNeighbors, &c->needCount));
 61:   PetscCall(PetscMalloc1(c->numNeighbors, &c->localOffsets));
 62:   PetscCall(PetscMalloc1(c->numNeighbors, &c->fillCount));
 63:   PetscCall(PetscMalloc1(c->numNeighbors, &c->remoteOffsets));
 64:   PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->request));
 65:   PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->status));
 66:   PetscFunctionReturn(PETSC_SUCCESS);
 67: }

 69: PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c)
 70: {
 71:   Characteristic_DA *da;

 73:   PetscFunctionBegin;
 74:   PetscCall(PetscNew(&da));
 75:   PetscCall(PetscMemzero(da, sizeof(Characteristic_DA)));
 76:   c->data = (void *)da;

 78:   c->ops->setup   = CharacteristicSetUp_DA;
 79:   c->ops->destroy = CharacteristicDestroy_DA;
 80:   c->ops->view    = CharacteristicView_DA;

 82:   da->dummy = 0;
 83:   PetscFunctionReturn(PETSC_SUCCESS);
 84: }

 86: /*@
 87:   DMDAMapCoordsToPeriodicDomain - Maps a `(x, y)` coordinate pair that lies outside a 2D `DMDA` domain back into the domain along any periodic boundaries

 89:   Not Collective

 91:   Input Parameter:
 92: . da - the 2D `DMDA` context

 94:   Input/Output Parameters:
 95: + x - the x coordinate; wrapped in place if the x boundary is periodic
 96: - y - the y coordinate; wrapped in place if the y boundary is periodic

 98:   Level: developer

100:   Note:
101:   Coordinates along non-periodic boundaries are left unchanged. This helper is used internally by the method-of-characteristics solver on structured grids.

103: .seealso: [](ch_ts), `Characteristic`, `DMDA`, `CharacteristicSolve()`
104: @*/
105: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
106: {
107:   DMBoundaryType bx, by;
108:   PetscInt       dim, gx, gy;

110:   PetscFunctionBegin;
111:   PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL));

113:   if (bx == DM_BOUNDARY_PERIODIC) {
114:     while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
115:     while (*x < 0.0) *x += (PetscScalar)gx;
116:   }
117:   if (by == DM_BOUNDARY_PERIODIC) {
118:     while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
119:     while (*y < 0.0) *y += (PetscScalar)gy;
120:   }
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }