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 { 7 Characteristic_DA *da = (Characteristic_DA*) c->data; 8 PetscBool iascii, isstring; 9 10 PetscFunctionBegin; 11 /* Pull out field names from DM */ 12 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii)); 13 PetscCall(PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERSTRING, &isstring)); 14 if (iascii) { 15 PetscCall(PetscViewerASCIIPrintf(viewer," DMDA: dummy=%D\n", da->dummy)); 16 } else if (isstring) { 17 PetscCall(PetscViewerStringSPrintf(viewer,"dummy %D", da->dummy)); 18 } 19 PetscFunctionReturn(0); 20 } 21 22 PetscErrorCode CharacteristicDestroy_DA(Characteristic c) 23 { 24 Characteristic_DA *da = (Characteristic_DA*) c->data; 25 26 PetscFunctionBegin; 27 PetscCall(PetscFree(da)); 28 PetscFunctionReturn(0); 29 } 30 31 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) c->numIds = dim; 40 else c->numIds = 3; 41 PetscCheckFalse(c->numFieldComp > MAX_COMPONENTS,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "The maximum number of fields allowed is %d, you have %d. 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; indices[0] = 0; oldtypes[0] = MPIU_INT; 46 blockLen[1] = numValues; indices[1] = (1+c->numIds)*sizeof(PetscInt); oldtypes[1] = MPIU_SCALAR; 47 PetscCallMPI(MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType)); 48 PetscCallMPI(MPI_Type_commit(&c->itemType)); 49 50 /* Initialize the local queue for char foot values */ 51 PetscCall(VecGetLocalSize(c->velocity, &c->queueMax)); 52 PetscCall(PetscMalloc1(c->queueMax, &c->queue)); 53 c->queueSize = 0; 54 55 /* Allocate communication structures */ 56 PetscCheckFalse(c->numNeighbors <= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. Call CharactersiticSetNeighbors() before setup.", c->numNeighbors); 57 PetscCall(PetscMalloc1(c->numNeighbors, &c->needCount)); 58 PetscCall(PetscMalloc1(c->numNeighbors, &c->localOffsets)); 59 PetscCall(PetscMalloc1(c->numNeighbors, &c->fillCount)); 60 PetscCall(PetscMalloc1(c->numNeighbors, &c->remoteOffsets)); 61 PetscCall(PetscMalloc1(c->numNeighbors-1, &c->request)); 62 PetscCall(PetscMalloc1(c->numNeighbors-1, &c->status)); 63 PetscFunctionReturn(0); 64 } 65 66 PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c) 67 { 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 { 90 DMBoundaryType bx, by; 91 PetscInt dim, gx, gy; 92 93 PetscFunctionBegin; 94 PetscCall(DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL)); 95 96 if (bx == DM_BOUNDARY_PERIODIC) { 97 while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx; 98 while (*x < 0.0) *x += (PetscScalar)gx; 99 } 100 if (by == DM_BOUNDARY_PERIODIC) { 101 while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy; 102 while (*y < 0.0) *y += (PetscScalar)gy; 103 } 104 PetscFunctionReturn(0); 105 } 106