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