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++)
108: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
109: }
111: VecRestoreArrayRead(U,&u);
112: VecRestoreArrayRead(Vl[0],&vl);
113: VecRestoreArrayRead(Vr,&vr);
114: VecRestoreArray(VHV[0],&vhv);
115: return 0;
116: }
118: static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
119: {
120: const PetscScalar *vl,*vr,*u;
121: PetscScalar *vhv;
122: PetscScalar dJdP[2][2][1]={{{0}}};
123: PetscInt i,j,k;
126: VecGetArrayRead(U,&u);
127: VecGetArrayRead(Vl[0],&vl);
128: VecGetArrayRead(Vr,&vr);
129: VecGetArray(VHV[0],&vhv);
131: dJdP[1][0][0] = -(1.+2.*u[0]*u[1]);
132: dJdP[1][1][0] = 1.-u[0]*u[0];
133: for (j=0; j<2; j++) {
134: vhv[j] = 0;
135: for (k=0; k<1; k++)
136: for (i=0; i<2; i++)
137: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
138: }
140: VecRestoreArrayRead(U,&u);
141: VecRestoreArrayRead(Vl[0],&vl);
142: VecRestoreArrayRead(Vr,&vr);
143: VecRestoreArray(VHV[0],&vhv);
144: return 0;
145: }
147: static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
148: {
149: const PetscScalar *vl,*vr,*u;
150: PetscScalar *vhv;
151: PetscScalar dJdU[2][1][2]={{{0}}};
152: PetscInt i,j,k;
155: VecGetArrayRead(U,&u);
156: VecGetArrayRead(Vl[0],&vl);
157: VecGetArrayRead(Vr,&vr);
158: VecGetArray(VHV[0],&vhv);
160: dJdU[1][0][0] = -1.-2.*u[1]*u[0];
161: dJdU[1][0][1] = 1.-u[0]*u[0];
162: for (j=0; j<1; j++) {
163: vhv[j] = 0;
164: for (k=0; k<2; k++)
165: for (i=0; i<2; i++)
166: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
167: }
169: VecRestoreArrayRead(U,&u);
170: VecRestoreArrayRead(Vl[0],&vl);
171: VecRestoreArrayRead(Vr,&vr);
172: VecRestoreArray(VHV[0],&vhv);
173: return 0;
174: }
176: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
177: {
179: return 0;
180: }
182: /* ----------------------- Implicit form of the ODE -------------------- */
184: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
185: {
186: User user = (User)ctx;
187: PetscScalar *f;
188: const PetscScalar *u,*udot;
191: VecGetArrayRead(U,&u);
192: VecGetArrayRead(Udot,&udot);
193: VecGetArray(F,&f);
195: f[0] = udot[0] - u[1];
196: f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
198: VecRestoreArrayRead(U,&u);
199: VecRestoreArrayRead(Udot,&udot);
200: VecRestoreArray(F,&f);
201: return 0;
202: }
204: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
205: {
206: User user = (User)ctx;
207: PetscInt rowcol[] = {0,1};
208: PetscScalar J[2][2];
209: const PetscScalar *u;
212: VecGetArrayRead(U,&u);
214: J[0][0] = a; J[0][1] = -1.0;
215: J[1][0] = user->mu*(1.0 + 2.0*u[0]*u[1]); J[1][1] = a - user->mu*(1.0-u[0]*u[0]);
216: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
217: VecRestoreArrayRead(U,&u);
218: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
219: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
220: if (A != B) {
221: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
222: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
223: }
225: VecRestoreArrayRead(U,&u);
226: return 0;
227: }
229: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
230: {
231: PetscInt row[] = {0,1},col[]={0};
232: PetscScalar J[2][1];
233: const PetscScalar *u;
236: VecGetArrayRead(U,&u);
238: J[0][0] = 0;
239: J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
240: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
241: VecRestoreArrayRead(U,&u);
242: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
243: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
245: VecRestoreArrayRead(U,&u);
246: return 0;
247: }
249: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
250: {
251: const PetscScalar *vl,*vr,*u;
252: PetscScalar *vhv;
253: PetscScalar dJdU[2][2][2]={{{0}}};
254: PetscInt i,j,k;
255: User user = (User)ctx;
258: VecGetArrayRead(U,&u);
259: VecGetArrayRead(Vl[0],&vl);
260: VecGetArrayRead(Vr,&vr);
261: VecGetArray(VHV[0],&vhv);
263: dJdU[1][0][0] = 2.*user->mu*u[1];
264: dJdU[1][1][0] = 2.*user->mu*u[0];
265: dJdU[1][0][1] = 2.*user->mu*u[0];
266: for (j=0; j<2; j++) {
267: vhv[j] = 0;
268: for (k=0; k<2; k++)
269: for (i=0; i<2; i++)
270: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
271: }
273: VecRestoreArrayRead(U,&u);
274: VecRestoreArrayRead(Vl[0],&vl);
275: VecRestoreArrayRead(Vr,&vr);
276: VecRestoreArray(VHV[0],&vhv);
277: return 0;
278: }
280: static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
281: {
282: const PetscScalar *vl,*vr,*u;
283: PetscScalar *vhv;
284: PetscScalar dJdP[2][2][1]={{{0}}};
285: PetscInt i,j,k;
288: VecGetArrayRead(U,&u);
289: VecGetArrayRead(Vl[0],&vl);
290: VecGetArrayRead(Vr,&vr);
291: VecGetArray(VHV[0],&vhv);
293: dJdP[1][0][0] = 1.+2.*u[0]*u[1];
294: dJdP[1][1][0] = u[0]*u[0]-1.;
295: for (j=0; j<2; j++) {
296: vhv[j] = 0;
297: for (k=0; k<1; k++)
298: for (i=0; i<2; i++)
299: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
300: }
302: VecRestoreArrayRead(U,&u);
303: VecRestoreArrayRead(Vl[0],&vl);
304: VecRestoreArrayRead(Vr,&vr);
305: VecRestoreArray(VHV[0],&vhv);
306: return 0;
307: }
309: static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
310: {
311: const PetscScalar *vl,*vr,*u;
312: PetscScalar *vhv;
313: PetscScalar dJdU[2][1][2]={{{0}}};
314: PetscInt i,j,k;
317: VecGetArrayRead(U,&u);
318: VecGetArrayRead(Vl[0],&vl);
319: VecGetArrayRead(Vr,&vr);
320: VecGetArray(VHV[0],&vhv);
322: dJdU[1][0][0] = 1.+2.*u[1]*u[0];
323: dJdU[1][0][1] = u[0]*u[0]-1.;
324: for (j=0; j<1; j++) {
325: vhv[j] = 0;
326: for (k=0; k<2; k++)
327: for (i=0; i<2; i++)
328: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
329: }
331: VecRestoreArrayRead(U,&u);
332: VecRestoreArrayRead(Vl[0],&vl);
333: VecRestoreArrayRead(Vr,&vr);
334: VecRestoreArray(VHV[0],&vhv);
335: return 0;
336: }
338: static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
339: {
341: return 0;
342: }
344: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
345: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
346: {
347: PetscErrorCode ierr;
348: const PetscScalar *x;
349: PetscReal tfinal, dt;
350: User user = (User)ctx;
351: Vec interpolatedX;
354: TSGetTimeStep(ts,&dt);
355: TSGetMaxTime(ts,&tfinal);
357: while (user->next_output <= t && user->next_output <= tfinal) {
358: VecDuplicate(X,&interpolatedX);
359: TSInterpolate(ts,user->next_output,interpolatedX);
360: VecGetArrayRead(interpolatedX,&x);
361: PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
362: (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
363: (double)PetscRealPart(x[1]));
364: VecRestoreArrayRead(interpolatedX,&x);
365: VecDestroy(&interpolatedX);
366: user->next_output += PetscRealConstant(0.1);
367: }
368: return 0;
369: }
371: int main(int argc,char **argv)
372: {
373: Vec P;
374: PetscBool monitor = PETSC_FALSE;
375: PetscScalar *x_ptr;
376: const PetscScalar *y_ptr;
377: PetscMPIInt size;
378: struct _n_User user;
379: Tao tao;
380: KSP ksp;
381: PC pc;
383: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
384: Initialize program
385: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
386: PetscInitialize(&argc,&argv,NULL,help);
387: MPI_Comm_size(PETSC_COMM_WORLD,&size);
390: /* Create TAO solver and set desired solution method */
391: TaoCreate(PETSC_COMM_WORLD,&tao);
392: TaoSetType(tao,TAOBQNLS);
394: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
395: Set runtime options
396: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
397: user.next_output = 0.0;
398: user.mu = PetscRealConstant(1.0e3);
399: user.ftime = PetscRealConstant(0.5);
400: user.implicitform = PETSC_TRUE;
402: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
403: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
404: PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);
406: /* Create necessary matrix and vectors, solve same ODE on every process */
407: MatCreate(PETSC_COMM_WORLD,&user.A);
408: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
409: MatSetFromOptions(user.A);
410: MatSetUp(user.A);
411: MatCreateVecs(user.A,&user.U,NULL);
412: MatCreateVecs(user.A,&user.Lambda[0],NULL);
413: MatCreateVecs(user.A,&user.Lambda2[0],NULL);
414: MatCreateVecs(user.A,&user.Ihp1[0],NULL);
415: MatCreateVecs(user.A,&user.Ihp2[0],NULL);
417: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
418: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
419: MatSetFromOptions(user.Jacp);
420: MatSetUp(user.Jacp);
421: MatCreateVecs(user.Jacp,&user.Dir,NULL);
422: MatCreateVecs(user.Jacp,&user.Mup[0],NULL);
423: MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);
424: MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);
425: MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);
427: /* Create timestepping solver context */
428: TSCreate(PETSC_COMM_WORLD,&user.ts);
429: TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
430: if (user.implicitform) {
431: TSSetIFunction(user.ts,NULL,IFunction,&user);
432: TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
433: TSSetType(user.ts,TSCN);
434: } else {
435: TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
436: TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
437: TSSetType(user.ts,TSRK);
438: }
439: TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
440: TSSetMaxTime(user.ts,user.ftime);
441: TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);
443: if (monitor) {
444: TSMonitorSet(user.ts,Monitor,&user,NULL);
445: }
447: /* Set ODE initial conditions */
448: VecGetArray(user.U,&x_ptr);
449: x_ptr[0] = 2.0;
450: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
451: VecRestoreArray(user.U,&x_ptr);
452: TSSetTimeStep(user.ts,PetscRealConstant(0.001));
454: /* Set runtime options */
455: TSSetFromOptions(user.ts);
457: TSSolve(user.ts,user.U);
458: VecGetArrayRead(user.U,&y_ptr);
459: user.ob[0] = y_ptr[0];
460: user.ob[1] = y_ptr[1];
461: VecRestoreArrayRead(user.U,&y_ptr);
463: /* Save trajectory of solution so that TSAdjointSolve() may be used.
464: Skip checkpointing for the first TSSolve since no adjoint run follows it.
465: */
466: TSSetSaveTrajectory(user.ts);
468: /* Optimization starts */
469: MatCreate(PETSC_COMM_WORLD,&user.H);
470: MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);
471: MatSetUp(user.H); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */
473: /* Set initial solution guess */
474: MatCreateVecs(user.Jacp,&P,NULL);
475: VecGetArray(P,&x_ptr);
476: x_ptr[0] = PetscRealConstant(1.2);
477: VecRestoreArray(P,&x_ptr);
478: TaoSetSolution(tao,P);
480: /* Set routine for function and gradient evaluation */
481: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
482: TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
484: /* Check for any TAO command line options */
485: TaoGetKSP(tao,&ksp);
486: if (ksp) {
487: KSPGetPC(ksp,&pc);
488: PCSetType(pc,PCNONE);
489: }
490: TaoSetFromOptions(tao);
492: TaoSolve(tao);
494: VecView(P,PETSC_VIEWER_STDOUT_WORLD);
495: /* Free TAO data structures */
496: TaoDestroy(&tao);
498: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
499: Free work space. All PETSc objects should be destroyed when they
500: are no longer needed.
501: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
502: MatDestroy(&user.H);
503: MatDestroy(&user.A);
504: VecDestroy(&user.U);
505: MatDestroy(&user.Jacp);
506: VecDestroy(&user.Lambda[0]);
507: VecDestroy(&user.Mup[0]);
508: VecDestroy(&user.Lambda2[0]);
509: VecDestroy(&user.Mup2[0]);
510: VecDestroy(&user.Ihp1[0]);
511: VecDestroy(&user.Ihp2[0]);
512: VecDestroy(&user.Ihp3[0]);
513: VecDestroy(&user.Ihp4[0]);
514: VecDestroy(&user.Dir);
515: TSDestroy(&user.ts);
516: VecDestroy(&P);
517: PetscFinalize();
518: return 0;
519: }
521: /* ------------------------------------------------------------------ */
522: /*
523: FormFunctionGradient - Evaluates the function and corresponding gradient.
525: Input Parameters:
526: tao - the Tao context
527: X - the input vector
528: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
530: Output Parameters:
531: f - the newly evaluated function
532: G - the newly evaluated gradient
533: */
534: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
535: {
536: User user_ptr = (User)ctx;
537: TS ts = user_ptr->ts;
538: PetscScalar *x_ptr,*g;
539: const PetscScalar *y_ptr;
542: VecGetArrayRead(P,&y_ptr);
543: user_ptr->mu = y_ptr[0];
544: VecRestoreArrayRead(P,&y_ptr);
546: TSSetTime(ts,0.0);
547: TSSetStepNumber(ts,0);
548: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
549: TSSetFromOptions(ts);
550: VecGetArray(user_ptr->U,&x_ptr);
551: x_ptr[0] = 2.0;
552: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
553: VecRestoreArray(user_ptr->U,&x_ptr);
555: TSSolve(ts,user_ptr->U);
557: VecGetArrayRead(user_ptr->U,&y_ptr);
558: *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]);
560: /* Reset initial conditions for the adjoint integration */
561: VecGetArray(user_ptr->Lambda[0],&x_ptr);
562: x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
563: x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
564: VecRestoreArrayRead(user_ptr->U,&y_ptr);
565: VecRestoreArray(user_ptr->Lambda[0],&x_ptr);
567: VecGetArray(user_ptr->Mup[0],&x_ptr);
568: x_ptr[0] = 0.0;
569: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
570: TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);
572: TSAdjointSolve(ts);
574: VecGetArray(user_ptr->Mup[0],&x_ptr);
575: VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
576: VecGetArray(G,&g);
577: 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];
578: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
579: VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
580: VecRestoreArray(G,&g);
581: return 0;
582: }
584: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
585: {
586: User user_ptr = (User)ctx;
587: PetscScalar harr[1];
588: const PetscInt rows[1] = {0};
589: PetscInt col = 0;
592: Adjoint2(P,harr,user_ptr);
593: MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);
595: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
596: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
597: if (H != Hpre) {
598: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
599: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
600: }
601: return 0;
602: }
604: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
605: {
606: TS ts = ctx->ts;
607: const PetscScalar *z_ptr;
608: PetscScalar *x_ptr,*y_ptr,dzdp,dzdp2;
609: Mat tlmsen;
612: /* Reset TSAdjoint so that AdjointSetUp will be called again */
613: TSAdjointReset(ts);
615: /* The directional vector should be 1 since it is one-dimensional */
616: VecGetArray(ctx->Dir,&x_ptr);
617: x_ptr[0] = 1.;
618: VecRestoreArray(ctx->Dir,&x_ptr);
620: VecGetArrayRead(P,&z_ptr);
621: ctx->mu = z_ptr[0];
622: VecRestoreArrayRead(P,&z_ptr);
624: dzdp = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
625: 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);
627: TSSetTime(ts,0.0);
628: TSSetStepNumber(ts,0);
629: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
630: TSSetFromOptions(ts);
631: TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);
633: MatZeroEntries(ctx->Jacp);
634: MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
635: MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
636: MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);
638: TSAdjointSetForward(ts,ctx->Jacp);
639: VecGetArray(ctx->U,&y_ptr);
640: y_ptr[0] = 2.0;
641: y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
642: VecRestoreArray(ctx->U,&y_ptr);
643: TSSolve(ts,ctx->U);
645: /* Set terminal conditions for first- and second-order adjonts */
646: VecGetArrayRead(ctx->U,&z_ptr);
647: VecGetArray(ctx->Lambda[0],&y_ptr);
648: y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
649: y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
650: VecRestoreArray(ctx->Lambda[0],&y_ptr);
651: VecRestoreArrayRead(ctx->U,&z_ptr);
652: VecGetArray(ctx->Mup[0],&y_ptr);
653: y_ptr[0] = 0.0;
654: VecRestoreArray(ctx->Mup[0],&y_ptr);
655: TSForwardGetSensitivities(ts,NULL,&tlmsen);
656: MatDenseGetColumn(tlmsen,0,&x_ptr);
657: VecGetArray(ctx->Lambda2[0],&y_ptr);
658: y_ptr[0] = 2.*x_ptr[0];
659: y_ptr[1] = 2.*x_ptr[1];
660: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
661: VecGetArray(ctx->Mup2[0],&y_ptr);
662: y_ptr[0] = 0.0;
663: VecRestoreArray(ctx->Mup2[0],&y_ptr);
664: MatDenseRestoreColumn(tlmsen,&x_ptr);
665: TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
666: if (ctx->implicitform) {
667: TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
668: } else {
669: TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
670: }
671: TSAdjointSolve(ts);
673: VecGetArray(ctx->Lambda[0],&x_ptr);
674: VecGetArray(ctx->Lambda2[0],&y_ptr);
675: VecGetArrayRead(ctx->Mup2[0],&z_ptr);
677: arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];
679: VecRestoreArray(ctx->Lambda2[0],&x_ptr);
680: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
681: VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);
683: /* Disable second-order adjoint mode */
684: TSAdjointReset(ts);
685: TSAdjointResetForward(ts);
686: return 0;
687: }
689: /*TEST
690: build:
691: requires: !complex !single
692: test:
693: args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
694: output_file: output/ex20opt_p_1.out
696: test:
697: suffix: 2
698: 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
699: output_file: output/ex20opt_p_2.out
701: test:
702: suffix: 3
703: args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
704: output_file: output/ex20opt_p_3.out
706: test:
707: suffix: 4
708: 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
709: output_file: output/ex20opt_p_4.out
711: TEST*/