1 #include <../src/ts/characteristic/impls/da/slda.h> /*I "petsccharacteristic.h" I*/ 2 #include <petscdmda.h> 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "CharacteristicView_DA" 6 PetscErrorCode CharacteristicView_DA(Characteristic c, PetscViewer viewer) 7 { 8 Characteristic_DA *da = (Characteristic_DA*) c->data; 9 PetscBool iascii, isstring; 10 PetscErrorCode ierr; 11 12 PetscFunctionBegin; 13 /* Pull out field names from DM */ 14 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr); 15 ierr = PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERSTRING, &isstring);CHKERRQ(ierr); 16 if (iascii) { 17 ierr = PetscViewerASCIIPrintf(viewer," DMDA: dummy=%D\n", da->dummy);CHKERRQ(ierr); 18 } else if (isstring) { 19 ierr = PetscViewerStringSPrintf(viewer,"dummy %D", da->dummy);CHKERRQ(ierr); 20 } 21 PetscFunctionReturn(0); 22 } 23 24 #undef __FUNCT__ 25 #define __FUNCT__ "CharacteristicDestroy_DA" 26 PetscErrorCode CharacteristicDestroy_DA(Characteristic c) 27 { 28 Characteristic_DA *da = (Characteristic_DA*) c->data; 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = PetscFree(da);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "CharacteristicSetUp_DA" 38 PetscErrorCode CharacteristicSetUp_DA(Characteristic c) 39 { 40 PetscMPIInt blockLen[2]; 41 MPI_Aint indices[2]; 42 MPI_Datatype oldtypes[2]; 43 PetscInt dim, numValues; 44 PetscErrorCode ierr; 45 46 ierr = DMDAGetInfo(c->velocityDA, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,0,0);CHKERRQ(ierr); 47 if (c->structured) c->numIds = dim; 48 else c->numIds = 3; 49 if (c->numFieldComp > MAX_COMPONENTS) SETERRQ2(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); 50 numValues = 4 + MAX_COMPONENTS; 51 52 /* Create new MPI datatype for communication of characteristic point structs */ 53 blockLen[0] = 1+c->numIds; indices[0] = 0; oldtypes[0] = MPIU_INT; 54 blockLen[1] = numValues; indices[1] = (1+c->numIds)*sizeof(PetscInt); oldtypes[1] = MPIU_SCALAR; 55 ierr = MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);CHKERRQ(ierr); 56 ierr = MPI_Type_commit(&c->itemType);CHKERRQ(ierr); 57 58 /* Initialize the local queue for char foot values */ 59 ierr = VecGetLocalSize(c->velocity, &c->queueMax);CHKERRQ(ierr); 60 ierr = PetscMalloc(c->queueMax * sizeof(CharacteristicPointDA2D), &c->queue);CHKERRQ(ierr); 61 c->queueSize = 0; 62 63 /* Allocate communication structures */ 64 if (c->numNeighbors <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. Call CharactersiticSetNeighbors() before setup.", c->numNeighbors); 65 ierr = PetscMalloc(c->numNeighbors * sizeof(PetscInt), &c->needCount);CHKERRQ(ierr); 66 ierr = PetscMalloc(c->numNeighbors * sizeof(PetscInt), &c->localOffsets);CHKERRQ(ierr); 67 ierr = PetscMalloc(c->numNeighbors * sizeof(PetscInt), &c->fillCount);CHKERRQ(ierr); 68 ierr = PetscMalloc(c->numNeighbors * sizeof(PetscInt), &c->remoteOffsets);CHKERRQ(ierr); 69 ierr = PetscMalloc((c->numNeighbors-1) * sizeof(MPI_Request), &c->request);CHKERRQ(ierr); 70 ierr = PetscMalloc((c->numNeighbors-1) * sizeof(MPI_Status), &c->status);CHKERRQ(ierr); 71 PetscFunctionReturn(0); 72 } 73 74 EXTERN_C_BEGIN 75 #undef __FUNCT__ 76 #define __FUNCT__ "CharacteristicCreate_DA" 77 PetscErrorCode CharacteristicCreate_DA(Characteristic c) 78 { 79 Characteristic_DA *da; 80 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 ierr = PetscNew(Characteristic_DA, &da);CHKERRQ(ierr); 84 ierr = PetscMemzero(da, sizeof(Characteristic_DA));CHKERRQ(ierr); 85 ierr = PetscLogObjectMemory(c, sizeof(Characteristic_DA));CHKERRQ(ierr); 86 c->data = (void*) da; 87 88 c->ops->setup = CharacteristicSetUp_DA; 89 c->ops->destroy = CharacteristicDestroy_DA; 90 c->ops->view = CharacteristicView_DA; 91 92 da->dummy = 0; 93 PetscFunctionReturn(0); 94 } 95 EXTERN_C_END 96 97 #undef __FUNCT__ 98 #define __FUNCT__ "DMDAMapCoordsToPeriodicDomain" 99 /* ----------------------------------------------------------------------------- 100 Checks for periodicity of a DM and Maps points outside of a domain back onto the domain 101 using appropriate periodicity. At the moment assumes only a 2-D DMDA. 102 ----------------------------------------------------------------------------------------*/ 103 PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y) 104 { 105 DMDABoundaryType bx, by; 106 PetscInt dim, gx, gy; 107 PetscErrorCode ierr; 108 109 PetscFunctionBegin; 110 ierr = DMDAGetInfo(da, &dim, &gx, &gy, 0, 0, 0, 0, 0, 0, &bx, &by, 0, 0); 111 112 if (bx == DMDA_BOUNDARY_PERIODIC) { 113 while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx; 114 while (*x < 0.0) *x += (PetscScalar)gx; 115 } 116 if (by == DMDA_BOUNDARY_PERIODIC) { 117 while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy; 118 while (*y < 0.0) *y += (PetscScalar)gy; 119 } 120 PetscFunctionReturn(ierr); 121 } 122