Actual source code: ex20adj.c
1: static char help[] = "Performs adjoint sensitivity analysis for the van der Pol equation.\n";
3: /*
4: Concepts: TS^time-dependent nonlinear problems
5: Concepts: TS^van der Pol equation DAE equivalent
6: Concepts: TS^adjoint sensitivity analysis
7: Processors: 1
8: */
9: /* ------------------------------------------------------------------------
11: This program solves the van der Pol DAE ODE equivalent
12: [ u_1' ] = [ u_2 ] (2)
13: [ u_2' ] [ \mu ((1 - u_1^2) u_2 - u_1) ]
14: on the domain 0 <= x <= 1, with the boundary conditions
15: u_1(0) = 2, u_2(0) = - 2/3 +10/(81*\mu) - 292/(2187*\mu^2),
16: and
17: \mu = 10^6 ( y'(0) ~ -0.6666665432100101).,
18: and computes the sensitivities of the final solution w.r.t. initial conditions and parameter \mu with the implicit theta method and its discrete adjoint.
20: Notes:
21: This code demonstrates the TSAdjoint interface to a DAE system.
23: The user provides the implicit right-hand-side function
24: [ F(u',u,t) ] = [u' - f(u,t)] = [ u_1'] - [ u_2 ]
25: [ u_2'] [ \mu ((1-u_1^2)u_2-u_1) ]
27: and the Jacobian of F (from the PETSc user manual)
29: dF dF
30: J(F) = a * -- + --
31: du' du
33: and the JacobianP of the explicit right-hand side of (2) f(u,t) ( which is equivalent to -F(0,u,t)).
34: df [ 0 ]
35: -- = [ ]
36: dp [ (1 - u_1^2) u_2 - u_1 ].
38: See ex20.c for more details on the Jacobian.
40: ------------------------------------------------------------------------- */
41: #include <petscts.h>
42: #include <petsctao.h>
44: typedef struct _n_User *User;
45: struct _n_User {
46: PetscReal mu;
47: PetscReal next_output;
49: /* Sensitivity analysis support */
50: PetscInt steps;
51: PetscReal ftime;
52: Mat A; /* Jacobian matrix */
53: Mat Jacp; /* JacobianP matrix */
54: Vec U,lambda[2],mup[2]; /* adjoint variables */
55: };
57: /* ----------------------- Explicit form of the ODE -------------------- */
59: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
60: {
61: User user = (User)ctx;
62: PetscScalar *f;
63: const PetscScalar *u;
66: VecGetArrayRead(U,&u);
67: VecGetArray(F,&f);
68: f[0] = u[1];
69: f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
70: VecRestoreArrayRead(U,&u);
71: VecRestoreArray(F,&f);
72: return 0;
73: }
75: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
76: {
77: User user = (User)ctx;
78: PetscReal mu = user->mu;
79: PetscInt rowcol[] = {0,1};
80: PetscScalar J[2][2];
81: const PetscScalar *u;
84: VecGetArrayRead(U,&u);
85: J[0][0] = 0;
86: J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
87: J[0][1] = 1.0;
88: J[1][1] = mu*(1.0-u[0]*u[0]);
89: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
90: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
92: if (A != B) {
93: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
94: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
95: }
96: VecRestoreArrayRead(U,&u);
97: return 0;
98: }
100: /* ----------------------- Implicit form of the ODE -------------------- */
102: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
103: {
104: User user = (User)ctx;
105: const PetscScalar *u,*udot;
106: PetscScalar *f;
109: VecGetArrayRead(U,&u);
110: VecGetArrayRead(Udot,&udot);
111: VecGetArray(F,&f);
112: f[0] = udot[0] - u[1];
113: f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]);
114: VecRestoreArrayRead(U,&u);
115: VecRestoreArrayRead(Udot,&udot);
116: VecRestoreArray(F,&f);
117: return 0;
118: }
120: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
121: {
122: User user = (User)ctx;
123: PetscInt rowcol[] = {0,1};
124: PetscScalar J[2][2];
125: const PetscScalar *u;
128: VecGetArrayRead(U,&u);
130: J[0][0] = a; J[0][1] = -1.0;
131: J[1][0] = user->mu*(2.0*u[0]*u[1] + 1.0); J[1][1] = a - user->mu*(1.0-u[0]*u[0]);
133: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
134: VecRestoreArrayRead(U,&u);
136: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
137: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
138: if (B && A != B) {
139: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
140: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
141: }
142: return 0;
143: }
145: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
146: {
147: PetscInt row[] = {0,1},col[]={0};
148: PetscScalar J[2][1];
149: const PetscScalar *u;
152: VecGetArrayRead(U,&u);
153: J[0][0] = 0;
154: J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
155: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
156: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
157: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
158: VecRestoreArrayRead(U,&u);
159: return 0;
160: }
162: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
163: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
164: {
165: const PetscScalar *u;
166: PetscReal tfinal, dt;
167: User user = (User)ctx;
168: Vec interpolatedU;
171: TSGetTimeStep(ts,&dt);
172: TSGetMaxTime(ts,&tfinal);
174: while (user->next_output <= t && user->next_output <= tfinal) {
175: VecDuplicate(U,&interpolatedU);
176: TSInterpolate(ts,user->next_output,interpolatedU);
177: VecGetArrayRead(interpolatedU,&u);
178: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
179: (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
180: (double)PetscRealPart(u[1])));
181: VecRestoreArrayRead(interpolatedU,&u);
182: VecDestroy(&interpolatedU);
183: user->next_output += 0.1;
184: }
185: return 0;
186: }
188: int main(int argc,char **argv)
189: {
190: TS ts;
191: PetscBool monitor = PETSC_FALSE,implicitform = PETSC_TRUE;
192: PetscScalar *x_ptr,*y_ptr,derp;
193: PetscMPIInt size;
194: struct _n_User user;
196: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
197: Initialize program
198: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
199: PetscInitialize(&argc,&argv,NULL,help);
200: MPI_Comm_size(PETSC_COMM_WORLD,&size);
203: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204: Set runtime options
205: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206: user.next_output = 0.0;
207: user.mu = 1.0e3;
208: user.steps = 0;
209: user.ftime = 0.5;
210: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
211: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
212: PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Create necessary matrix and vectors, solve same ODE on every process
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: MatCreate(PETSC_COMM_WORLD,&user.A);
218: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
219: MatSetFromOptions(user.A);
220: MatSetUp(user.A);
221: MatCreateVecs(user.A,&user.U,NULL);
223: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
224: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
225: MatSetFromOptions(user.Jacp);
226: MatSetUp(user.Jacp);
228: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
229: Create timestepping solver context
230: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
231: TSCreate(PETSC_COMM_WORLD,&ts);
232: TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
233: if (implicitform) {
234: TSSetIFunction(ts,NULL,IFunction,&user);
235: TSSetIJacobian(ts,user.A,user.A,IJacobian,&user);
236: TSSetType(ts,TSCN);
237: } else {
238: TSSetRHSFunction(ts,NULL,RHSFunction,&user);
239: TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user);
240: TSSetType(ts,TSRK);
241: }
242: TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user);
243: TSSetMaxTime(ts,user.ftime);
244: TSSetTimeStep(ts,0.001);
245: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
246: if (monitor) TSMonitorSet(ts,Monitor,&user,NULL);
248: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
249: Set initial conditions
250: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
251: VecGetArray(user.U,&x_ptr);
252: x_ptr[0] = 2.0;
253: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
254: VecRestoreArray(user.U,&x_ptr);
255: TSSetTimeStep(ts,0.001);
257: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: Save trajectory of solution so that TSAdjointSolve() may be used
259: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260: TSSetSaveTrajectory(ts);
262: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
263: Set runtime options
264: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
265: TSSetFromOptions(ts);
267: TSSolve(ts,user.U);
268: TSGetSolveTime(ts,&user.ftime);
269: TSGetStepNumber(ts,&user.steps);
271: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272: Adjoint model starts here
273: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: MatCreateVecs(user.A,&user.lambda[0],NULL);
275: /* Set initial conditions for the adjoint integration */
276: VecGetArray(user.lambda[0],&y_ptr);
277: y_ptr[0] = 1.0; y_ptr[1] = 0.0;
278: VecRestoreArray(user.lambda[0],&y_ptr);
279: MatCreateVecs(user.A,&user.lambda[1],NULL);
280: VecGetArray(user.lambda[1],&y_ptr);
281: y_ptr[0] = 0.0; y_ptr[1] = 1.0;
282: VecRestoreArray(user.lambda[1],&y_ptr);
284: MatCreateVecs(user.Jacp,&user.mup[0],NULL);
285: VecGetArray(user.mup[0],&x_ptr);
286: x_ptr[0] = 0.0;
287: VecRestoreArray(user.mup[0],&x_ptr);
288: MatCreateVecs(user.Jacp,&user.mup[1],NULL);
289: VecGetArray(user.mup[1],&x_ptr);
290: x_ptr[0] = 0.0;
291: VecRestoreArray(user.mup[1],&x_ptr);
293: TSSetCostGradients(ts,2,user.lambda,user.mup);
295: TSAdjointSolve(ts);
297: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[y(tf)]/d[y0] d[y(tf)]/d[z0]\n");
298: VecView(user.lambda[0],PETSC_VIEWER_STDOUT_WORLD);
299: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[z(tf)]/d[y0] d[z(tf)]/d[z0]\n");
300: VecView(user.lambda[1],PETSC_VIEWER_STDOUT_WORLD);
302: VecGetArray(user.mup[0],&x_ptr);
303: VecGetArray(user.lambda[0],&y_ptr);
304: derp = y_ptr[1]*(-10.0/(81.0*user.mu*user.mu)+2.0*292.0/(2187.0*user.mu*user.mu*user.mu))+x_ptr[0];
305: VecRestoreArray(user.mup[0],&x_ptr);
306: VecRestoreArray(user.lambda[0],&y_ptr);
307: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameters: d[y(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp));
309: VecGetArray(user.mup[1],&x_ptr);
310: VecGetArray(user.lambda[1],&y_ptr);
311: derp = y_ptr[1]*(-10.0/(81.0*user.mu*user.mu)+2.0*292.0/(2187.0*user.mu*user.mu*user.mu))+x_ptr[0];
312: VecRestoreArray(user.mup[1],&x_ptr);
313: VecRestoreArray(user.lambda[1],&y_ptr);
314: PetscPrintf(PETSC_COMM_WORLD,"\n sensivitity wrt parameters: d[z(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp));
316: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
317: Free work space. All PETSc objects should be destroyed when they
318: are no longer needed.
319: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320: MatDestroy(&user.A);
321: MatDestroy(&user.Jacp);
322: VecDestroy(&user.U);
323: VecDestroy(&user.lambda[0]);
324: VecDestroy(&user.lambda[1]);
325: VecDestroy(&user.mup[0]);
326: VecDestroy(&user.mup[1]);
327: TSDestroy(&ts);
329: PetscFinalize();
330: return 0;
331: }
333: /*TEST
335: test:
336: requires: revolve
337: args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -viewer_binary_skip_info -ts_dt 0.001 -mu 100000
339: test:
340: suffix: 2
341: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only
343: test:
344: suffix: 3
345: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only 0
346: output_file: output/ex20adj_2.out
348: test:
349: suffix: 4
350: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
351: output_file: output/ex20adj_2.out
353: test:
354: suffix: 5
355: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
356: output_file: output/ex20adj_2.out
358: test:
359: suffix: 6
360: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack 0
361: output_file: output/ex20adj_2.out
363: test:
364: suffix: 7
365: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
366: output_file: output/ex20adj_2.out
368: test:
369: suffix: 8
370: requires: revolve !cams
371: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_solution_only -ts_trajectory_monitor
372: output_file: output/ex20adj_3.out
374: test:
375: suffix: 9
376: requires: revolve !cams
377: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_solution_only 0 -ts_trajectory_monitor
378: output_file: output/ex20adj_4.out
380: test:
381: requires: revolve
382: suffix: 10
383: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_revolve_online -ts_trajectory_solution_only
384: output_file: output/ex20adj_2.out
386: test:
387: requires: revolve
388: suffix: 11
389: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 5 -ts_trajectory_revolve_online -ts_trajectory_solution_only 0
390: output_file: output/ex20adj_2.out
392: test:
393: suffix: 12
394: requires: revolve
395: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only
396: output_file: output/ex20adj_2.out
398: test:
399: suffix: 13
400: requires: revolve
401: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only 0
402: output_file: output/ex20adj_2.out
404: test:
405: suffix: 14
406: requires: revolve
407: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
408: output_file: output/ex20adj_2.out
410: test:
411: suffix: 15
412: requires: revolve
413: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack 0
414: output_file: output/ex20adj_2.out
416: test:
417: suffix: 16
418: requires: revolve
419: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
420: output_file: output/ex20adj_2.out
422: test:
423: suffix: 17
424: requires: revolve
425: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
426: output_file: output/ex20adj_2.out
428: test:
429: suffix: 18
430: requires: revolve
431: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only -ts_trajectory_save_stack
432: output_file: output/ex20adj_2.out
434: test:
435: suffix: 19
436: requires: revolve
437: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack
438: output_file: output/ex20adj_2.out
440: test:
441: suffix: 20
442: requires: revolve
443: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_solution_only 0
444: output_file: output/ex20adj_2.out
446: test:
447: suffix: 21
448: requires: revolve
449: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_cps_ram 3 -ts_trajectory_max_cps_disk 8 -ts_trajectory_stride 5 -ts_trajectory_solution_only 0 -ts_trajectory_save_stack 0
450: output_file: output/ex20adj_2.out
452: test:
453: suffix: 22
454: args: -ts_type beuler -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only
455: output_file: output/ex20adj_2.out
457: test:
458: suffix: 23
459: requires: cams
460: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_units_ram 5 -ts_trajectory_solution_only -ts_trajectory_monitor -ts_trajectory_memory_type cams
461: output_file: output/ex20adj_5.out
463: test:
464: suffix: 24
465: requires: cams
466: args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_max_units_ram 5 -ts_trajectory_solution_only 0 -ts_trajectory_monitor -ts_trajectory_memory_type cams
467: output_file: output/ex20adj_6.out
469: TEST*/