1 #include <../src/ts/characteristic/impls/da/slda.h> /*I "petsccharacteristic.h" I*/ 2 #include <petscdmda.h> 3 #include <petscviewer.h> 4 5 static PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer) 6 { 7 Characteristic_DA *da = (Characteristic_DA *)c->data; 8 PetscBool isascii, isstring; 9 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 } 21 22 static PetscErrorCode CharacteristicDestroy_DA(Characteristic c) 23 { 24 Characteristic_DA *da = (Characteristic_DA *)c->data; 25 26 PetscFunctionBegin; 27 PetscCall(PetscFree(da)); 28 PetscFunctionReturn(PETSC_SUCCESS); 29 } 30 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; 37 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; 43 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)); 53 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; 58 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 } 69 70 PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c) 71 { 72 Characteristic_DA *da; 73 74 PetscFunctionBegin; 75 PetscCall(PetscNew(&da)); 76 PetscCall(PetscMemzero(da, sizeof(Characteristic_DA))); 77 c->data = (void *)da; 78 79 c->ops->setup = CharacteristicSetUp_DA; 80 c->ops->destroy = CharacteristicDestroy_DA; 81 c->ops->view = CharacteristicView_DA; 82 83 da->dummy = 0; 84 PetscFunctionReturn(PETSC_SUCCESS); 85 } 86 87 /* ----------------------------------------------------------------------------- 88 Checks for periodicity of a DM and Maps points outside of a domain back onto the domain 89 using appropriate periodicity. At the moment assumes only a 2-D DMDA. 90 ----------------------------------------------------------------------------------------*/ 91 PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y) 92 { 93 DMBoundaryType bx, by; 94 PetscInt dim, gx, gy; 95 96 PetscFunctionBegin; 97 PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 98 99 if (bx == DM_BOUNDARY_PERIODIC) { 100 while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx; 101 while (*x < 0.0) *x += (PetscScalar)gx; 102 } 103 if (by == DM_BOUNDARY_PERIODIC) { 104 while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy; 105 while (*y < 0.0) *y += (PetscScalar)gy; 106 } 107 PetscFunctionReturn(PETSC_SUCCESS); 108 } 109