1c4762a1bSJed Brown #include <petscdmda.h> 2c4762a1bSJed Brown #include <adolc/adalloc.h> 3c4762a1bSJed Brown 4c4762a1bSJed Brown /* 5c4762a1bSJed Brown REQUIRES configuration of PETSc with option --download-adolc. 6c4762a1bSJed Brown 7c4762a1bSJed Brown For documentation on ADOL-C, see 8c4762a1bSJed Brown $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 9c4762a1bSJed Brown */ 10c4762a1bSJed Brown 11c4762a1bSJed Brown /* 12c4762a1bSJed Brown Wrapper function for allocating contiguous memory in a 2d array 13c4762a1bSJed Brown 14c4762a1bSJed Brown Input parameters: 15c4762a1bSJed Brown m,n - number of rows and columns of array, respectively 16c4762a1bSJed Brown 17c4762a1bSJed Brown Output parameter: 18c4762a1bSJed Brown A - pointer to array for which memory is allocated 19c4762a1bSJed Brown 20c4762a1bSJed Brown Note: Only arrays of doubles are currently accounted for in ADOL-C's myalloc2 function. 21c4762a1bSJed Brown */ 229371c9d4SSatish Balay template <class T> 23d71ae5a4SJacob Faibussowitsch PetscErrorCode AdolcMalloc2(PetscInt m, PetscInt n, T **A[]) 24d71ae5a4SJacob Faibussowitsch { 25c4762a1bSJed Brown PetscFunctionBegin; 26c4762a1bSJed Brown *A = myalloc2(m, n); 27*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 28c4762a1bSJed Brown } 29c4762a1bSJed Brown 30c4762a1bSJed Brown /* 31c4762a1bSJed Brown Wrapper function for freeing memory allocated with AdolcMalloc2 32c4762a1bSJed Brown 33c4762a1bSJed Brown Input parameter: 34c4762a1bSJed Brown A - array to free memory of 35c4762a1bSJed Brown 36c4762a1bSJed Brown Note: Only arrays of doubles are currently accounted for in ADOL-C's myfree2 function. 37c4762a1bSJed Brown */ 389371c9d4SSatish Balay template <class T> 39d71ae5a4SJacob Faibussowitsch PetscErrorCode AdolcFree2(T **A) 40d71ae5a4SJacob Faibussowitsch { 41c4762a1bSJed Brown PetscFunctionBegin; 42c4762a1bSJed Brown myfree2(A); 43*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 44c4762a1bSJed Brown } 45c4762a1bSJed Brown 46c4762a1bSJed Brown /* 47c4762a1bSJed Brown Shift indices in an array of type T to endow it with ghost points. 48c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 49c4762a1bSJed Brown 50c4762a1bSJed Brown Input parameters: 51c4762a1bSJed Brown da - distributed array upon which variables are defined 52c4762a1bSJed Brown cgs - contiguously allocated 1-array with as many entries as there are 53c4762a1bSJed Brown interior and ghost points, in total 54c4762a1bSJed Brown 55c4762a1bSJed Brown Output parameter: 56c4762a1bSJed Brown array - contiguously allocated array of the appropriate dimension with 57c4762a1bSJed Brown ghost points, pointing to the 1-array 58c4762a1bSJed Brown */ 599371c9d4SSatish Balay template <class T> 60d71ae5a4SJacob Faibussowitsch PetscErrorCode GiveGhostPoints(DM da, T *cgs, void *array) 61d71ae5a4SJacob Faibussowitsch { 62c4762a1bSJed Brown PetscInt dim; 63c4762a1bSJed Brown 64c4762a1bSJed Brown PetscFunctionBegin; 659566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, &dim, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)); 66c4762a1bSJed Brown if (dim == 1) { 679566063dSJacob Faibussowitsch PetscCall(GiveGhostPoints1d(da, (T **)array)); 68c4762a1bSJed Brown } else if (dim == 2) { 699566063dSJacob Faibussowitsch PetscCall(GiveGhostPoints2d(da, cgs, (T ***)array)); 7008401ef6SPierre Jolivet } else PetscCheck(dim != 3, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "GiveGhostPoints3d not yet implemented"); // TODO 71*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 72c4762a1bSJed Brown } 73c4762a1bSJed Brown 74c4762a1bSJed Brown /* 75c4762a1bSJed Brown Shift indices in a 1-array of type T to endow it with ghost points. 76c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 77c4762a1bSJed Brown 78c4762a1bSJed Brown Input parameters: 79c4762a1bSJed Brown da - distributed array upon which variables are defined 80c4762a1bSJed Brown 81c4762a1bSJed Brown Output parameter: 82c4762a1bSJed Brown a1d - contiguously allocated 1-array 83c4762a1bSJed Brown */ 849371c9d4SSatish Balay template <class T> 85d71ae5a4SJacob Faibussowitsch PetscErrorCode GiveGhostPoints1d(DM da, T *a1d[]) 86d71ae5a4SJacob Faibussowitsch { 87c4762a1bSJed Brown PetscInt gxs; 88c4762a1bSJed Brown 89c4762a1bSJed Brown PetscFunctionBegin; 909566063dSJacob Faibussowitsch PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, NULL, NULL, NULL)); 91c4762a1bSJed Brown *a1d -= gxs; 92*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 93c4762a1bSJed Brown } 94c4762a1bSJed Brown 95c4762a1bSJed Brown /* 96c4762a1bSJed Brown Shift indices in a 2-array of type T to endow it with ghost points. 97c4762a1bSJed Brown (e.g. This works for arrays of adoubles or arrays (of structs) thereof.) 98c4762a1bSJed Brown 99c4762a1bSJed Brown Input parameters: 100c4762a1bSJed Brown da - distributed array upon which variables are defined 101c4762a1bSJed Brown cgs - contiguously allocated 1-array with as many entries as there are 102c4762a1bSJed Brown interior and ghost points, in total 103c4762a1bSJed Brown 104c4762a1bSJed Brown Output parameter: 105c4762a1bSJed Brown a2d - contiguously allocated 2-array with ghost points, pointing to the 106c4762a1bSJed Brown 1-array 107c4762a1bSJed Brown */ 1089371c9d4SSatish Balay template <class T> 109d71ae5a4SJacob Faibussowitsch PetscErrorCode GiveGhostPoints2d(DM da, T *cgs, T **a2d[]) 110d71ae5a4SJacob Faibussowitsch { 1115f80ce2aSJacob Faibussowitsch PetscInt gxs, gys, gxm, gym; 112c4762a1bSJed Brown 113c4762a1bSJed Brown PetscFunctionBegin; 1149566063dSJacob Faibussowitsch PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gxm, &gym, NULL)); 1155f80ce2aSJacob Faibussowitsch for (PetscInt j = 0; j < gym; j++) (*a2d)[j] = cgs + j * gxm - gxs; 116c4762a1bSJed Brown *a2d -= gys; 117*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 118c4762a1bSJed Brown } 119c4762a1bSJed Brown 120c4762a1bSJed Brown /* 121c4762a1bSJed Brown Create a rectangular sub-identity of the m x m identity matrix, as an array. 122c4762a1bSJed Brown 123c4762a1bSJed Brown Input parameters: 124c4762a1bSJed Brown n - number of (adjacent) rows to take in slice 125c4762a1bSJed Brown s - starting row index 126c4762a1bSJed Brown 127c4762a1bSJed Brown Output parameter: 128c4762a1bSJed Brown S - resulting n x m submatrix 129c4762a1bSJed Brown */ 1309371c9d4SSatish Balay template <class T> 131d71ae5a4SJacob Faibussowitsch PetscErrorCode Subidentity(PetscInt n, PetscInt s, T **S) 132d71ae5a4SJacob Faibussowitsch { 133c4762a1bSJed Brown PetscInt i; 134c4762a1bSJed Brown 135c4762a1bSJed Brown PetscFunctionBegin; 136ad540459SPierre Jolivet for (i = 0; i < n; i++) S[i][i + s] = 1.; 137*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 138c4762a1bSJed Brown } 139c4762a1bSJed Brown 140c4762a1bSJed Brown /* 141c4762a1bSJed Brown Create an identity matrix, as an array. 142c4762a1bSJed Brown 143c4762a1bSJed Brown Input parameter: 144c4762a1bSJed Brown n - number of rows/columns 145c4762a1bSJed Brown I - n x n array with memory pre-allocated 146c4762a1bSJed Brown */ 1479371c9d4SSatish Balay template <class T> 148d71ae5a4SJacob Faibussowitsch PetscErrorCode Identity(PetscInt n, T **I) 149d71ae5a4SJacob Faibussowitsch { 150c4762a1bSJed Brown PetscFunctionBegin; 1519566063dSJacob Faibussowitsch PetscCall(Subidentity(n, 0, I)); 152*3ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 153c4762a1bSJed Brown } 154