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