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