Actual source code: ex20opt_p.c
2: static char help[] = "Solves the van der Pol equation.\n\
3: Input parameters include:\n";
5: /*
6: Concepts: TS^time-dependent nonlinear problems
7: Concepts: TS^van der Pol equation DAE equivalent
8: Concepts: TS^Optimization using adjoint sensitivity analysis
9: Processors: 1
10: */
11: /* ------------------------------------------------------------------------
13: Notes:
14: This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
15: The nonlinear problem is written in a DAE equivalent form.
16: The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
17: The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
18: ------------------------------------------------------------------------- */
19: #include <petsctao.h>
20: #include <petscts.h>
22: typedef struct _n_User *User;
23: struct _n_User {
24: TS ts;
25: PetscReal mu;
26: PetscReal next_output;
28: /* Sensitivity analysis support */
29: PetscReal ftime;
30: Mat A; /* Jacobian matrix */
31: Mat Jacp; /* JacobianP matrix */
32: Mat H; /* Hessian matrix for optimization */
33: Vec U,Lambda[1],Mup[1]; /* adjoint variables */
34: Vec Lambda2[1],Mup2[1]; /* second-order adjoint variables */
35: Vec Ihp1[1]; /* working space for Hessian evaluations */
36: Vec Ihp2[1]; /* working space for Hessian evaluations */
37: Vec Ihp3[1]; /* working space for Hessian evaluations */
38: Vec Ihp4[1]; /* working space for Hessian evaluations */
39: Vec Dir; /* direction vector */
40: PetscReal ob[2]; /* observation used by the cost function */
41: PetscBool implicitform; /* implicit ODE? */
42: };
44: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
45: PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
46: PetscErrorCode Adjoint2(Vec,PetscScalar[],User);
48: /* ----------------------- Explicit form of the ODE -------------------- */
50: static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx)
51: {
52: User user = (User)ctx;
53: PetscScalar *f;
54: const PetscScalar *u;
57: VecGetArrayRead(U,&u);
58: VecGetArray(F,&f);
59: f[0] = u[1];
60: f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
61: VecRestoreArrayRead(U,&u);
62: VecRestoreArray(F,&f);
63: return 0;
64: }
66: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
67: {
68: User user = (User)ctx;
69: PetscReal mu = user->mu;
70: PetscInt rowcol[] = {0,1};
71: PetscScalar J[2][2];
72: const PetscScalar *u;
75: VecGetArrayRead(U,&u);
77: J[0][0] = 0;
78: J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
79: J[0][1] = 1.0;
80: J[1][1] = mu*(1.0-u[0]*u[0]);
81: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
82: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
83: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
84: if (B && A != B) {
85: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
87: }
89: VecRestoreArrayRead(U,&u);
90: return 0;
91: }
93: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
94: {
95: const PetscScalar *vl,*vr,*u;
96: PetscScalar *vhv;
97: PetscScalar dJdU[2][2][2]={{{0}}};
98: PetscInt i,j,k;
99: User user = (User)ctx;
102: VecGetArrayRead(U,&u);
103: VecGetArrayRead(Vl[0],&vl);
104: VecGetArrayRead(Vr,&vr);
105: VecGetArray(VHV[0],&vhv);
107: dJdU[1][0][0] = -2.*user->mu*u[1];
108: dJdU[1][1][0] = -2.*user->mu*u[0];
109: dJdU[1][0][1] = -2.*user->mu*u[0];
110: for (j=0; j<2; j++) {
111: vhv[j] = 0;
112: for (k=0; k<2; k++)
113: for (i=0; i<2; i++)
114: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
115: }
117: VecRestoreArrayRead(U,&u);
118: VecRestoreArrayRead(Vl[0],&vl);
119: VecRestoreArrayRead(Vr,&vr);
120: VecRestoreArray(VHV[0],&vhv);
121: return 0;
122: }
124: static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
125: {
126: const PetscScalar *vl,*vr,*u;
127: PetscScalar *vhv;
128: PetscScalar dJdP[2][2][1]={{{0}}};
129: PetscInt i,j,k;
132: VecGetArrayRead(U,&u);
133: VecGetArrayRead(Vl[0],&vl);
134: VecGetArrayRead(Vr,&vr);
135: VecGetArray(VHV[0],&vhv);
137: dJdP[1][0][0] = -(1.+2.*u[0]*u[1]);
138: dJdP[1][1][0] = 1.-u[0]*u[0];
139: for (j=0; j<2; j++) {
140: vhv[j] = 0;
141: for (k=0; k<1; k++)
142: for (i=0; i<2; i++)
143: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
144: }
146: VecRestoreArrayRead(U,&u);
147: VecRestoreArrayRead(Vl[0],&vl);
148: VecRestoreArrayRead(Vr,&vr);
149: VecRestoreArray(VHV[0],&vhv);
150: return 0;
151: }
153: static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
154: {
155: const PetscScalar *vl,*vr,*u;
156: PetscScalar *vhv;
157: PetscScalar dJdU[2][1][2]={{{0}}};
158: PetscInt i,j,k;
161: VecGetArrayRead(U,&u);
162: VecGetArrayRead(Vl[0],&vl);
163: VecGetArrayRead(Vr,&vr);
164: VecGetArray(VHV[0],&vhv);
166: dJdU[1][0][0] = -1.-2.*u[1]*u[0];
167: dJdU[1][0][1] = 1.-u[0]*u[0];
168: for (j=0; j<1; j++) {
169: vhv[j] = 0;
170: for (k=0; k<2; k++)
171: for (i=0; i<2; i++)
172: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
173: }
175: VecRestoreArrayRead(U,&u);
176: VecRestoreArrayRead(Vl[0],&vl);
177: VecRestoreArrayRead(Vr,&vr);
178: VecRestoreArray(VHV[0],&vhv);
179: return 0;
180: }
182: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
183: {
185: return 0;
186: }
188: /* ----------------------- Implicit form of the ODE -------------------- */
190: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
191: {
192: User user = (User)ctx;
193: PetscScalar *f;
194: const PetscScalar *u,*udot;
197: VecGetArrayRead(U,&u);
198: VecGetArrayRead(Udot,&udot);
199: VecGetArray(F,&f);
201: f[0] = udot[0] - u[1];
202: f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
204: VecRestoreArrayRead(U,&u);
205: VecRestoreArrayRead(Udot,&udot);
206: VecRestoreArray(F,&f);
207: return 0;
208: }
210: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
211: {
212: User user = (User)ctx;
213: PetscInt rowcol[] = {0,1};
214: PetscScalar J[2][2];
215: const PetscScalar *u;
218: VecGetArrayRead(U,&u);
220: J[0][0] = a; J[0][1] = -1.0;
221: 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]);
222: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
223: VecRestoreArrayRead(U,&u);
224: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
225: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
226: if (A != B) {
227: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
228: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
229: }
231: VecRestoreArrayRead(U,&u);
232: return 0;
233: }
235: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
236: {
237: PetscInt row[] = {0,1},col[]={0};
238: PetscScalar J[2][1];
239: const PetscScalar *u;
242: VecGetArrayRead(U,&u);
244: J[0][0] = 0;
245: J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
246: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
247: VecRestoreArrayRead(U,&u);
248: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
249: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
251: VecRestoreArrayRead(U,&u);
252: return 0;
253: }
255: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
256: {
257: const PetscScalar *vl,*vr,*u;
258: PetscScalar *vhv;
259: PetscScalar dJdU[2][2][2]={{{0}}};
260: PetscInt i,j,k;
261: User user = (User)ctx;
264: VecGetArrayRead(U,&u);
265: VecGetArrayRead(Vl[0],&vl);
266: VecGetArrayRead(Vr,&vr);
267: VecGetArray(VHV[0],&vhv);
269: dJdU[1][0][0] = 2.*user->mu*u[1];
270: dJdU[1][1][0] = 2.*user->mu*u[0];
271: dJdU[1][0][1] = 2.*user->mu*u[0];
272: for (j=0; j<2; j++) {
273: vhv[j] = 0;
274: for (k=0; k<2; k++)
275: for (i=0; i<2; i++)
276: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
277: }
279: VecRestoreArrayRead(U,&u);
280: VecRestoreArrayRead(Vl[0],&vl);
281: VecRestoreArrayRead(Vr,&vr);
282: VecRestoreArray(VHV[0],&vhv);
283: return 0;
284: }
286: static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
287: {
288: const PetscScalar *vl,*vr,*u;
289: PetscScalar *vhv;
290: PetscScalar dJdP[2][2][1]={{{0}}};
291: PetscInt i,j,k;
294: VecGetArrayRead(U,&u);
295: VecGetArrayRead(Vl[0],&vl);
296: VecGetArrayRead(Vr,&vr);
297: VecGetArray(VHV[0],&vhv);
299: dJdP[1][0][0] = 1.+2.*u[0]*u[1];
300: dJdP[1][1][0] = u[0]*u[0]-1.;
301: for (j=0; j<2; j++) {
302: vhv[j] = 0;
303: for (k=0; k<1; k++)
304: for (i=0; i<2; i++)
305: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
306: }
308: VecRestoreArrayRead(U,&u);
309: VecRestoreArrayRead(Vl[0],&vl);
310: VecRestoreArrayRead(Vr,&vr);
311: VecRestoreArray(VHV[0],&vhv);
312: return 0;
313: }
315: static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
316: {
317: const PetscScalar *vl,*vr,*u;
318: PetscScalar *vhv;
319: PetscScalar dJdU[2][1][2]={{{0}}};
320: PetscInt i,j,k;
323: VecGetArrayRead(U,&u);
324: VecGetArrayRead(Vl[0],&vl);
325: VecGetArrayRead(Vr,&vr);
326: VecGetArray(VHV[0],&vhv);
328: dJdU[1][0][0] = 1.+2.*u[1]*u[0];
329: dJdU[1][0][1] = u[0]*u[0]-1.;
330: for (j=0; j<1; j++) {
331: vhv[j] = 0;
332: for (k=0; k<2; k++)
333: for (i=0; i<2; i++)
334: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
335: }
337: VecRestoreArrayRead(U,&u);
338: VecRestoreArrayRead(Vl[0],&vl);
339: VecRestoreArrayRead(Vr,&vr);
340: VecRestoreArray(VHV[0],&vhv);
341: return 0;
342: }
344: static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
345: {
347: return 0;
348: }
350: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
351: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
352: {
353: PetscErrorCode ierr;
354: const PetscScalar *x;
355: PetscReal tfinal, dt;
356: User user = (User)ctx;
357: Vec interpolatedX;
360: TSGetTimeStep(ts,&dt);
361: TSGetMaxTime(ts,&tfinal);
363: while (user->next_output <= t && user->next_output <= tfinal) {
364: VecDuplicate(X,&interpolatedX);
365: TSInterpolate(ts,user->next_output,interpolatedX);
366: VecGetArrayRead(interpolatedX,&x);
367: PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
368: (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
369: (double)PetscRealPart(x[1]));
370: VecRestoreArrayRead(interpolatedX,&x);
371: VecDestroy(&interpolatedX);
372: user->next_output += PetscRealConstant(0.1);
373: }
374: return 0;
375: }
377: int main(int argc,char **argv)
378: {
379: Vec P;
380: PetscBool monitor = PETSC_FALSE;
381: PetscScalar *x_ptr;
382: const PetscScalar *y_ptr;
383: PetscMPIInt size;
384: struct _n_User user;
385: Tao tao;
386: KSP ksp;
387: PC pc;
389: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
390: Initialize program
391: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
392: PetscInitialize(&argc,&argv,NULL,help);
393: MPI_Comm_size(PETSC_COMM_WORLD,&size);
396: /* Create TAO solver and set desired solution method */
397: TaoCreate(PETSC_COMM_WORLD,&tao);
398: TaoSetType(tao,TAOBQNLS);
400: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
401: Set runtime options
402: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
403: user.next_output = 0.0;
404: user.mu = PetscRealConstant(1.0e3);
405: user.ftime = PetscRealConstant(0.5);
406: user.implicitform = PETSC_TRUE;
408: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
409: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
410: PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);
412: /* Create necessary matrix and vectors, solve same ODE on every process */
413: MatCreate(PETSC_COMM_WORLD,&user.A);
414: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
415: MatSetFromOptions(user.A);
416: MatSetUp(user.A);
417: MatCreateVecs(user.A,&user.U,NULL);
418: MatCreateVecs(user.A,&user.Lambda[0],NULL);
419: MatCreateVecs(user.A,&user.Lambda2[0],NULL);
420: MatCreateVecs(user.A,&user.Ihp1[0],NULL);
421: MatCreateVecs(user.A,&user.Ihp2[0],NULL);
423: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
424: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
425: MatSetFromOptions(user.Jacp);
426: MatSetUp(user.Jacp);
427: MatCreateVecs(user.Jacp,&user.Dir,NULL);
428: MatCreateVecs(user.Jacp,&user.Mup[0],NULL);
429: MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);
430: MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);
431: MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);
433: /* Create timestepping solver context */
434: TSCreate(PETSC_COMM_WORLD,&user.ts);
435: TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
436: if (user.implicitform) {
437: TSSetIFunction(user.ts,NULL,IFunction,&user);
438: TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
439: TSSetType(user.ts,TSCN);
440: } else {
441: TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
442: TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
443: TSSetType(user.ts,TSRK);
444: }
445: TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
446: TSSetMaxTime(user.ts,user.ftime);
447: TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);
449: if (monitor) {
450: TSMonitorSet(user.ts,Monitor,&user,NULL);
451: }
453: /* Set ODE initial conditions */
454: VecGetArray(user.U,&x_ptr);
455: x_ptr[0] = 2.0;
456: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
457: VecRestoreArray(user.U,&x_ptr);
458: TSSetTimeStep(user.ts,PetscRealConstant(0.001));
460: /* Set runtime options */
461: TSSetFromOptions(user.ts);
463: TSSolve(user.ts,user.U);
464: VecGetArrayRead(user.U,&y_ptr);
465: user.ob[0] = y_ptr[0];
466: user.ob[1] = y_ptr[1];
467: VecRestoreArrayRead(user.U,&y_ptr);
469: /* Save trajectory of solution so that TSAdjointSolve() may be used.
470: Skip checkpointing for the first TSSolve since no adjoint run follows it.
471: */
472: TSSetSaveTrajectory(user.ts);
474: /* Optimization starts */
475: MatCreate(PETSC_COMM_WORLD,&user.H);
476: MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);
477: MatSetUp(user.H); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */
479: /* Set initial solution guess */
480: MatCreateVecs(user.Jacp,&P,NULL);
481: VecGetArray(P,&x_ptr);
482: x_ptr[0] = PetscRealConstant(1.2);
483: VecRestoreArray(P,&x_ptr);
484: TaoSetSolution(tao,P);
486: /* Set routine for function and gradient evaluation */
487: TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user);
488: TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user);
490: /* Check for any TAO command line options */
491: TaoGetKSP(tao,&ksp);
492: if (ksp) {
493: KSPGetPC(ksp,&pc);
494: PCSetType(pc,PCNONE);
495: }
496: TaoSetFromOptions(tao);
498: TaoSolve(tao);
500: VecView(P,PETSC_VIEWER_STDOUT_WORLD);
501: /* Free TAO data structures */
502: TaoDestroy(&tao);
504: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
505: Free work space. All PETSc objects should be destroyed when they
506: are no longer needed.
507: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
508: MatDestroy(&user.H);
509: MatDestroy(&user.A);
510: VecDestroy(&user.U);
511: MatDestroy(&user.Jacp);
512: VecDestroy(&user.Lambda[0]);
513: VecDestroy(&user.Mup[0]);
514: VecDestroy(&user.Lambda2[0]);
515: VecDestroy(&user.Mup2[0]);
516: VecDestroy(&user.Ihp1[0]);
517: VecDestroy(&user.Ihp2[0]);
518: VecDestroy(&user.Ihp3[0]);
519: VecDestroy(&user.Ihp4[0]);
520: VecDestroy(&user.Dir);
521: TSDestroy(&user.ts);
522: VecDestroy(&P);
523: PetscFinalize();
524: return 0;
525: }
527: /* ------------------------------------------------------------------ */
528: /*
529: FormFunctionGradient - Evaluates the function and corresponding gradient.
531: Input Parameters:
532: tao - the Tao context
533: X - the input vector
534: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
536: Output Parameters:
537: f - the newly evaluated function
538: G - the newly evaluated gradient
539: */
540: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
541: {
542: User user_ptr = (User)ctx;
543: TS ts = user_ptr->ts;
544: PetscScalar *x_ptr,*g;
545: const PetscScalar *y_ptr;
548: VecGetArrayRead(P,&y_ptr);
549: user_ptr->mu = y_ptr[0];
550: VecRestoreArrayRead(P,&y_ptr);
552: TSSetTime(ts,0.0);
553: TSSetStepNumber(ts,0);
554: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
555: TSSetFromOptions(ts);
556: VecGetArray(user_ptr->U,&x_ptr);
557: x_ptr[0] = 2.0;
558: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
559: VecRestoreArray(user_ptr->U,&x_ptr);
561: TSSolve(ts,user_ptr->U);
563: VecGetArrayRead(user_ptr->U,&y_ptr);
564: *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]);
566: /* Reset initial conditions for the adjoint integration */
567: VecGetArray(user_ptr->Lambda[0],&x_ptr);
568: x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
569: x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
570: VecRestoreArrayRead(user_ptr->U,&y_ptr);
571: VecRestoreArray(user_ptr->Lambda[0],&x_ptr);
573: VecGetArray(user_ptr->Mup[0],&x_ptr);
574: x_ptr[0] = 0.0;
575: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
576: TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);
578: TSAdjointSolve(ts);
580: VecGetArray(user_ptr->Mup[0],&x_ptr);
581: VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
582: VecGetArray(G,&g);
583: 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];
584: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
585: VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
586: VecRestoreArray(G,&g);
587: return 0;
588: }
590: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
591: {
592: User user_ptr = (User)ctx;
593: PetscScalar harr[1];
594: const PetscInt rows[1] = {0};
595: PetscInt col = 0;
598: Adjoint2(P,harr,user_ptr);
599: MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);
601: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
602: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
603: if (H != Hpre) {
604: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
605: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
606: }
607: return 0;
608: }
610: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
611: {
612: TS ts = ctx->ts;
613: const PetscScalar *z_ptr;
614: PetscScalar *x_ptr,*y_ptr,dzdp,dzdp2;
615: Mat tlmsen;
618: /* Reset TSAdjoint so that AdjointSetUp will be called again */
619: TSAdjointReset(ts);
621: /* The directional vector should be 1 since it is one-dimensional */
622: VecGetArray(ctx->Dir,&x_ptr);
623: x_ptr[0] = 1.;
624: VecRestoreArray(ctx->Dir,&x_ptr);
626: VecGetArrayRead(P,&z_ptr);
627: ctx->mu = z_ptr[0];
628: VecRestoreArrayRead(P,&z_ptr);
630: dzdp = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
631: 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);
633: TSSetTime(ts,0.0);
634: TSSetStepNumber(ts,0);
635: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
636: TSSetFromOptions(ts);
637: TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);
639: MatZeroEntries(ctx->Jacp);
640: MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
641: MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
642: MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);
644: TSAdjointSetForward(ts,ctx->Jacp);
645: VecGetArray(ctx->U,&y_ptr);
646: y_ptr[0] = 2.0;
647: y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
648: VecRestoreArray(ctx->U,&y_ptr);
649: TSSolve(ts,ctx->U);
651: /* Set terminal conditions for first- and second-order adjonts */
652: VecGetArrayRead(ctx->U,&z_ptr);
653: VecGetArray(ctx->Lambda[0],&y_ptr);
654: y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
655: y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
656: VecRestoreArray(ctx->Lambda[0],&y_ptr);
657: VecRestoreArrayRead(ctx->U,&z_ptr);
658: VecGetArray(ctx->Mup[0],&y_ptr);
659: y_ptr[0] = 0.0;
660: VecRestoreArray(ctx->Mup[0],&y_ptr);
661: TSForwardGetSensitivities(ts,NULL,&tlmsen);
662: MatDenseGetColumn(tlmsen,0,&x_ptr);
663: VecGetArray(ctx->Lambda2[0],&y_ptr);
664: y_ptr[0] = 2.*x_ptr[0];
665: y_ptr[1] = 2.*x_ptr[1];
666: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
667: VecGetArray(ctx->Mup2[0],&y_ptr);
668: y_ptr[0] = 0.0;
669: VecRestoreArray(ctx->Mup2[0],&y_ptr);
670: MatDenseRestoreColumn(tlmsen,&x_ptr);
671: TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
672: if (ctx->implicitform) {
673: TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
674: } else {
675: TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
676: }
677: TSAdjointSolve(ts);
679: VecGetArray(ctx->Lambda[0],&x_ptr);
680: VecGetArray(ctx->Lambda2[0],&y_ptr);
681: VecGetArrayRead(ctx->Mup2[0],&z_ptr);
683: arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];
685: VecRestoreArray(ctx->Lambda2[0],&x_ptr);
686: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
687: VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);
689: /* Disable second-order adjoint mode */
690: TSAdjointReset(ts);
691: TSAdjointResetForward(ts);
692: return 0;
693: }
695: /*TEST
696: build:
697: requires: !complex !single
698: test:
699: args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
700: output_file: output/ex20opt_p_1.out
702: test:
703: suffix: 2
704: 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
705: output_file: output/ex20opt_p_2.out
707: test:
708: suffix: 3
709: args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
710: output_file: output/ex20opt_p_3.out
712: test:
713: suffix: 4
714: 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
715: output_file: output/ex20opt_p_4.out
717: TEST*/