Actual source code: ex20adj.c

petsc-master 2020-11-25
Report Typos and Errors
  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:   PetscErrorCode    ierr;
62:   User              user = (User)ctx;
63:   PetscScalar       *f;
64:   const PetscScalar *u;

68:   VecGetArray(F,&f);
69:   f[0] = u[1];
70:   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
72:   VecRestoreArray(F,&f);
73:   return(0);
74: }

76: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
77: {
78:   PetscErrorCode    ierr;
79:   User              user = (User)ctx;
80:   PetscReal         mu   = user->mu;
81:   PetscInt          rowcol[] = {0,1};
82:   PetscScalar       J[2][2];
83:   const PetscScalar *u;

87:   J[0][0] = 0;
88:   J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
89:   J[0][1] = 1.0;
90:   J[1][1] = mu*(1.0-u[0]*u[0]);
91:   MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
92:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
93:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
94:   if (A != B) {
95:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
96:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
97:   }
99:   return(0);
100: }

102: /* ----------------------- Implicit form of the ODE  -------------------- */

104: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
105: {
106:   PetscErrorCode    ierr;
107:   User              user = (User)ctx;
108:   const PetscScalar *u,*udot;
109:   PetscScalar       *f;

114:   VecGetArray(F,&f);
115:   f[0] = udot[0] - u[1];
116:   f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]);
119:   VecRestoreArray(F,&f);
120:   return(0);
121: }

123: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
124: {
125:   PetscErrorCode    ierr;
126:   User              user     = (User)ctx;
127:   PetscInt          rowcol[] = {0,1};
128:   PetscScalar       J[2][2];
129:   const PetscScalar *u;

134:   J[0][0] = a;     J[0][1] =  -1.0;
135:   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]);

137:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);

140:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
141:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
142:   if (B && A != B) {
143:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
144:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
145:   }
146:   return(0);
147: }

149: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
150: {
151:   PetscErrorCode    ierr;
152:   PetscInt          row[] = {0,1},col[]={0};
153:   PetscScalar       J[2][1];
154:   const PetscScalar *u;

158:   J[0][0] = 0;
159:   J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
160:   MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
161:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
162:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
164:   return(0);
165: }

167: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
168: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx)
169: {
170:   PetscErrorCode    ierr;
171:   const PetscScalar *u;
172:   PetscReal         tfinal, dt;
173:   User              user = (User)ctx;
174:   Vec               interpolatedU;

177:   TSGetTimeStep(ts,&dt);
178:   TSGetMaxTime(ts,&tfinal);

180:   while (user->next_output <= t && user->next_output <= tfinal) {
181:     VecDuplicate(U,&interpolatedU);
182:     TSInterpolate(ts,user->next_output,interpolatedU);
184:     PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
185:                        (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]),
186:                        (double)PetscRealPart(u[1]));
188:     VecDestroy(&interpolatedU);
189:     user->next_output += 0.1;
190:   }
191:   return(0);
192: }

194: int main(int argc,char **argv)
195: {
196:   TS             ts;
197:   PetscBool      monitor = PETSC_FALSE,implicitform = PETSC_TRUE;
198:   PetscScalar    *x_ptr,*y_ptr,derp;
199:   PetscMPIInt    size;
200:   struct _n_User user;

203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204:      Initialize program
205:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206:   PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
207:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
208:   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");

210:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
211:     Set runtime options
212:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
213:   user.next_output = 0.0;
214:   user.mu          = 1.0e3;
215:   user.steps       = 0;
216:   user.ftime       = 0.5;
217:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
218:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
219:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);

221:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
222:     Create necessary matrix and vectors, solve same ODE on every process
223:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224:   MatCreate(PETSC_COMM_WORLD,&user.A);
225:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
226:   MatSetFromOptions(user.A);
227:   MatSetUp(user.A);
228:   MatCreateVecs(user.A,&user.U,NULL);

230:   MatCreate(PETSC_COMM_WORLD,&user.Jacp);
231:   MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
232:   MatSetFromOptions(user.Jacp);
233:   MatSetUp(user.Jacp);

235:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
236:      Create timestepping solver context
237:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238:   TSCreate(PETSC_COMM_WORLD,&ts);
239:   TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
240:   if (implicitform) {
241:     TSSetIFunction(ts,NULL,IFunction,&user);
242:     TSSetIJacobian(ts,user.A,user.A,IJacobian,&user);
243:     TSSetType(ts,TSCN);
244:   } else {
245:     TSSetRHSFunction(ts,NULL,RHSFunction,&user);
246:     TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user);
247:     TSSetType(ts,TSRK);
248:   }
249:   TSSetRHSJacobianP(ts,user.Jacp,RHSJacobianP,&user);
250:   TSSetMaxTime(ts,user.ftime);
251:   TSSetTimeStep(ts,0.001);
252:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
253:   if (monitor) {
254:     TSMonitorSet(ts,Monitor,&user,NULL);
255:   }

257:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258:      Set initial conditions
259:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260:   VecGetArray(user.U,&x_ptr);
261:   x_ptr[0] = 2.0;
262:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
263:   VecRestoreArray(user.U,&x_ptr);
264:   TSSetTimeStep(ts,0.001);

266:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
267:     Save trajectory of solution so that TSAdjointSolve() may be used
268:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
269:   TSSetSaveTrajectory(ts);

271:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272:      Set runtime options
273:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274:   TSSetFromOptions(ts);

276:   TSSolve(ts,user.U);
277:   TSGetSolveTime(ts,&user.ftime);
278:   TSGetStepNumber(ts,&user.steps);

280:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
281:      Adjoint model starts here
282:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
283:   MatCreateVecs(user.A,&user.lambda[0],NULL);
284:   /* Set initial conditions for the adjoint integration */
285:   VecGetArray(user.lambda[0],&y_ptr);
286:   y_ptr[0] = 1.0; y_ptr[1] = 0.0;
287:   VecRestoreArray(user.lambda[0],&y_ptr);
288:   MatCreateVecs(user.A,&user.lambda[1],NULL);
289:   VecGetArray(user.lambda[1],&y_ptr);
290:   y_ptr[0] = 0.0; y_ptr[1] = 1.0;
291:   VecRestoreArray(user.lambda[1],&y_ptr);

293:   MatCreateVecs(user.Jacp,&user.mup[0],NULL);
294:   VecGetArray(user.mup[0],&x_ptr);
295:   x_ptr[0] = 0.0;
296:   VecRestoreArray(user.mup[0],&x_ptr);
297:   MatCreateVecs(user.Jacp,&user.mup[1],NULL);
298:   VecGetArray(user.mup[1],&x_ptr);
299:   x_ptr[0] = 0.0;
300:   VecRestoreArray(user.mup[1],&x_ptr);

306:   PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[y(tf)]/d[y0]  d[y(tf)]/d[z0]\n");
307:   VecView(user.lambda[0],PETSC_VIEWER_STDOUT_WORLD);
308:   PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt initial conditions: d[z(tf)]/d[y0]  d[z(tf)]/d[z0]\n");
309:   VecView(user.lambda[1],PETSC_VIEWER_STDOUT_WORLD);

311:   VecGetArray(user.mup[0],&x_ptr);
312:   VecGetArray(user.lambda[0],&y_ptr);
313:   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];
314:   VecRestoreArray(user.mup[0],&x_ptr);
315:   VecRestoreArray(user.lambda[0],&y_ptr);
316:   PetscPrintf(PETSC_COMM_WORLD,"\n sensitivity wrt parameters: d[y(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp));

318:   VecGetArray(user.mup[1],&x_ptr);
319:   VecGetArray(user.lambda[1],&y_ptr);
320:   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];
321:   VecRestoreArray(user.mup[1],&x_ptr);
322:   VecRestoreArray(user.lambda[1],&y_ptr);
323:   PetscPrintf(PETSC_COMM_WORLD,"\n sensivitity wrt parameters: d[z(tf)]/d[mu]\n%g\n",(double)PetscRealPart(derp));

325:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
326:      Free work space.  All PETSc objects should be destroyed when they
327:      are no longer needed.
328:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
329:   MatDestroy(&user.A);
330:   MatDestroy(&user.Jacp);
331:   VecDestroy(&user.U);
332:   VecDestroy(&user.lambda[0]);
333:   VecDestroy(&user.lambda[1]);
334:   VecDestroy(&user.mup[0]);
335:   VecDestroy(&user.mup[1]);
336:   TSDestroy(&ts);

338:   PetscFinalize();
339:   return(ierr);
340: }

342: /*TEST

344:     test:
345:       requires: revolve
346:       args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -viewer_binary_skip_info -ts_dt 0.001 -mu 100000

348:     test:
349:       suffix: 2
350:       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only

352:     test:
353:       suffix: 3
354:       args: -ts_type cn -ts_dt 0.001 -mu 100000 -ts_max_steps 15 -ts_trajectory_type memory -ts_trajectory_solution_only 0

357:     test:
358:       suffix: 4
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 -ts_trajectory_save_stack

362:     test:
363:       suffix: 5
364:       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

367:     test:
368:       suffix: 6
369:       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

372:     test:
373:       suffix: 7
374:       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

377:     test:
378:       suffix: 8
379:       requires: revolve
380:       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

383:     test:
384:       suffix: 9
385:       requires: revolve
386:       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

389:     test:
390:       requires: revolve
391:       suffix: 10
392:       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

395:     test:
396:       requires: revolve
397:       suffix: 11
398:       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

401:     test:
402:       suffix: 12
403:       requires: revolve
404:       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

407:     test:
408:       suffix: 13
409:       requires: revolve
410:       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

413:     test:
414:       suffix: 14
415:       requires: revolve
416:       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

419:     test:
420:       suffix: 15
421:       requires: revolve
422:       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

425:     test:
426:       suffix: 16
427:       requires: revolve
428:       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

431:     test:
432:       suffix: 17
433:       requires: revolve
434:       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

437:     test:
438:       suffix: 18
439:       requires: revolve
440:       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

443:     test:
444:       suffix: 19
445:       requires: revolve
446:       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

449:     test:
450:       suffix: 20
451:       requires: revolve
452:       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