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=%" PetscInt_FMT "\n", da->dummy);
 15:   } else if (isstring) {
 16:     PetscViewerStringSPrintf(viewer, "dummy %" PetscInt_FMT, 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;
 44:   indices[0]  = 0;
 45:   oldtypes[0] = MPIU_INT;
 46:   blockLen[1] = numValues;
 47:   indices[1]  = (1 + c->numIds) * sizeof(PetscInt);
 48:   oldtypes[1] = MPIU_SCALAR;
 49:   MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);
 50:   MPI_Type_commit(&c->itemType);

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

 57:   /* Allocate communication structures */
 59:   PetscMalloc1(c->numNeighbors, &c->needCount);
 60:   PetscMalloc1(c->numNeighbors, &c->localOffsets);
 61:   PetscMalloc1(c->numNeighbors, &c->fillCount);
 62:   PetscMalloc1(c->numNeighbors, &c->remoteOffsets);
 63:   PetscMalloc1(c->numNeighbors - 1, &c->request);
 64:   PetscMalloc1(c->numNeighbors - 1, &c->status);
 65:   return 0;
 66: }

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

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

 76:   c->ops->setup   = CharacteristicSetUp_DA;
 77:   c->ops->destroy = CharacteristicDestroy_DA;
 78:   c->ops->view    = CharacteristicView_DA;

 80:   da->dummy = 0;
 81:   return 0;
 82: }

 84: /* -----------------------------------------------------------------------------
 85:    Checks for periodicity of a DM and Maps points outside of a domain back onto the domain
 86:    using appropriate periodicity. At the moment assumes only a 2-D DMDA.
 87:    ----------------------------------------------------------------------------------------*/
 88: PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y)
 89: {
 90:   DMBoundaryType bx, by;
 91:   PetscInt       dim, gx, gy;

 93:   DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);

 95:   if (bx == DM_BOUNDARY_PERIODIC) {
 96:     while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx;
 97:     while (*x < 0.0) *x += (PetscScalar)gx;
 98:   }
 99:   if (by == DM_BOUNDARY_PERIODIC) {
100:     while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy;
101:     while (*y < 0.0) *y += (PetscScalar)gy;
102:   }
103:   return 0;
104: }