Actual source code: ex20opt_p.c


  2: static char help[] = "Solves the van der Pol equation.\n\
  3: Input parameters include:\n";

  5: /* ------------------------------------------------------------------------

  7:   Notes:
  8:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
  9:   The nonlinear problem is written in a DAE equivalent form.
 10:   The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
 11:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
 12:   ------------------------------------------------------------------------- */
 13: #include <petsctao.h>
 14: #include <petscts.h>

 16: typedef struct _n_User *User;
 17: struct _n_User {
 18:   TS        ts;
 19:   PetscReal mu;
 20:   PetscReal next_output;

 22:   /* Sensitivity analysis support */
 23:   PetscReal ftime;
 24:   Mat       A;                    /* Jacobian matrix */
 25:   Mat       Jacp;                 /* JacobianP matrix */
 26:   Mat       H;                    /* Hessian matrix for optimization */
 27:   Vec       U, Lambda[1], Mup[1]; /* adjoint variables */
 28:   Vec       Lambda2[1], Mup2[1];  /* second-order adjoint variables */
 29:   Vec       Ihp1[1];              /* working space for Hessian evaluations */
 30:   Vec       Ihp2[1];              /* working space for Hessian evaluations */
 31:   Vec       Ihp3[1];              /* working space for Hessian evaluations */
 32:   Vec       Ihp4[1];              /* working space for Hessian evaluations */
 33:   Vec       Dir;                  /* direction vector */
 34:   PetscReal ob[2];                /* observation used by the cost function */
 35:   PetscBool implicitform;         /* implicit ODE? */
 36: };

 38: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 39: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 40: PetscErrorCode Adjoint2(Vec, PetscScalar[], User);

 42: /* ----------------------- Explicit form of the ODE  -------------------- */

 44: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
 45: {
 46:   User               user = (User)ctx;
 47:   PetscScalar       *f;
 48:   const PetscScalar *u;

 51:   VecGetArrayRead(U, &u);
 52:   VecGetArray(F, &f);
 53:   f[0] = u[1];
 54:   f[1] = user->mu * ((1. - u[0] * u[0]) * u[1] - u[0]);
 55:   VecRestoreArrayRead(U, &u);
 56:   VecRestoreArray(F, &f);
 57:   return 0;
 58: }

 60: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
 61: {
 62:   User               user     = (User)ctx;
 63:   PetscReal          mu       = user->mu;
 64:   PetscInt           rowcol[] = {0, 1};
 65:   PetscScalar        J[2][2];
 66:   const PetscScalar *u;

 69:   VecGetArrayRead(U, &u);

 71:   J[0][0] = 0;
 72:   J[1][0] = -mu * (2.0 * u[1] * u[0] + 1.);
 73:   J[0][1] = 1.0;
 74:   J[1][1] = mu * (1.0 - u[0] * u[0]);
 75:   MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
 76:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
 77:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 78:   if (B && A != B) {
 79:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
 80:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
 81:   }

 83:   VecRestoreArrayRead(U, &u);
 84:   return 0;
 85: }

 87: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 88: {
 89:   const PetscScalar *vl, *vr, *u;
 90:   PetscScalar       *vhv;
 91:   PetscScalar        dJdU[2][2][2] = {{{0}}};
 92:   PetscInt           i, j, k;
 93:   User               user = (User)ctx;

 96:   VecGetArrayRead(U, &u);
 97:   VecGetArrayRead(Vl[0], &vl);
 98:   VecGetArrayRead(Vr, &vr);
 99:   VecGetArray(VHV[0], &vhv);

101:   dJdU[1][0][0] = -2. * user->mu * u[1];
102:   dJdU[1][1][0] = -2. * user->mu * u[0];
103:   dJdU[1][0][1] = -2. * user->mu * u[0];
104:   for (j = 0; j < 2; j++) {
105:     vhv[j] = 0;
106:     for (k = 0; k < 2; k++)
107:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
108:   }

110:   VecRestoreArrayRead(U, &u);
111:   VecRestoreArrayRead(Vl[0], &vl);
112:   VecRestoreArrayRead(Vr, &vr);
113:   VecRestoreArray(VHV[0], &vhv);
114:   return 0;
115: }

