1 #include <../src/ts/characteristic/impls/da/slda.h> /*I "petsccharacteristic.h" I*/ 2 #include <petscdmda.h> 3 #include <petscviewer.h> 4 5 PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer) { 6 Characteristic_DA *da = (Characteristic_DA *)c->data; 7 PetscBool iascii, isstring; 8 9 PetscFunctionBegin; 10 /* Pull out field names from DM */ 11 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 12 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring)); 13 if (iascii) { 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(0); 19 } 20 21 PetscErrorCode CharacteristicDestroy_DA(Characteristic c) { 22 Characteristic_DA *da = (Characteristic_DA *)c->data; 23 24 PetscFunctionBegin; 25 PetscCall(PetscFree(da)); 26 PetscFunctionReturn(0); 27 } 28 29 PetscErrorCode CharacteristicSetUp_DA(Characteristic c) { 30 PetscMPIInt blockLen[2]; 31 MPI_Aint indices[2]; 32 MPI_Datatype oldtypes[2]; 33 PetscInt dim, numValues; 34 35 PetscCall(DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL)); 36 if (c->structured) c->numIds = dim; 37 else c->numIds = 3; 38 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); 39 numValues = 4 + MAX_COMPONENTS; 40 41 /* Create new MPI datatype for communication of characteristic point structs */ 42 blockLen[0] = 1 + c->numIds; 43 indices[0] = 0; 44 oldtypes[0] = MPIU_INT; 45 blockLen[1] = numValues; 46 indices[1] = (1 + c->numIds) * sizeof(PetscInt); 47 oldtypes[1] = MPIU_SCALAR; 48 PetscCallMPI(MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType)); 49 PetscCallMPI(MPI_Type_commit(&c->itemType)); 50 51 /* Initialize the local queue for char foot values */ 52 PetscCall(VecGetLocalSize(c->velocity, &c->queueMax)); 53 PetscCall(PetscMalloc1(c->queueMax, &c->queue)); 54 c->queueSize = 0; 55 56 /* Allocate communication structures */ 57 PetscCheck(c->numNeighbors > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %" PetscInt_FMT ". Call CharactersiticSetNeighbors() before setup.", c->numNeighbors); 58 PetscCall(PetscMalloc1(c->numNeighbors, &c->needCount)); 59 PetscCall(PetscMalloc1(c->numNeighbors, &c->localOffsets)); 60 PetscCall(PetscMalloc1(c->numNeighbors, &c->fillCount)); 61 PetscCall(PetscMalloc1(c->numNeighbors, &c->remoteOffsets)); 62 PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->request)); 63 PetscCall(PetscMalloc1(c->numNeighbors - 1, &c->status)); 64 PetscFunctionReturn(0); 65 } 66 67 PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c) { 68 Characteristic_DA *da; 69 70 PetscFunctionBegin; 71 PetscCall(PetscNew(&da)); 72 PetscCall(PetscMemzero(da, sizeof(Characteristic_DA))); 73 PetscCall(PetscLogObjectMemory((PetscObject)c, sizeof(Characteristic_DA))); 74 c->data = (void *)da; 75 76 c->ops->setup = CharacteristicSetUp_DA; 77 c->ops->destroy = CharacteristicDestroy_DA; 78 c->ops->view = CharacteristicView_DA; 79 80 da->dummy = 0; 81 PetscFunctionReturn(0); 82 } 83 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 DMBoundaryType bx, by; 90 PetscInt dim, gx, gy; 91 92 PetscFunctionBegin; 93 PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 94 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 PetscFunctionReturn(0); 104 } 105