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