117: static PetscErrorCode RHSHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
118: {
119:   const PetscScalar *vl, *vr, *u;
120:   PetscScalar       *vhv;
121:   PetscScalar        dJdP[2][2][1] = {{{0}}};
122:   PetscInt           i, j, k;

125:   VecGetArrayRead(U, &u);
126:   VecGetArrayRead(Vl[0], &vl);
127:   VecGetArrayRead(Vr, &vr);
128:   VecGetArray(VHV[0], &vhv);

130:   dJdP[1][0][0] = -(1. + 2. * u[0] * u[1]);
131:   dJdP[1][1][0] = 1. - u[0] * u[0];
132:   for (j = 0; j < 2; j++) {
133:     vhv[j] = 0;
134:     for (k = 0; k < 1; k++)
135:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
136:   }

138:   VecRestoreArrayRead(U, &u);
139:   VecRestoreArrayRead(Vl[0], &vl);
140:   VecRestoreArrayRead(Vr, &vr);
141:   VecRestoreArray(VHV[0], &vhv);
142:   return 0;
143: }

145: static PetscErrorCode RHSHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
146: {
147:   const PetscScalar *vl, *vr, *u;
148:   PetscScalar       *vhv;
149:   PetscScalar        dJdU[2][1][2] = {{{0}}};
150:   PetscInt           i, j, k;

153:   VecGetArrayRead(U, &u);
154:   VecGetArrayRead(Vl[0], &vl);
155:   VecGetArrayRead(Vr, &vr);
156:   VecGetArray(VHV[0], &vhv);

158:   dJdU[1][0][0] = -1. - 2. * u[1] * u[0];
159:   dJdU[1][0][1] = 1. - u[0] * u[0];
160:   for (j = 0; j < 1; j++) {
161:     vhv[j] = 0;
162:     for (k = 0; k < 2; k++)
163:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
164:   }

166:   VecRestoreArrayRead(U, &u);
167:   VecRestoreArrayRead(Vl[0], &vl);
168:   VecRestoreArrayRead(Vr, &vr);
169:   VecRestoreArray(VHV[0], &vhv);
170:   return 0;
171: }

173: static PetscErrorCode RHSHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
174: {
176:   return 0;
177: }

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

181: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
182: {
183:   User               user = (User)ctx;
184:   PetscScalar       *f;
185:   const PetscScalar *u, *udot;

188:   VecGetArrayRead(U, &u);
189:   VecGetArrayRead(Udot, &udot);
190:   VecGetArray(F, &f);

192:   f[0] = udot[0] - u[1];
193:   f[1] = udot[1] - user->mu * ((1.0 - u[0] * u[0]) * u[1] - u[0]);

195:   VecRestoreArrayRead(U, &u);
196:   VecRestoreArrayRead(Udot, &udot);
197:   VecRestoreArray(F, &f);
198:   return 0;
199: }

201: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
202: {
203:   User               user     = (User)ctx;
204:   PetscInt           rowcol[] = {0, 1};
205:   PetscScalar        J[2][2];
206:   const PetscScalar *u;

209:   VecGetArrayRead(U, &u);

211:   J[0][0] = a;
212:   J[0][1] = -1.0;
213:   J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
214:   J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
215:   MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES);
216:   VecRestoreArrayRead(U, &u);
217:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
218:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
219:   if (A != B) {
220:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
221:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
222:   }

224:   VecRestoreArrayRead(U, &u);
225:   return 0;
226: }

228: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec U, Mat A, void *ctx)
229: {
230:   PetscInt           row[] = {0, 1}, col[] = {0};
231:   PetscScalar        J[2][1];
232:   const PetscScalar *u;

235:   VecGetArrayRead(U, &u);

237:   J[0][0] = 0;
238:   J[1][0] = (1. - u[0] * u[0]) * u[1] - u[0];
239:   MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES);
240:   VecRestoreArrayRead(U, &u);
241:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
242:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

