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 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 } 20 PetscFunctionReturn(0); 21 } 22 23 PetscErrorCode CharacteristicDestroy_DA(Characteristic c) 24 { 25 Characteristic_DA *da = (Characteristic_DA*) c->data; 26 PetscErrorCode ierr; 27 28 PetscFunctionBegin; 29 ierr = PetscFree(da);CHKERRQ(ierr); 30 PetscFunctionReturn(0); 31 } 32 33 PetscErrorCode CharacteristicSetUp_DA(Characteristic c) 34 { 35 PetscMPIInt blockLen[2]; 36 MPI_Aint indices[2]; 37 MPI_Datatype oldtypes[2]; 38 PetscInt dim, numValues; 39 PetscErrorCode ierr; 40 41 ierr = DMDAGetInfo(c->velocityDA, &dim, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);CHKERRQ(ierr); 42 if (c->structured) c->numIds = dim; 43 else c->numIds = 3; 44 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); 45 numValues = 4 + MAX_COMPONENTS; 46 47 /* Create new MPI datatype for communication of characteristic point structs */ 48 blockLen[0] = 1+c->numIds; indices[0] = 0; oldtypes[0] = MPIU_INT; 49 blockLen[1] = numValues; indices[1] = (1+c->numIds)*sizeof(PetscInt); oldtypes[1] = MPIU_SCALAR; 50 ierr = MPI_Type_create_struct(2, blockLen, indices, oldtypes, &c->itemType);CHKERRMPI(ierr); 51 ierr = MPI_Type_commit(&c->itemType);CHKERRMPI(ierr); 52 53 /* Initialize the local queue for char foot values */ 54 ierr = VecGetLocalSize(c->velocity, &c->queueMax);CHKERRQ(ierr); 55 ierr = PetscMalloc1(c->queueMax, &c->queue);CHKERRQ(ierr); 56 c->queueSize = 0; 57 58 /* Allocate communication structures */ 59 if (c->numNeighbors <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE, "Invalid number of neighbors %d. Call CharactersiticSetNeighbors() before setup.", c->numNeighbors); 60 ierr = PetscMalloc1(c->numNeighbors, &c->needCount);CHKERRQ(ierr); 61 ierr = PetscMalloc1(c->numNeighbors, &c->localOffsets);CHKERRQ(ierr); 62 ierr = PetscMalloc1(c->numNeighbors, &c->fillCount);CHKERRQ(ierr); 63 ierr = PetscMalloc1(c->numNeighbors, &c->remoteOffsets);CHKERRQ(ierr); 64 ierr = PetscMalloc1(c->numNeighbors-1, &c->request);CHKERRQ(ierr); 65 ierr = PetscMalloc1(c->numNeighbors-1, &c->status);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 PETSC_EXTERN PetscErrorCode CharacteristicCreate_DA(Characteristic c) 70 { 71 Characteristic_DA *da; 72 PetscErrorCode ierr; 73 74 PetscFunctionBegin; 75 ierr = PetscNew(&da);CHKERRQ(ierr); 76 ierr = PetscMemzero(da, sizeof(Characteristic_DA));CHKERRQ(ierr); 77 ierr = PetscLogObjectMemory((PetscObject)c, sizeof(Characteristic_DA));CHKERRQ(ierr); 78 c->data = (void*) da; 79 80 c->ops->setup = CharacteristicSetUp_DA; 81 c->ops->destroy = CharacteristicDestroy_DA; 82 c->ops->view = CharacteristicView_DA; 83 84 da->dummy = 0; 85 PetscFunctionReturn(0); 86 } 87 88 /* ----------------------------------------------------------------------------- 89 Checks for periodicity of a DM and Maps points outside of a domain back onto the domain 90 using appropriate periodicity. At the moment assumes only a 2-D DMDA. 91 ----------------------------------------------------------------------------------------*/ 92 PetscErrorCode DMDAMapCoordsToPeriodicDomain(DM da, PetscScalar *x, PetscScalar *y) 93 { 94 DMBoundaryType bx, by; 95 PetscInt dim, gx, gy; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 ierr = DMDAGetInfo(da, &dim, &gx, &gy, NULL, NULL, NULL, NULL, NULL, NULL, &bx, &by, NULL, NULL);CHKERRQ(ierr); 100 101 if (bx == DM_BOUNDARY_PERIODIC) { 102 while (*x >= (PetscScalar)gx) *x -= (PetscScalar)gx; 103 while (*x < 0.0) *x += (PetscScalar)gx; 104 } 105 if (by == DM_BOUNDARY_PERIODIC) { 106 while (*y >= (PetscScalar)gy) *y -= (PetscScalar)gy; 107 while (*y < 0.0) *y += (PetscScalar)gy; 108 } 109 PetscFunctionReturn(0); 110 } 111