Actual source code: drivers.cxx
1: #include "contexts.cxx"
2: #include "sparse.cxx"
3: #include "init.cxx"
4: #include <adolc/drivers/drivers.h>
5: #include <adolc/interfaces.h>
7: /*
8: REQUIRES configuration of PETSc with option --download-adolc.
10: For documentation on ADOL-C, see
11: $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
12: */
14: /* --------------------------------------------------------------------------------
15: Drivers for RHSJacobian and IJacobian
16: ----------------------------------------------------------------------------- */
18: /*
19: Compute Jacobian for explicit TS in compressed format and recover from this, using
20: precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
21: assembled (not recommended for non-toy problems!).
23: Input parameters:
24: tag - tape identifier
25: u_vec - vector at which to evaluate Jacobian
26: ctx - ADOL-C context, as defined above
28: Output parameter:
29: A - Mat object corresponding to Jacobian
30: */
31: PetscErrorCode PetscAdolcComputeRHSJacobian(PetscInt tag, Mat A, const PetscScalar *u_vec, void *ctx)
32: {
33: AdolcCtx *adctx = (AdolcCtx *)ctx;
34: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
35: PetscScalar **J;
37: PetscFunctionBegin;
38: PetscCall(AdolcMalloc2(m, p, &J));
39: if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
40: else jacobian(tag, m, n, u_vec, J);
41: if (adctx->sparse) {
42: PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
43: } else {
44: for (i = 0; i < m; i++) {
45: for (j = 0; j < n; j++) {
46: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
47: }
48: }
49: }
50: PetscCall(AdolcFree2(J));
51: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
52: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: /*
57: Compute Jacobian for explicit TS in compressed format and recover from this, using
58: precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
59: assembled (not recommended for non-toy problems!).
61: Input parameters:
62: tag - tape identifier
63: u_vec - vector at which to evaluate Jacobian
64: ctx - ADOL-C context, as defined above
66: Output parameter:
67: A - Mat object corresponding to Jacobian
68: */
69: PetscErrorCode PetscAdolcComputeRHSJacobianLocal(PetscInt tag, Mat A, const PetscScalar *u_vec, void *ctx)
70: {
71: AdolcCtx *adctx = (AdolcCtx *)ctx;
72: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
73: PetscScalar **J;
75: PetscFunctionBegin;
76: PetscCall(AdolcMalloc2(m, p, &J));
77: if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
78: else jacobian(tag, m, n, u_vec, J);
79: if (adctx->sparse) {
80: PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
81: } else {
82: for (i = 0; i < m; i++) {
83: for (j = 0; j < n; j++) {
84: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
85: }
86: }
87: }
88: PetscCall(AdolcFree2(J));
89: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
90: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /*
95: Compute Jacobian for implicit TS in compressed format and recover from this, using
96: precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
97: assembled (not recommended for non-toy problems!).
99: Input parameters:
100: tag1 - tape identifier for df/dx part
101: tag2 - tape identifier for df/d(xdot) part
102: u_vec - vector at which to evaluate Jacobian
103: ctx - ADOL-C context, as defined above
105: Output parameter:
106: A - Mat object corresponding to Jacobian
107: */
108: PetscErrorCode PetscAdolcComputeIJacobian(PetscInt tag1, PetscInt tag2, Mat A, const PetscScalar *u_vec, PetscReal a, void *ctx)
109: {
110: AdolcCtx *adctx = (AdolcCtx *)ctx;
111: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
112: PetscScalar **J;
114: PetscFunctionBegin;
115: PetscCall(AdolcMalloc2(m, p, &J));
117: /* dF/dx part */
118: if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
119: else jacobian(tag1, m, n, u_vec, J);
120: PetscCall(MatZeroEntries(A));
121: if (adctx->sparse) {
122: PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
123: } else {
124: for (i = 0; i < m; i++) {
125: for (j = 0; j < n; j++) {
126: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
127: }
128: }
129: }
130: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
131: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
133: /* a * dF/d(xdot) part */
134: if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
135: else jacobian(tag2, m, n, u_vec, J);
136: if (adctx->sparse) {
137: PetscCall(RecoverJacobian(A, ADD_VALUES, m, p, adctx->Rec, J, &a));
138: } else {
139: for (i = 0; i < m; i++) {
140: for (j = 0; j < n; j++) {
141: if (fabs(J[i][j]) > 1.e-16) {
142: J[i][j] *= a;
143: PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], ADD_VALUES));
144: }
145: }
146: }
147: }
148: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
149: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
150: PetscCall(AdolcFree2(J));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: /*
155: Compute Jacobian for implicit TS in the special case where it is
156: known that the mass matrix is simply the identity. i.e. We have
157: a problem of the form
158: du/dt = F(u).
160: Input parameters:
161: tag - tape identifier for df/dx part
162: u_vec - vector at which to evaluate Jacobian
163: ctx - ADOL-C context, as defined above
165: Output parameter:
166: A - Mat object corresponding to Jacobian
167: */
168: PetscErrorCode PetscAdolcComputeIJacobianIDMass(PetscInt tag, Mat A, PetscScalar *u_vec, PetscReal a, void *ctx)
169: {
170: AdolcCtx *adctx = (AdolcCtx *)ctx;
171: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
172: PetscScalar **J;
174: PetscFunctionBegin;
175: PetscCall(AdolcMalloc2(m, p, &J));
177: /* dF/dx part */
178: if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
179: else jacobian(tag, m, n, u_vec, J);
180: PetscCall(MatZeroEntries(A));
181: if (adctx->sparse) {
182: PetscCall(RecoverJacobian(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
183: } else {
184: for (i = 0; i < m; i++) {
185: for (j = 0; j < n; j++) {
186: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
187: }
188: }
189: }
190: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
191: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
192: PetscCall(AdolcFree2(J));
194: /* a * dF/d(xdot) part */
195: PetscCall(MatShift(A, a));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*
200: Compute local portion of Jacobian for implicit TS in compressed format and recover from this, using
201: precomputed seed and recovery matrices. If sparse mode is not used, full Jacobian is
202: assembled (not recommended for non-toy problems!).
204: Input parameters:
205: tag1 - tape identifier for df/dx part
206: tag2 - tape identifier for df/d(xdot) part
207: u_vec - vector at which to evaluate Jacobian
208: ctx - ADOL-C context, as defined above
210: Output parameter:
211: A - Mat object corresponding to Jacobian
212: */
213: PetscErrorCode PetscAdolcComputeIJacobianLocal(PetscInt tag1, PetscInt tag2, Mat A, PetscScalar *u_vec, PetscReal a, void *ctx)
214: {
215: AdolcCtx *adctx = (AdolcCtx *)ctx;
216: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
217: PetscScalar **J;
219: PetscFunctionBegin;
220: PetscCall(AdolcMalloc2(m, p, &J));
222: /* dF/dx part */
223: if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
224: else jacobian(tag1, m, n, u_vec, J);
225: if (adctx->sparse) {
226: PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
227: } else {
228: for (i = 0; i < m; i++) {
229: for (j = 0; j < n; j++) {
230: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
231: }
232: }
233: }
234: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
235: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
237: /* a * dF/d(xdot) part */
238: if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
239: else jacobian(tag2, m, n, u_vec, J);
240: if (adctx->sparse) {
241: PetscCall(RecoverJacobianLocal(A, ADD_VALUES, m, p, adctx->Rec, J, &a));
242: } else {
243: for (i = 0; i < m; i++) {
244: for (j = 0; j < n; j++) {
245: if (fabs(J[i][j]) > 1.e-16) {
246: J[i][j] *= a;
247: PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], ADD_VALUES));
248: }
249: }
250: }
251: }
252: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
253: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
254: PetscCall(AdolcFree2(J));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*
259: Compute local portion of Jacobian for implicit TS in the special case where it is
260: known that the mass matrix is simply the identity. i.e. We have
261: a problem of the form
262: du/dt = F(u).
264: Input parameters:
265: tag - tape identifier for df/dx part
266: u_vec - vector at which to evaluate Jacobian
267: ctx - ADOL-C context, as defined above
269: Output parameter:
270: A - Mat object corresponding to Jacobian
271: */
272: PetscErrorCode PetscAdolcComputeIJacobianLocalIDMass(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscReal a, void *ctx)
273: {
274: AdolcCtx *adctx = (AdolcCtx *)ctx;
275: PetscInt i, j, m = adctx->m, n = adctx->n, p = adctx->p;
276: PetscScalar **J;
278: PetscFunctionBegin;
279: PetscCall(AdolcMalloc2(m, p, &J));
281: /* dF/dx part */
282: if (adctx->Seed) fov_forward(tag, m, n, p, u_vec, adctx->Seed, NULL, J);
283: else jacobian(tag, m, n, u_vec, J);
284: if (adctx->sparse) {
285: PetscCall(RecoverJacobianLocal(A, INSERT_VALUES, m, p, adctx->Rec, J, NULL));
286: } else {
287: for (i = 0; i < m; i++) {
288: for (j = 0; j < n; j++) {
289: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
290: }
291: }
292: }
293: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
294: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
295: PetscCall(AdolcFree2(J));
297: /* a * dF/d(xdot) part */
298: PetscCall(MatShift(A, a));
299: PetscFunctionReturn(PETSC_SUCCESS);
300: }
302: /* --------------------------------------------------------------------------------
303: Drivers for Jacobian w.r.t. a parameter
304: ----------------------------------------------------------------------------- */
306: /*
307: Compute Jacobian w.r.t a parameter for explicit TS.
309: Input parameters:
310: tag - tape identifier
311: u_vec - vector at which to evaluate Jacobian
312: params - the parameters w.r.t. which we differentiate
313: ctx - ADOL-C context, as defined above
315: Output parameter:
316: A - Mat object corresponding to Jacobian
317: */
318: PetscErrorCode PetscAdolcComputeRHSJacobianP(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscScalar *params, void *ctx)
319: {
320: AdolcCtx *adctx = (AdolcCtx *)ctx;
321: PetscInt i, j = 0, m = adctx->m, n = adctx->n, p = adctx->num_params;
322: PetscScalar **J, *concat, **S;
324: PetscFunctionBegin;
325: /* Allocate memory and concatenate independent variable values with parameter */
326: PetscCall(AdolcMalloc2(m, p, &J));
327: PetscCall(PetscMalloc1(n + p, &concat));
328: PetscCall(AdolcMalloc2(n + p, p, &S));
329: PetscCall(Subidentity(p, n, S));
330: for (i = 0; i < n; i++) concat[i] = u_vec[i];
331: for (i = 0; i < p; i++) concat[n + i] = params[i];
333: /* Propagate the appropriate seed matrix through the forward mode of AD */
334: fov_forward(tag, m, n + p, p, concat, S, NULL, J);
335: PetscCall(AdolcFree2(S));
336: PetscCall(PetscFree(concat));
338: /* Set matrix values */
339: for (i = 0; i < m; i++) {
340: for (j = 0; j < p; j++) {
341: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
342: }
343: }
344: PetscCall(AdolcFree2(J));
345: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
346: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: /*
351: Compute local portion of Jacobian w.r.t a parameter for explicit TS.
353: Input parameters:
354: tag - tape identifier
355: u_vec - vector at which to evaluate Jacobian
356: params - the parameters w.r.t. which we differentiate
357: ctx - ADOL-C context, as defined above
359: Output parameter:
360: A - Mat object corresponding to Jacobian
361: */
362: PetscErrorCode PetscAdolcComputeRHSJacobianPLocal(PetscInt tag, Mat A, const PetscScalar *u_vec, PetscScalar *params, void *ctx)
363: {
364: AdolcCtx *adctx = (AdolcCtx *)ctx;
365: PetscInt i, j = 0, m = adctx->m, n = adctx->n, p = adctx->num_params;
366: PetscScalar **J, *concat, **S;
368: PetscFunctionBegin;
369: /* Allocate memory and concatenate independent variable values with parameter */
370: PetscCall(AdolcMalloc2(m, p, &J));
371: PetscCall(PetscMalloc1(n + p, &concat));
372: PetscCall(AdolcMalloc2(n + p, p, &S));
373: PetscCall(Subidentity(p, n, S));
374: for (i = 0; i < n; i++) concat[i] = u_vec[i];
375: for (i = 0; i < p; i++) concat[n + i] = params[i];
377: /* Propagate the appropriate seed matrix through the forward mode of AD */
378: fov_forward(tag, m, n + p, p, concat, S, NULL, J);
379: PetscCall(AdolcFree2(S));
380: PetscCall(PetscFree(concat));
382: /* Set matrix values */
383: for (i = 0; i < m; i++) {
384: for (j = 0; j < p; j++) {
385: if (fabs(J[i][j]) > 1.e-16) PetscCall(MatSetValuesLocal(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
386: }
387: }
388: PetscCall(AdolcFree2(J));
389: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
390: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /* --------------------------------------------------------------------------------
395: Drivers for Jacobian diagonal
396: ----------------------------------------------------------------------------- */
398: /*
399: Compute local portion of Jacobian diagonal for implicit TS in compressed format and recover
400: from this, using precomputed seed matrix and recovery vector.
402: Input parameters:
403: tag1 - tape identifier for df/dx part
404: tag2 - tape identifier for df/d(xdot) part
405: u_vec - vector at which to evaluate Jacobian
406: ctx - ADOL-C context, as defined above
408: Output parameter:
409: diag - Vec object corresponding to Jacobian diagonal
410: */
411: PetscErrorCode PetscAdolcComputeIJacobianAndDiagonalLocal(PetscInt tag1, PetscInt tag2, Vec diag, PetscScalar *u_vec, PetscReal a, void *ctx)
412: {
413: AdolcCtx *adctx = (AdolcCtx *)ctx;
414: PetscInt i, m = adctx->m, n = adctx->n, p = adctx->p;
415: PetscScalar **J;
417: PetscFunctionBegin;
418: PetscCall(AdolcMalloc2(m, p, &J));
420: /* dF/dx part */
421: if (adctx->Seed) fov_forward(tag1, m, n, p, u_vec, adctx->Seed, NULL, J);
422: else jacobian(tag1, m, n, u_vec, J);
423: if (adctx->sparse) {
424: PetscCall(RecoverDiagonalLocal(diag, INSERT_VALUES, m, adctx->rec, J, NULL));
425: } else {
426: for (i = 0; i < m; i++) {
427: if (fabs(J[i][i]) > 1.e-16) PetscCall(VecSetValuesLocal(diag, 1, &i, &J[i][i], INSERT_VALUES));
428: }
429: }
430: PetscCall(VecAssemblyBegin(diag));
431: PetscCall(VecAssemblyEnd(diag));
433: /* a * dF/d(xdot) part */
434: if (adctx->Seed) fov_forward(tag2, m, n, p, u_vec, adctx->Seed, NULL, J);
435: else jacobian(tag2, m, n, u_vec, J);
436: if (adctx->sparse) {
437: PetscCall(RecoverDiagonalLocal(diag, ADD_VALUES, m, adctx->rec, J, NULL));
438: } else {
439: for (i = 0; i < m; i++) {
440: if (fabs(J[i][i]) > 1.e-16) {
441: J[i][i] *= a;
442: PetscCall(VecSetValuesLocal(diag, 1, &i, &J[i][i], ADD_VALUES));
443: }
444: }
445: }
446: PetscCall(VecAssemblyBegin(diag));
447: PetscCall(VecAssemblyEnd(diag));
448: PetscCall(AdolcFree2(J));
449: PetscFunctionReturn(PETSC_SUCCESS);
450: }