244:   VecRestoreArrayRead(U, &u);
245:   return 0;
246: }

248: static PetscErrorCode IHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
249: {
250:   const PetscScalar *vl, *vr, *u;
251:   PetscScalar       *vhv;
252:   PetscScalar        dJdU[2][2][2] = {{{0}}};
253:   PetscInt           i, j, k;
254:   User               user = (User)ctx;

257:   VecGetArrayRead(U, &u);
258:   VecGetArrayRead(Vl[0], &vl);
259:   VecGetArrayRead(Vr, &vr);
260:   VecGetArray(VHV[0], &vhv);

262:   dJdU[1][0][0] = 2. * user->mu * u[1];
263:   dJdU[1][1][0] = 2. * user->mu * u[0];
264:   dJdU[1][0][1] = 2. * user->mu * u[0];
265:   for (j = 0; j < 2; j++) {
266:     vhv[j] = 0;
267:     for (k = 0; k < 2; k++)
268:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
269:   }

271:   VecRestoreArrayRead(U, &u);
272:   VecRestoreArrayRead(Vl[0], &vl);
273:   VecRestoreArrayRead(Vr, &vr);
274:   VecRestoreArray(VHV[0], &vhv);
275:   return 0;
276: }

278: static PetscErrorCode IHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
279: {
280:   const PetscScalar *vl, *vr, *u;
281:   PetscScalar       *vhv;
282:   PetscScalar        dJdP[2][2][1] = {{{0}}};
283:   PetscInt           i, j, k;

286:   VecGetArrayRead(U, &u);
287:   VecGetArrayRead(Vl[0], &vl);
288:   VecGetArrayRead(Vr, &vr);
289:   VecGetArray(VHV[0], &vhv);

291:   dJdP[1][0][0] = 1. + 2. * u[0] * u[1];
292:   dJdP[1][1][0] = u[0] * u[0] - 1.;
293:   for (j = 0; j < 2; j++) {
294:     vhv[j] = 0;
295:     for (k = 0; k < 1; k++)
296:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
297:   }

299:   VecRestoreArrayRead(U, &u);
300:   VecRestoreArrayRead(Vl[0], &vl);
301:   VecRestoreArrayRead(Vr, &vr);
302:   VecRestoreArray(VHV[0], &vhv);
303:   return 0;
304: }

306: static PetscErrorCode IHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
307: {
308:   const PetscScalar *vl, *vr, *u;
309:   PetscScalar       *vhv;
310:   PetscScalar        dJdU[2][1][2] = {{{0}}};
311:   PetscInt           i, j, k;

314:   VecGetArrayRead(U, &u);
315:   VecGetArrayRead(Vl[0], &vl);
316:   VecGetArrayRead(Vr, &vr);
317:   VecGetArray(VHV[0], &vhv);

319:   dJdU[1][0][0] = 1. + 2. * u[1] * u[0];
320:   dJdU[1][0][1] = u[0] * u[0] - 1.;
321:   for (j = 0; j < 1; j++) {
322:     vhv[j] = 0;
323:     for (k = 0; k < 2; k++)
324:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
325:   }

327:   VecRestoreArrayRead(U, &u);
328:   VecRestoreArrayRead(Vl[0], &vl);
329:   VecRestoreArrayRead(Vr, &vr);
330:   VecRestoreArray(VHV[0], &vhv);
331:   return 0;
332: }

334: static PetscErrorCode IHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
335: {
337:   return 0;
338: }

