Actual source code: ex20td.c
1: static char help[] = "Performs adjoint sensitivity analysis for a van der Pol like \n\
2: equation with time dependent parameters using two approaches : \n\
3: track : track only local sensitivities at each adjoint step \n\
4: and accumulate them in a global array \n\
5: global : track parameters at all timesteps together \n\
6: Choose one of the two at runtime by -sa_method {track,global}. \n";
8: /*
9: Concepts: TS^adjoint for time dependent parameters
10: Concepts: TS^Customized adjoint monitor based sensitivity tracking
11: Concepts: TS^All at once approach to sensitivity tracking
12: Processors: 1
13: */
15: /*
16: Simple example to demonstrate TSAdjoint capabilities for time dependent params
17: without integral cost terms using either a tracking or global method.
19: Modify the Van Der Pol Eq to :
20: [u1'] = [mu1(t)*u1]
21: [u2'] = [mu2(t)*((1-u1^2)*u2-u1)]
22: (with initial conditions & params independent)
24: Define uref to be solution with initail conditions (2,-2/3), mu=(1,1e3)
25: - u_ref : (1.5967,-1.02969)
27: Define const function as cost = 2-norm(u - u_ref);
29: Initialization for the adjoint TS :
30: - dcost/dy|final_time = 2*(u-u_ref)|final_time
31: - dcost/dp|final_time = 0
33: The tracking method only tracks local sensitivity at each time step
34: and accumulates these sensitivities in a global array. Since the structure
35: of the equations being solved at each time step does not change, the jacobian
36: wrt parameters is defined analogous to constant RHSJacobian for a liner
37: TSSolve and the size of the jacP is independent of the number of time
38: steps. Enable this mode of adjoint analysis by -sa_method track.
40: The global method combines the parameters at all timesteps and tracks them
41: together. Thus, the columns of the jacP matrix are filled dependent upon the
42: time step. Also, the dimensions of the jacP matrix now depend upon the number
43: of time steps. Enable this mode of adjoint analysis by -sa_method global.
45: Since the equations here have parameters at predefined time steps, this
46: example should be run with non adaptive time stepping solvers only. This
47: can be ensured by -ts_adapt_type none (which is the default behavior only
48: for certain TS solvers like TSCN. If using an explicit TS like TSRK,
49: please be sure to add the aforementioned option to disable adaptive
50: timestepping.)
51: */
53: /*
54: Include "petscts.h" so that we can use TS solvers. Note that this file
55: automatically includes:
56: petscsys.h - base PETSc routines petscvec.h - vectors
57: petscmat.h - matrices
58: petscis.h - index sets petscksp.h - Krylov subspace methods
59: petscviewer.h - viewers petscpc.h - preconditioners
60: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
61: */
62: #include <petscts.h>
64: extern PetscErrorCode RHSFunction(TS ,PetscReal ,Vec ,Vec ,void *);
65: extern PetscErrorCode RHSJacobian(TS ,PetscReal ,Vec ,Mat ,Mat ,void *);
66: extern PetscErrorCode RHSJacobianP_track(TS ,PetscReal ,Vec ,Mat ,void *);
67: extern PetscErrorCode RHSJacobianP_global(TS ,PetscReal ,Vec ,Mat ,void *);
68: extern PetscErrorCode Monitor(TS ,PetscInt ,PetscReal ,Vec ,void *);
69: extern PetscErrorCode AdjointMonitor(TS ,PetscInt ,PetscReal ,Vec ,PetscInt ,Vec *, Vec *,void *);
71: /*
72: User-defined application context - contains data needed by the
73: application-provided call-back routines.
74: */
76: typedef struct {
77: /*------------- Forward solve data structures --------------*/
78: PetscInt max_steps; /* number of steps to run ts for */
79: PetscReal final_time; /* final time to integrate to*/
80: PetscReal time_step; /* ts integration time step */
81: Vec mu1; /* time dependent params */
82: Vec mu2; /* time dependent params */
83: Vec U; /* solution vector */
84: Mat A; /* Jacobian matrix */
86: /*------------- Adjoint solve data structures --------------*/
87: Mat Jacp; /* JacobianP matrix */
88: Vec lambda; /* adjoint variable */
89: Vec mup; /* adjoint variable */
91: /*------------- Global accumation vecs for monitor based tracking --------------*/
92: Vec sens_mu1; /* global sensitivity accumulation */
93: Vec sens_mu2; /* global sensitivity accumulation */
94: PetscInt adj_idx; /* to keep track of adjoint solve index */
95: } AppCtx;
97: typedef enum {SA_TRACK, SA_GLOBAL} SAMethod;
98: static const char *const SAMethods[] = {"TRACK","GLOBAL","SAMethod","SA_",0};
100: /* ----------------------- Explicit form of the ODE -------------------- */
102: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
103: {
104: AppCtx *user = (AppCtx*) ctx;
105: PetscScalar *f;
106: PetscInt curr_step;
107: const PetscScalar *u;
108: const PetscScalar *mu1;
109: const PetscScalar *mu2;
112: TSGetStepNumber(ts,&curr_step);
113: VecGetArrayRead(U,&u);
114: VecGetArrayRead(user->mu1,&mu1);
115: VecGetArrayRead(user->mu2,&mu2);
116: VecGetArray(F,&f);
117: f[0] = mu1[curr_step]*u[1];
118: f[1] = mu2[curr_step]*((1.-u[0]*u[0])*u[1]-u[0]);
119: VecRestoreArrayRead(U,&u);
120: VecRestoreArrayRead(user->mu1,&mu1);
121: VecRestoreArrayRead(user->mu2,&mu2);
122: VecRestoreArray(F,&f);
123: return 0;
124: }
126: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
127: {
128: AppCtx *user = (AppCtx*) ctx;
129: PetscInt rowcol[] = {0,1};
130: PetscScalar J[2][2];
131: PetscInt curr_step;
132: const PetscScalar *u;
133: const PetscScalar *mu1;
134: const PetscScalar *mu2;
137: TSGetStepNumber(ts,&curr_step);
138: VecGetArrayRead(user->mu1,&mu1);
139: VecGetArrayRead(user->mu2,&mu2);
140: VecGetArrayRead(U,&u);
141: J[0][0] = 0;
142: J[1][0] = -mu2[curr_step]*(2.0*u[1]*u[0]+1.);
143: J[0][1] = mu1[curr_step];
144: J[1][1] = mu2[curr_step]*(1.0-u[0]*u[0]);
145: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
146: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
147: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
148: VecRestoreArrayRead(U,&u);
149: VecRestoreArrayRead(user->mu1,&mu1);
150: VecRestoreArrayRead(user->mu2,&mu2);
151: return 0;
152: }
154: /* ------------------ Jacobian wrt parameters for tracking method ------------------ */
156: PetscErrorCode RHSJacobianP_track(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
157: {
158: PetscInt row[] = {0,1},col[] = {0,1};
159: PetscScalar J[2][2];
160: const PetscScalar *u;
163: VecGetArrayRead(U,&u);
164: J[0][0] = u[1];
165: J[1][0] = 0;
166: J[0][1] = 0;
167: J[1][1] = (1.-u[0]*u[0])*u[1]-u[0];
168: MatSetValues(A,2,row,2,col,&J[0][0],INSERT_VALUES);
169: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
170: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
171: VecRestoreArrayRead(U,&u);
172: return 0;
173: }
175: /* ------------------ Jacobian wrt parameters for global method ------------------ */
177: PetscErrorCode RHSJacobianP_global(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
178: {
179: PetscInt row[] = {0,1},col[] = {0,1};
180: PetscScalar J[2][2];
181: const PetscScalar *u;
182: PetscInt curr_step;
185: TSGetStepNumber(ts,&curr_step);
186: VecGetArrayRead(U,&u);
187: J[0][0] = u[1];
188: J[1][0] = 0;
189: J[0][1] = 0;
190: J[1][1] = (1.-u[0]*u[0])*u[1]-u[0];
191: col[0] = (curr_step)*2;
192: col[1] = (curr_step)*2+1;
193: MatSetValues(A,2,row,2,col,&J[0][0],INSERT_VALUES);
194: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
195: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
196: VecRestoreArrayRead(U,&u);
197: return 0;
198: }
200: /* Dump solution to console if called */
201: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
202: {
204: PetscPrintf(PETSC_COMM_WORLD,"\n Solution at time %e is \n", t);
205: VecView(U,PETSC_VIEWER_STDOUT_WORLD);
206: return 0;
207: }
209: /* Customized adjoint monitor to keep track of local
210: sensitivities by storing them in a global sensitivity array.
211: Note : This routine is only used for the tracking method. */
212: PetscErrorCode AdjointMonitor(TS ts,PetscInt steps,PetscReal time,Vec u,PetscInt numcost,Vec *lambda, Vec *mu,void *ctx)
213: {
214: AppCtx *user = (AppCtx*) ctx;
215: PetscInt curr_step;
216: PetscScalar *sensmu1_glob;
217: PetscScalar *sensmu2_glob;
218: const PetscScalar *sensmu_loc;
221: TSGetStepNumber(ts,&curr_step);
222: /* Note that we skip the first call to the monitor in the adjoint
223: solve since the sensitivities are already set (during
224: initialization of adjoint vectors).
225: We also note that each indvidial TSAdjointSolve calls the monitor
226: twice, once at the step it is integrating from and once at the step
227: it integrates to. Only the second call is useful for transferring
228: local sensitivities to the global array. */
229: if (curr_step == user->adj_idx) {
230: return 0;
231: } else {
232: VecGetArrayRead(*mu,&sensmu_loc);
233: VecGetArray(user->sens_mu1,&sensmu1_glob);
234: VecGetArray(user->sens_mu2,&sensmu2_glob);
235: sensmu1_glob[curr_step] = sensmu_loc[0];
236: sensmu2_glob[curr_step] = sensmu_loc[1];
237: VecRestoreArray(user->sens_mu1,&sensmu1_glob);
238: VecRestoreArray(user->sens_mu2,&sensmu2_glob);
239: VecRestoreArrayRead(*mu,&sensmu_loc);
240: return 0;
241: }
242: }
244: int main(int argc,char **argv)
245: {
246: TS ts;
247: AppCtx user;
248: PetscScalar *x_ptr,*y_ptr,*u_ptr;
249: PetscMPIInt size;
250: PetscBool monitor = PETSC_FALSE;
251: SAMethod sa = SA_GLOBAL;
254: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
255: Initialize program
256: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
257: PetscInitialize(&argc,&argv,NULL,help);
258: MPI_Comm_size(PETSC_COMM_WORLD,&size);
261: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
262: Set runtime options
263: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
264: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"SA analysis options.","");{
265: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
266: PetscOptionsEnum("-sa_method","Sensitivity analysis method (track or global)","",SAMethods,(PetscEnum)sa,(PetscEnum*)&sa,NULL);
267: }
268: PetscOptionsEnd();
270: user.final_time = 0.1;
271: user.max_steps = 5;
272: user.time_step = user.final_time/user.max_steps;
274: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: Create necessary matrix and vectors for forward solve.
276: Create Jacp matrix for adjoint solve.
277: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278: VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.mu1);
279: VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.mu2);
280: VecSet(user.mu1,1.25);
281: VecSet(user.mu2,1.0e2);
283: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
284: For tracking method : create the global sensitivity array to
285: accumulate sensitivity with respect to parameters at each step.
286: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
287: if (sa == SA_TRACK) {
288: VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.sens_mu1);
289: VecCreateSeq(PETSC_COMM_WORLD,user.max_steps,&user.sens_mu2);
290: }
292: MatCreate(PETSC_COMM_WORLD,&user.A);
293: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
294: MatSetFromOptions(user.A);
295: MatSetUp(user.A);
296: MatCreateVecs(user.A,&user.U,NULL);
298: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
299: Note that the dimensions of the Jacp matrix depend upon the
300: sensitivity analysis method being used !
301: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
302: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
303: if (sa == SA_TRACK) {
304: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2);
305: }
306: if (sa == SA_GLOBAL) {
307: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,user.max_steps*2);
308: }
309: MatSetFromOptions(user.Jacp);
310: MatSetUp(user.Jacp);
312: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313: Create timestepping solver context
314: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
315: TSCreate(PETSC_COMM_WORLD,&ts);
316: TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT);
317: TSSetType(ts,TSCN);
319: TSSetRHSFunction(ts,NULL,RHSFunction,&user);
320: TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user);
321: if (sa == SA_TRACK) {
322: TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP_track,&user);
323: }
324: if (sa == SA_GLOBAL) {
325: TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP_global,&user);
326: }
328: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
329: TSSetMaxTime(ts,user.final_time);
330: TSSetTimeStep(ts,user.final_time/user.max_steps);
332: if (monitor) {
333: TSMonitorSet(ts,Monitor,&user,NULL);
334: }
335: if (sa == SA_TRACK) {
336: TSAdjointMonitorSet(ts,AdjointMonitor,&user,NULL);
337: }
339: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
340: Set initial conditions
341: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
342: VecGetArray(user.U,&x_ptr);
343: x_ptr[0] = 2.0;
344: x_ptr[1] = -2.0/3.0;
345: VecRestoreArray(user.U,&x_ptr);
347: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
348: Save trajectory of solution so that TSAdjointSolve() may be used
349: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
350: TSSetSaveTrajectory(ts);
352: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
353: Set runtime options
354: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
355: TSSetFromOptions(ts);
357: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
358: Execute forward model and print solution.
359: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
360: TSSolve(ts,user.U);
361: PetscPrintf(PETSC_COMM_WORLD,"\n Solution of forward TS :\n");
362: VecView(user.U,PETSC_VIEWER_STDOUT_WORLD);
363: PetscPrintf(PETSC_COMM_WORLD,"\n Forward TS solve successfull! Adjoint run begins!\n");
365: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
366: Adjoint model starts here! Create adjoint vectors.
367: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
368: MatCreateVecs(user.A,&user.lambda,NULL);
369: MatCreateVecs(user.Jacp,&user.mup,NULL);
371: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
372: Set initial conditions for the adjoint vector
373: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
374: VecGetArray(user.U,&u_ptr);
375: VecGetArray(user.lambda,&y_ptr);
376: y_ptr[0] = 2*(u_ptr[0] - 1.5967);
377: y_ptr[1] = 2*(u_ptr[1] - -(1.02969));
378: VecRestoreArray(user.lambda,&y_ptr);
379: VecRestoreArray(user.U,&y_ptr);
380: VecSet(user.mup,0);
382: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
383: Set number of cost functions.
384: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
385: TSSetCostGradients(ts,1,&user.lambda,&user.mup);
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: The adjoint vector mup has to be reset for each adjoint step when
389: using the tracking method as we want to treat the parameters at each
390: time step one at a time and prevent accumulation of the sensitivities
391: from parameters at previous time steps.
392: This is not necessary for the global method as each time dependent
393: parameter is treated as an independent parameter.
394: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
395: if (sa == SA_TRACK) {
396: for (user.adj_idx=user.max_steps; user.adj_idx>0; user.adj_idx--) {
397: VecSet(user.mup,0);
398: TSAdjointSetSteps(ts, 1);
399: TSAdjointSolve(ts);
400: }
401: }
402: if (sa == SA_GLOBAL) {
403: TSAdjointSolve(ts);
404: }
406: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407: Dispaly adjoint sensitivities wrt parameters and initial conditions
408: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
409: if (sa == SA_TRACK) {
410: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt mu1: d[cost]/d[mu1]\n");
411: VecView(user.sens_mu1,PETSC_VIEWER_STDOUT_WORLD);
412: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt mu2: d[cost]/d[mu2]\n");
413: VecView(user.sens_mu2,PETSC_VIEWER_STDOUT_WORLD);
414: }
416: if (sa == SA_GLOBAL) {
417: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt params: d[cost]/d[p], where p refers to \n\
418: the interlaced vector made by combining mu1,mu2\n");
419: VecView(user.mup,PETSC_VIEWER_STDOUT_WORLD);
420: }
422: PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[cost]/d[u(t=0)]\n");
423: VecView(user.lambda,PETSC_VIEWER_STDOUT_WORLD);
425: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
426: Free work space!
427: All PETSc objects should be destroyed when they are no longer needed.
428: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
429: MatDestroy(&user.A);
430: MatDestroy(&user.Jacp);
431: VecDestroy(&user.U);
432: VecDestroy(&user.lambda);
433: VecDestroy(&user.mup);
434: VecDestroy(&user.mu1);
435: VecDestroy(&user.mu2);
436: if (sa == SA_TRACK) {
437: VecDestroy(&user.sens_mu1);
438: VecDestroy(&user.sens_mu2);
439: }
440: TSDestroy(&ts);
441: PetscFinalize();
442: return(ierr);
443: }
445: /*TEST
447: test:
448: requires: !complex
449: suffix : track
450: args : -sa_method track
452: test:
453: requires: !complex
454: suffix : global
455: args : -sa_method global
457: TEST*/