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          isascii, isstring;

 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: }

 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) 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;

 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));

 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 %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: }

 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:   DMDAMapCoordsToPeriodicDomain - Maps a `(x, y)` coordinate pair that lies outside a 2D `DMDA` domain back into the domain along any periodic boundaries

 90:   Not Collective

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

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

 99:   Level: developer

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

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

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

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