340: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
341: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
342: {
343:   const PetscScalar *x;
344:   PetscReal          tfinal, dt;
345:   User               user = (User)ctx;
346:   Vec                interpolatedX;

349:   TSGetTimeStep(ts, &dt);
350:   TSGetMaxTime(ts, &tfinal);

352:   while (user->next_output <= t && user->next_output <= tfinal) {
353:     VecDuplicate(X, &interpolatedX);
354:     TSInterpolate(ts, user->next_output, interpolatedX);
355:     VecGetArrayRead(interpolatedX, &x);
356:     PetscPrintf(PETSC_COMM_WORLD, "[%g] %" PetscInt_FMT " TS %g (dt = %g) X %g %g\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1]));
357:     VecRestoreArrayRead(interpolatedX, &x);
358:     VecDestroy(&interpolatedX);
359:     user->next_output += PetscRealConstant(0.1);
360:   }
361:   return 0;
362: }

364: int main(int argc, char **argv)
365: {
366:   Vec                P;
367:   PetscBool          monitor = PETSC_FALSE;
368:   PetscScalar       *x_ptr;
369:   const PetscScalar *y_ptr;
370:   PetscMPIInt        size;
371:   struct _n_User     user;
372:   Tao                tao;
373:   KSP                ksp;
374:   PC                 pc;

376:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377:      Initialize program
378:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
380:   PetscInitialize(&argc, &argv, NULL, help);
381:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

384:   /* Create TAO solver and set desired solution method */
385:   TaoCreate(PETSC_COMM_WORLD, &tao);
386:   TaoSetType(tao, TAOBQNLS);

388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389:     Set runtime options
390:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391:   user.next_output  = 0.0;
392:   user.mu           = PetscRealConstant(1.0e3);
393:   user.ftime        = PetscRealConstant(0.5);
394:   user.implicitform = PETSC_TRUE;

396:   PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL);
397:   PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL);
398:   PetscOptionsGetBool(NULL, NULL, "-implicitform", &user.implicitform, NULL);

400:   /* Create necessary matrix and vectors, solve same ODE on every process */
401:   MatCreate(PETSC_COMM_WORLD, &user.A);
402:   MatSetSizes(user.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2);
403:   MatSetFromOptions(user.A);
404:   MatSetUp(user.A);
405:   MatCreateVecs(user.A, &user.U, NULL);
406:   MatCreateVecs(user.A, &user.Lambda[0], NULL);
407:   MatCreateVecs(user.A, &user.Lambda2[0], NULL);
408:   MatCreateVecs(user.A, &user.Ihp1[0], NULL);
409:   MatCreateVecs(user.A, &user.Ihp2[0], NULL);

411:   MatCreate(PETSC_COMM_WORLD, &user.Jacp);
412:   MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1);
413:   MatSetFromOptions(user.Jacp);
414:   MatSetUp(user.Jacp);
415:   MatCreateVecs(user.Jacp, &user.Dir, NULL);
416:   MatCreateVecs(user.Jacp, &user.Mup[0], NULL);
417:   MatCreateVecs(user.Jacp, &user.Mup2[0], NULL);
418:   MatCreateVecs(user.Jacp, &user.Ihp3[0], NULL);
419:   MatCreateVecs(user.Jacp, &user.Ihp4[0], NULL);

421:   /* Create timestepping solver context */
422:   TSCreate(PETSC_COMM_WORLD, &user.ts);
423:   TSSetEquationType(user.ts, TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
424:   if (user.implicitform) {
425:     TSSetIFunction(user.ts, NULL, IFunction, &user);
426:     TSSetIJacobian(user.ts, user.A, user.A, IJacobian, &user);
427:     TSSetType(user.ts, TSCN);
428:   } else {
429:     TSSetRHSFunction(user.ts, NULL, RHSFunction, &user);
430:     TSSetRHSJacobian(user.ts, user.A, user.A, RHSJacobian, &user);
431:     TSSetType(user.ts, TSRK);
432:   }
433:   TSSetRHSJacobianP(user.ts, user.Jacp, RHSJacobianP, &user);
434:   TSSetMaxTime(user.ts, user.ftime);
435:   TSSetExactFinalTime(user.ts, TS_EXACTFINALTIME_MATCHSTEP);

437:   if (monitor) TSMonitorSet(user.ts, Monitor, &user, NULL);

439:   /* Set ODE initial conditions */
440:   VecGetArray(user.U, &x_ptr);
441:   x_ptr[0] = 2.0;
442:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
443:   VecRestoreArray(user.U, &x_ptr);
444:   TSSetTimeStep(user.ts, PetscRealConstant(0.001));

446:   /* Set runtime options */
447:   TSSetFromOptions(user.ts);

449:   TSSolve(user.ts, user.U);
450:   VecGetArrayRead(user.U, &y_ptr);
451:   user.ob[0] = y_ptr[0];
452:   user.ob[1] = y_ptr[1];
453:   VecRestoreArrayRead(user.U, &y_ptr);

455:   /* Save trajectory of solution so that TSAdjointSolve() may be used.
456:      Skip checkpointing for the first TSSolve since no adjoint run follows it.
457:    */
458:   TSSetSaveTrajectory(user.ts);

460:   /* Optimization starts */
461:   MatCreate(PETSC_COMM_WORLD, &user.H);
462:   MatSetSizes(user.H, PETSC_DECIDE, PETSC_DECIDE, 1, 1);
463:   MatSetUp(user.H); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */

465:   /* Set initial solution guess */
466:   MatCreateVecs(user.Jacp, &P, NULL);
467:   VecGetArray(P, &x_ptr);
468:   x_ptr[0] = PetscRealConstant(1.2);
469:   VecRestoreArray(P, &x_ptr);
470:   TaoSetSolution(tao, P);

472:   /* Set routine for function and gradient evaluation */
473:   TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user);
474:   TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user);

