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