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*/