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