476:   /* Check for any TAO command line options */
477:   TaoGetKSP(tao, &ksp);
478:   if (ksp) {
479:     KSPGetPC(ksp, &pc);
480:     PCSetType(pc, PCNONE);
481:   }
482:   TaoSetFromOptions(tao);

484:   TaoSolve(tao);

486:   VecView(P, PETSC_VIEWER_STDOUT_WORLD);
487:   /* Free TAO data structures */
488:   TaoDestroy(&tao);

490:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
491:      Free work space.  All PETSc objects should be destroyed when they
492:      are no longer needed.
493:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
494:   MatDestroy(&user.H);
495:   MatDestroy(&user.A);
496:   VecDestroy(&user.U);
497:   MatDestroy(&user.Jacp);
498:   VecDestroy(&user.Lambda[0]);
499:   VecDestroy(&user.Mup[0]);
500:   VecDestroy(&user.Lambda2[0]);
501:   VecDestroy(&user.Mup2[0]);
502:   VecDestroy(&user.Ihp1[0]);
503:   VecDestroy(&user.Ihp2[0]);
504:   VecDestroy(&user.Ihp3[0]);
505:   VecDestroy(&user.Ihp4[0]);
506:   VecDestroy(&user.Dir);
507:   TSDestroy(&user.ts);
508:   VecDestroy(&P);
509:   PetscFinalize();
510:   return 0;
511: }

513: /* ------------------------------------------------------------------ */
514: /*
515:    FormFunctionGradient - Evaluates the function and corresponding gradient.

517:    Input Parameters:
518:    tao - the Tao context
519:    X   - the input vector
520:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

522:    Output Parameters:
523:    f   - the newly evaluated function
524:    G   - the newly evaluated gradient
525: */
526: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx)
527: {
528:   User               user_ptr = (User)ctx;
529:   TS                 ts       = user_ptr->ts;
530:   PetscScalar       *x_ptr, *g;
531:   const PetscScalar *y_ptr;

534:   VecGetArrayRead(P, &y_ptr);
535:   user_ptr->mu = y_ptr[0];
536:   VecRestoreArrayRead(P, &y_ptr);

538:   TSSetTime(ts, 0.0);
539:   TSSetStepNumber(ts, 0);
540:   TSSetTimeStep(ts, PetscRealConstant(0.001)); /* can be overwritten by command line options */
541:   TSSetFromOptions(ts);
542:   VecGetArray(user_ptr->U, &x_ptr);
543:   x_ptr[0] = 2.0;
544:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user_ptr->mu) - 292.0 / (2187.0 * user_ptr->mu * user_ptr->mu);
545:   VecRestoreArray(user_ptr->U, &x_ptr);

547:   TSSolve(ts, user_ptr->U);

