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>
AdolcMalloc2(PetscInt m,PetscInt n,T ** A[])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>
AdolcFree2(T ** A)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>
GiveGhostPoints(DM da,T * cgs,void * array)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>
GiveGhostPoints1d(DM da,T * a1d[])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>
GiveGhostPoints2d(DM da,T * cgs,T ** a2d[])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>
Subidentity(PetscInt n,PetscInt s,T ** S)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>
Identity(PetscInt n,T ** I)148 PetscErrorCode Identity(PetscInt n, T **I)
149 {
150 PetscFunctionBegin;
151 PetscCall(Subidentity(n, 0, I));
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154