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