549:   VecGetArrayRead(user_ptr->U, &y_ptr);
550:   *f = (y_ptr[0] - user_ptr->ob[0]) * (y_ptr[0] - user_ptr->ob[0]) + (y_ptr[1] - user_ptr->ob[1]) * (y_ptr[1] - user_ptr->ob[1]);

552:   /*   Reset initial conditions for the adjoint integration */
553:   VecGetArray(user_ptr->Lambda[0], &x_ptr);
554:   x_ptr[0] = 2. * (y_ptr[0] - user_ptr->ob[0]);
555:   x_ptr[1] = 2. * (y_ptr[1] - user_ptr->ob[1]);
556:   VecRestoreArrayRead(user_ptr->U, &y_ptr);
557:   VecRestoreArray(user_ptr->Lambda[0], &x_ptr);

559:   VecGetArray(user_ptr->Mup[0], &x_ptr);
560:   x_ptr[0] = 0.0;
561:   VecRestoreArray(user_ptr->Mup[0], &x_ptr);
562:   TSSetCostGradients(ts, 1, user_ptr->Lambda, user_ptr->Mup);

564:   TSAdjointSolve(ts);

566:   VecGetArray(user_ptr->Mup[0], &x_ptr);
567:   VecGetArrayRead(user_ptr->Lambda[0], &y_ptr);
568:   VecGetArray(G, &g);
569:   g[0] = y_ptr[1] * (-10.0 / (81.0 * user_ptr->mu * user_ptr->mu) + 2.0 * 292.0 / (2187.0 * user_ptr->mu * user_ptr->mu * user_ptr->mu)) + x_ptr[0];
570:   VecRestoreArray(user_ptr->Mup[0], &x_ptr);
571:   VecRestoreArrayRead(user_ptr->Lambda[0], &y_ptr);
572:   VecRestoreArray(G, &g);
573:   return 0;
574: }

576: PetscErrorCode FormHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
577: {
578:   User           user_ptr = (User)ctx;
579:   PetscScalar    harr[1];
580:   const PetscInt rows[1] = {0};
581:   PetscInt       col     = 0;

584:   Adjoint2(P, harr, user_ptr);
585:   MatSetValues(H, 1, rows, 1, &col, harr, INSERT_VALUES);

587:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
588:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
589:   if (H != Hpre) {
590:     MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY);
591:     MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY);
592:   }
593:   return 0;
594: }

596: PetscErrorCode Adjoint2(Vec P, PetscScalar arr[], User ctx)
597: {
598:   TS                 ts = ctx->ts;
599:   const PetscScalar *z_ptr;
600:   PetscScalar       *x_ptr, *y_ptr, dzdp, dzdp2;
601:   Mat                tlmsen;

604:   /* Reset TSAdjoint so that AdjointSetUp will be called again */
605:   TSAdjointReset(ts);

607:   /* The directional vector should be 1 since it is one-dimensional */
608:   VecGetArray(ctx->Dir, &x_ptr);
609:   x_ptr[0] = 1.;
610:   VecRestoreArray(ctx->Dir, &x_ptr);

612:   VecGetArrayRead(P, &z_ptr);
613:   ctx->mu = z_ptr[0];
614:   VecRestoreArrayRead(P, &z_ptr);

616:   dzdp  = -10.0 / (81.0 * ctx->mu * ctx->mu) + 2.0 * 292.0 / (2187.0 * ctx->mu * ctx->mu * ctx->mu);
617:   dzdp2 = 2. * 10.0 / (81.0 * ctx->mu * ctx->mu * ctx->mu) - 3.0 * 2.0 * 292.0 / (2187.0 * ctx->mu * ctx->mu * ctx->mu * ctx->mu);

619:   TSSetTime(ts, 0.0);
620:   TSSetStepNumber(ts, 0);
621:   TSSetTimeStep(ts, PetscRealConstant(0.001)); /* can be overwritten by command line options */
622:   TSSetFromOptions(ts);
623:   TSSetCostHessianProducts(ts, 1, ctx->Lambda2, ctx->Mup2, ctx->Dir);

625:   MatZeroEntries(ctx->Jacp);
626:   MatSetValue(ctx->Jacp, 1, 0, dzdp, INSERT_VALUES);
627:   MatAssemblyBegin(ctx->Jacp, MAT_FINAL_ASSEMBLY);
628:   MatAssemblyEnd(ctx->Jacp, MAT_FINAL_ASSEMBLY);

630:   TSAdjointSetForward(ts, ctx->Jacp);
631:   VecGetArray(ctx->U, &y_ptr);
632:   y_ptr[0] = 2.0;
633:   y_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * ctx->mu) - 292.0 / (2187.0 * ctx->mu * ctx->mu);
634:   VecRestoreArray(ctx->U, &y_ptr);
635:   TSSolve(ts, ctx->U);

