1 static char help[] = "Demonstrates tapeless automatic Jacobian generation using ADOL-C for an adjoint sensitivity analysis of the van der Pol equation.\n\
2 Input parameters include:\n\
3 -mu : stiffness parameter\n\n";
4
5 /*
6 REQUIRES configuration of PETSc with option --download-adolc.
7
8 For documentation on ADOL-C, see
9 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf
10 */
11 /* ------------------------------------------------------------------------
12 See ex16adj for a description of the problem being solved.
13 ------------------------------------------------------------------------- */
14
15 #include <petscts.h>
16 #include <petscmat.h>
17
18 #define ADOLC_TAPELESS
19 #define NUMBER_DIRECTIONS 3
20 #include "adolc-utils/drivers.cxx"
21 #include <adolc/adtl.h>
22 using namespace adtl;
23
24 typedef struct _n_User *User;
25 struct _n_User {
26 PetscReal mu;
27 PetscReal next_output;
28 PetscReal tprev;
29
30 /* Automatic differentiation support */
31 AdolcCtx *adctx;
32 Vec F;
33 };
34
35 /*
36 Residual evaluation templated, so as to allow for PetscScalar or adouble
37 arguments.
38 */
39 template <class T>
EvaluateResidual(const T * x,T mu,T * f)40 PetscErrorCode EvaluateResidual(const T *x, T mu, T *f)
41 {
42 PetscFunctionBegin;
43 f[0] = x[1];
44 f[1] = mu * (1. - x[0] * x[0]) * x[1] - x[0];
45 PetscFunctionReturn(PETSC_SUCCESS);
46 }
47
48 /*
49 'Passive' RHS function, used in residual evaluations during the time integration.
50 */
RHSFunctionPassive(TS ts,PetscReal t,Vec X,Vec F,PetscCtx ctx)51 static PetscErrorCode RHSFunctionPassive(TS ts, PetscReal t, Vec X, Vec F, PetscCtx ctx)
52 {
53 User user = (User)ctx;
54 PetscScalar *f;
55 const PetscScalar *x;
56
57 PetscFunctionBeginUser;
58 PetscCall(VecGetArrayRead(X, &x));
59 PetscCall(VecGetArray(F, &f));
60 PetscCall(EvaluateResidual(x, user->mu, f));
61 PetscCall(VecRestoreArrayRead(X, &x));
62 PetscCall(VecRestoreArray(F, &f));
63 PetscFunctionReturn(PETSC_SUCCESS);
64 }
65
66 /*
67 Compute the Jacobian w.r.t. x using tapeless mode of ADOL-C.
68 */
RHSJacobian(TS ts,PetscReal t,Vec X,Mat A,Mat B,PetscCtx ctx)69 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat A, Mat B, PetscCtx ctx)
70 {
71 User user = (User)ctx;
72 PetscScalar **J;
73 const PetscScalar *x;
74 adouble f_a[2]; /* 'active' double for dependent variables */
75 adouble x_a[2], mu_a; /* 'active' doubles for independent variables */
76 PetscInt i, j;
77
78 PetscFunctionBeginUser;
79 /* Set values for independent variables and parameters */
80 PetscCall(VecGetArrayRead(X, &x));
81 x_a[0].setValue(x[0]);
82 x_a[1].setValue(x[1]);
83 mu_a.setValue(user->mu);
84 PetscCall(VecRestoreArrayRead(X, &x));
85
86 /* Set seed matrix as 3x3 identity matrix */
87 x_a[0].setADValue(0, 1.);
88 x_a[0].setADValue(1, 0.);
89 x_a[0].setADValue(2, 0.);
90 x_a[1].setADValue(0, 0.);
91 x_a[1].setADValue(1, 1.);
92 x_a[1].setADValue(2, 0.);
93 mu_a.setADValue(0, 0.);
94 mu_a.setADValue(1, 0.);
95 mu_a.setADValue(2, 1.);
96
97 /* Evaluate residual (on active variables) */
98 PetscCall(EvaluateResidual(x_a, mu_a, f_a));
99
100 /* Extract derivatives */
101 PetscCall(PetscMalloc1(user->adctx->n, &J));
102 J[0] = (PetscScalar *)f_a[0].getADValue();
103 J[1] = (PetscScalar *)f_a[1].getADValue();
104
105 /* Set matrix values */
106 for (i = 0; i < user->adctx->m; i++) {
107 for (j = 0; j < user->adctx->n; j++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES));
108 }
109 PetscCall(PetscFree(J));
110 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
111 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
112 if (A != B) {
113 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
114 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
115 }
116 PetscFunctionReturn(PETSC_SUCCESS);
117 }
118
119 /*
120 Compute the Jacobian w.r.t. mu using tapeless mode of ADOL-C.
121 */
RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,PetscCtx ctx)122 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx)
123 {
124 User user = (User)ctx;
125 PetscScalar **J;
126 PetscScalar *x;
127 adouble f_a[2]; /* 'active' double for dependent variables */
128 adouble x_a[2], mu_a; /* 'active' doubles for independent variables */
129 PetscInt i, j = 0;
130
131 PetscFunctionBeginUser;
132 /* Set values for independent variables and parameters */
133 PetscCall(VecGetArray(X, &x));
134 x_a[0].setValue(x[0]);
135 x_a[1].setValue(x[1]);
136 mu_a.setValue(user->mu);
137 PetscCall(VecRestoreArray(X, &x));
138
139 /* Set seed matrix as 3x3 identity matrix */
140 x_a[0].setADValue(0, 1.);
141 x_a[0].setADValue(1, 0.);
142 x_a[0].setADValue(2, 0.);
143 x_a[1].setADValue(0, 0.);
144 x_a[1].setADValue(1, 1.);
145 x_a[1].setADValue(2, 0.);
146 mu_a.setADValue(0, 0.);
147 mu_a.setADValue(1, 0.);
148 mu_a.setADValue(2, 1.);
149
150 /* Evaluate residual (on active variables) */
151 PetscCall(EvaluateResidual(x_a, mu_a, f_a));
152
153 /* Extract derivatives */
154 PetscCall(PetscMalloc1(2, &J));
155 J[0] = (PetscScalar *)f_a[0].getADValue();
156 J[1] = (PetscScalar *)f_a[1].getADValue();
157
158 /* Set matrix values */
159 for (i = 0; i < user->adctx->m; i++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][user->adctx->n], INSERT_VALUES));
160 PetscCall(PetscFree(J));
161 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
162 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
163 PetscFunctionReturn(PETSC_SUCCESS);
164 }
165
166 /*
167 Monitor timesteps and use interpolation to output at integer multiples of 0.1
168 */
Monitor(TS ts,PetscInt step,PetscReal t,Vec X,PetscCtx ctx)169 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, PetscCtx ctx)
170 {
171 const PetscScalar *x;
172 PetscReal tfinal, dt, tprev;
173 User user = (User)ctx;
174
175 PetscFunctionBeginUser;
176 PetscCall(TSGetTimeStep(ts, &dt));
177 PetscCall(TSGetMaxTime(ts, &tfinal));
178 PetscCall(TSGetPrevTime(ts, &tprev));
179 PetscCall(VecGetArrayRead(X, &x));
180 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
181 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "t %.6f (tprev = %.6f) \n", (double)t, (double)tprev));
182 PetscCall(VecRestoreArrayRead(X, &x));
183 PetscFunctionReturn(PETSC_SUCCESS);
184 }
185
main(int argc,char ** argv)186 int main(int argc, char **argv)
187 {
188 TS ts; /* nonlinear solver */
189 Vec x; /* solution, residual vectors */
190 Mat A; /* Jacobian matrix */
191 Mat Jacp; /* JacobianP matrix */
192 PetscInt steps;
193 PetscReal ftime = 0.5;
194 PetscBool monitor = PETSC_FALSE;
195 PetscScalar *x_ptr;
196 PetscMPIInt size;
197 struct _n_User user;
198 AdolcCtx *adctx;
199 Vec lambda[2], mu[2];
200
201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202 Initialize program
203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204 PetscFunctionBeginUser;
205 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
206 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
207 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
208
209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210 Set runtime options and create AdolcCtx
211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212 PetscCall(PetscNew(&adctx));
213 user.mu = 1;
214 user.next_output = 0.0;
215 adctx->m = 2;
216 adctx->n = 2;
217 adctx->p = 2;
218 user.adctx = adctx;
219 adtl::setNumDir(adctx->n + 1); /* #indep. variables, plus parameters */
220
221 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
222 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
223
224 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 Create necessary matrix and vectors, solve same ODE on every process
226 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
228 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
229 PetscCall(MatSetFromOptions(A));
230 PetscCall(MatSetUp(A));
231 PetscCall(MatCreateVecs(A, &x, NULL));
232
233 PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
234 PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
235 PetscCall(MatSetFromOptions(Jacp));
236 PetscCall(MatSetUp(Jacp));
237
238 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
239 Create timestepping solver context
240 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
241 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
242 PetscCall(TSSetType(ts, TSRK));
243 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionPassive, &user));
244
245 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
246 Set initial conditions
247 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
248 PetscCall(VecGetArray(x, &x_ptr));
249 x_ptr[0] = 2;
250 x_ptr[1] = 0.66666654321;
251 PetscCall(VecRestoreArray(x, &x_ptr));
252
253 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
254 Set RHS Jacobian for the adjoint integration
255 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
256 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, &user));
257 PetscCall(TSSetMaxTime(ts, ftime));
258 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
259 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
260 PetscCall(TSSetTimeStep(ts, .001));
261
262 /*
263 Have the TS save its trajectory so that TSAdjointSolve() may be used
264 */
265 PetscCall(TSSetSaveTrajectory(ts));
266
267 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
268 Set runtime options
269 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
270 PetscCall(TSSetFromOptions(ts));
271
272 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273 Solve nonlinear system
274 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275 PetscCall(TSSolve(ts, x));
276 PetscCall(TSGetSolveTime(ts, &ftime));
277 PetscCall(TSGetStepNumber(ts, &steps));
278 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, steps, (double)ftime));
279 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD));
280
281 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
282 Start the Adjoint model
283 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
284 PetscCall(MatCreateVecs(A, &lambda[0], NULL));
285 PetscCall(MatCreateVecs(A, &lambda[1], NULL));
286 /* Reset initial conditions for the adjoint integration */
287 PetscCall(VecGetArray(lambda[0], &x_ptr));
288 x_ptr[0] = 1.0;
289 x_ptr[1] = 0.0;
290 PetscCall(VecRestoreArray(lambda[0], &x_ptr));
291 PetscCall(VecGetArray(lambda[1], &x_ptr));
292 x_ptr[0] = 0.0;
293 x_ptr[1] = 1.0;
294 PetscCall(VecRestoreArray(lambda[1], &x_ptr));
295
296 PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
297 PetscCall(MatCreateVecs(Jacp, &mu[1], NULL));
298 PetscCall(VecGetArray(mu[0], &x_ptr));
299 x_ptr[0] = 0.0;
300 PetscCall(VecRestoreArray(mu[0], &x_ptr));
301 PetscCall(VecGetArray(mu[1], &x_ptr));
302 x_ptr[0] = 0.0;
303 PetscCall(VecRestoreArray(mu[1], &x_ptr));
304 PetscCall(TSSetCostGradients(ts, 2, lambda, mu));
305
306 /* Set RHS JacobianP */
307 PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &user));
308
309 PetscCall(TSAdjointSolve(ts));
310
311 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
312 PetscCall(VecView(lambda[1], PETSC_VIEWER_STDOUT_WORLD));
313 PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
314 PetscCall(VecView(mu[1], PETSC_VIEWER_STDOUT_WORLD));
315
316 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
317 Free work space. All PETSc objects should be destroyed when they
318 are no longer needed.
319 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320 PetscCall(MatDestroy(&A));
321 PetscCall(MatDestroy(&Jacp));
322 PetscCall(VecDestroy(&x));
323 PetscCall(VecDestroy(&lambda[0]));
324 PetscCall(VecDestroy(&lambda[1]));
325 PetscCall(VecDestroy(&mu[0]));
326 PetscCall(VecDestroy(&mu[1]));
327 PetscCall(TSDestroy(&ts));
328 PetscCall(PetscFree(adctx));
329 PetscCall(PetscFinalize());
330 return 0;
331 }
332
333 /*TEST
334
335 build:
336 requires: double !complex adolc
337
338 test:
339 suffix: 1
340 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor
341 output_file: output/ex16adj_tl_1.out
342
343 test:
344 suffix: 2
345 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -mu 5
346 output_file: output/ex16adj_tl_2.out
347
348 test:
349 suffix: 3
350 args: -ts_max_steps 10 -monitor
351 output_file: output/ex16adj_tl_3.out
352
353 test:
354 suffix: 4
355 args: -ts_max_steps 10 -monitor -mu 5
356 output_file: output/ex16adj_tl_4.out
357
358 TEST*/
359