637:   /* Set terminal conditions for first- and second-order adjonts */
638:   VecGetArrayRead(ctx->U, &z_ptr);
639:   VecGetArray(ctx->Lambda[0], &y_ptr);
640:   y_ptr[0] = 2. * (z_ptr[0] - ctx->ob[0]);
641:   y_ptr[1] = 2. * (z_ptr[1] - ctx->ob[1]);
642:   VecRestoreArray(ctx->Lambda[0], &y_ptr);
643:   VecRestoreArrayRead(ctx->U, &z_ptr);
644:   VecGetArray(ctx->Mup[0], &y_ptr);
645:   y_ptr[0] = 0.0;
646:   VecRestoreArray(ctx->Mup[0], &y_ptr);
647:   TSForwardGetSensitivities(ts, NULL, &tlmsen);
648:   MatDenseGetColumn(tlmsen, 0, &x_ptr);
649:   VecGetArray(ctx->Lambda2[0], &y_ptr);
650:   y_ptr[0] = 2. * x_ptr[0];
651:   y_ptr[1] = 2. * x_ptr[1];
652:   VecRestoreArray(ctx->Lambda2[0], &y_ptr);
653:   VecGetArray(ctx->Mup2[0], &y_ptr);
654:   y_ptr[0] = 0.0;
655:   VecRestoreArray(ctx->Mup2[0], &y_ptr);
656:   MatDenseRestoreColumn(tlmsen, &x_ptr);
657:   TSSetCostGradients(ts, 1, ctx->Lambda, ctx->Mup);
658:   if (ctx->implicitform) {
659:     TSSetIHessianProduct(ts, ctx->Ihp1, IHessianProductUU, ctx->Ihp2, IHessianProductUP, ctx->Ihp3, IHessianProductPU, ctx->Ihp4, IHessianProductPP, ctx);
660:   } else {
661:     TSSetRHSHessianProduct(ts, ctx->Ihp1, RHSHessianProductUU, ctx->Ihp2, RHSHessianProductUP, ctx->Ihp3, RHSHessianProductPU, ctx->Ihp4, RHSHessianProductPP, ctx);
662:   }
663:   TSAdjointSolve(ts);

665:   VecGetArray(ctx->Lambda[0], &x_ptr);
666:   VecGetArray(ctx->Lambda2[0], &y_ptr);
667:   VecGetArrayRead(ctx->Mup2[0], &z_ptr);

669:   arr[0] = x_ptr[1] * dzdp2 + y_ptr[1] * dzdp2 + z_ptr[0];

671:   VecRestoreArray(ctx->Lambda2[0], &x_ptr);
672:   VecRestoreArray(ctx->Lambda2[0], &y_ptr);
673:   VecRestoreArrayRead(ctx->Mup2[0], &z_ptr);

675:   /* Disable second-order adjoint mode */
676:   TSAdjointReset(ts);
677:   TSAdjointResetForward(ts);
678:   return 0;
679: }

681: /*TEST
682:     build:
683:       requires: !complex !single
684:     test:
685:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
686:       output_file: output/ex20opt_p_1.out

688:     test:
689:       suffix: 2
690:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none
691:       output_file: output/ex20opt_p_2.out

693:     test:
694:       suffix: 3
695:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
696:       output_file: output/ex20opt_p_3.out

698:     test:
699:       suffix: 4
700:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none
701:       output_file: output/ex20opt_p_4.out

703: TEST*/