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: PetscErrorCode ierr;
53: User user = (User)ctx;
54: PetscScalar *f;
55: const PetscScalar *u;
58: VecGetArrayRead(U,&u);
59: VecGetArray(F,&f);
60: f[0] = u[1];
61: f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
62: VecRestoreArrayRead(U,&u);
63: VecRestoreArray(F,&f);
64: return(0);
65: }
67: static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx)
68: {
69: PetscErrorCode ierr;
70: User user = (User)ctx;
71: PetscReal mu = user->mu;
72: PetscInt rowcol[] = {0,1};
73: PetscScalar J[2][2];
74: const PetscScalar *u;
77: VecGetArrayRead(U,&u);
79: J[0][0] = 0;
80: J[1][0] = -mu*(2.0*u[1]*u[0]+1.);
81: J[0][1] = 1.0;
82: J[1][1] = mu*(1.0-u[0]*u[0]);
83: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
84: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
85: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
86: if (B && A != B) {
87: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
88: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
89: }
91: VecRestoreArrayRead(U,&u);
92: return(0);
93: }
95: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
96: {
97: const PetscScalar *vl,*vr,*u;
98: PetscScalar *vhv;
99: PetscScalar dJdU[2][2][2]={{{0}}};
100: PetscInt i,j,k;
101: User user = (User)ctx;
102: PetscErrorCode ierr;
105: VecGetArrayRead(U,&u);
106: VecGetArrayRead(Vl[0],&vl);
107: VecGetArrayRead(Vr,&vr);
108: VecGetArray(VHV[0],&vhv);
110: dJdU[1][0][0] = -2.*user->mu*u[1];
111: dJdU[1][1][0] = -2.*user->mu*u[0];
112: dJdU[1][0][1] = -2.*user->mu*u[0];
113: for (j=0; j<2; j++) {
114: vhv[j] = 0;
115: for (k=0; k<2; k++)
116: for (i=0; i<2; i++)
117: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
118: }
120: VecRestoreArrayRead(U,&u);
121: VecRestoreArrayRead(Vl[0],&vl);
122: VecRestoreArrayRead(Vr,&vr);
123: VecRestoreArray(VHV[0],&vhv);
124: return(0);
125: }
127: static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
128: {
129: const PetscScalar *vl,*vr,*u;
130: PetscScalar *vhv;
131: PetscScalar dJdP[2][2][1]={{{0}}};
132: PetscInt i,j,k;
133: PetscErrorCode ierr;
136: VecGetArrayRead(U,&u);
137: VecGetArrayRead(Vl[0],&vl);
138: VecGetArrayRead(Vr,&vr);
139: VecGetArray(VHV[0],&vhv);
141: dJdP[1][0][0] = -(1.+2.*u[0]*u[1]);
142: dJdP[1][1][0] = 1.-u[0]*u[0];
143: for (j=0; j<2; j++) {
144: vhv[j] = 0;
145: for (k=0; k<1; k++)
146: for (i=0; i<2; i++)
147: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
148: }
150: VecRestoreArrayRead(U,&u);
151: VecRestoreArrayRead(Vl[0],&vl);
152: VecRestoreArrayRead(Vr,&vr);
153: VecRestoreArray(VHV[0],&vhv);
154: return(0);
155: }
157: static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
158: {
159: const PetscScalar *vl,*vr,*u;
160: PetscScalar *vhv;
161: PetscScalar dJdU[2][1][2]={{{0}}};
162: PetscInt i,j,k;
163: PetscErrorCode ierr;
166: VecGetArrayRead(U,&u);
167: VecGetArrayRead(Vl[0],&vl);
168: VecGetArrayRead(Vr,&vr);
169: VecGetArray(VHV[0],&vhv);
171: dJdU[1][0][0] = -1.-2.*u[1]*u[0];
172: dJdU[1][0][1] = 1.-u[0]*u[0];
173: for (j=0; j<1; j++) {
174: vhv[j] = 0;
175: for (k=0; k<2; k++)
176: for (i=0; i<2; i++)
177: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
178: }
180: VecRestoreArrayRead(U,&u);
181: VecRestoreArrayRead(Vl[0],&vl);
182: VecRestoreArrayRead(Vr,&vr);
183: VecRestoreArray(VHV[0],&vhv);
184: return(0);
185: }
187: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
188: {
190: return(0);
191: }
193: /* ----------------------- Implicit form of the ODE -------------------- */
195: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx)
196: {
197: PetscErrorCode ierr;
198: User user = (User)ctx;
199: PetscScalar *f;
200: const PetscScalar *u,*udot;
203: VecGetArrayRead(U,&u);
204: VecGetArrayRead(Udot,&udot);
205: VecGetArray(F,&f);
207: f[0] = udot[0] - u[1];
208: f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
210: VecRestoreArrayRead(U,&u);
211: VecRestoreArrayRead(Udot,&udot);
212: VecRestoreArray(F,&f);
213: return(0);
214: }
216: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx)
217: {
218: PetscErrorCode ierr;
219: User user = (User)ctx;
220: PetscInt rowcol[] = {0,1};
221: PetscScalar J[2][2];
222: const PetscScalar *u;
225: VecGetArrayRead(U,&u);
227: J[0][0] = a; J[0][1] = -1.0;
228: 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]);
229: MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
230: VecRestoreArrayRead(U,&u);
231: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
232: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
233: if (A != B) {
234: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
235: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
236: }
238: VecRestoreArrayRead(U,&u);
239: return(0);
240: }
242: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
243: {
244: PetscErrorCode ierr;
245: PetscInt row[] = {0,1},col[]={0};
246: PetscScalar J[2][1];
247: const PetscScalar *u;
250: VecGetArrayRead(U,&u);
252: J[0][0] = 0;
253: J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
254: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
255: VecRestoreArrayRead(U,&u);
256: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
257: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
259: VecRestoreArrayRead(U,&u);
260: return(0);
261: }
263: static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
264: {
265: const PetscScalar *vl,*vr,*u;
266: PetscScalar *vhv;
267: PetscScalar dJdU[2][2][2]={{{0}}};
268: PetscInt i,j,k;
269: User user = (User)ctx;
270: PetscErrorCode ierr;
273: VecGetArrayRead(U,&u);
274: VecGetArrayRead(Vl[0],&vl);
275: VecGetArrayRead(Vr,&vr);
276: VecGetArray(VHV[0],&vhv);
278: dJdU[1][0][0] = 2.*user->mu*u[1];
279: dJdU[1][1][0] = 2.*user->mu*u[0];
280: dJdU[1][0][1] = 2.*user->mu*u[0];
281: for (j=0; j<2; j++) {
282: vhv[j] = 0;
283: for (k=0; k<2; k++)
284: for (i=0; i<2; i++)
285: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
286: }
288: VecRestoreArrayRead(U,&u);
289: VecRestoreArrayRead(Vl[0],&vl);
290: VecRestoreArrayRead(Vr,&vr);
291: VecRestoreArray(VHV[0],&vhv);
292: return(0);
293: }
295: static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
296: {
297: const PetscScalar *vl,*vr,*u;
298: PetscScalar *vhv;
299: PetscScalar dJdP[2][2][1]={{{0}}};
300: PetscInt i,j,k;
301: PetscErrorCode ierr;
304: VecGetArrayRead(U,&u);
305: VecGetArrayRead(Vl[0],&vl);
306: VecGetArrayRead(Vr,&vr);
307: VecGetArray(VHV[0],&vhv);
309: dJdP[1][0][0] = 1.+2.*u[0]*u[1];
310: dJdP[1][1][0] = u[0]*u[0]-1.;
311: for (j=0; j<2; j++) {
312: vhv[j] = 0;
313: for (k=0; k<1; k++)
314: for (i=0; i<2; i++)
315: vhv[j] += vl[i]*dJdP[i][j][k]*vr[k];
316: }
318: VecRestoreArrayRead(U,&u);
319: VecRestoreArrayRead(Vl[0],&vl);
320: VecRestoreArrayRead(Vr,&vr);
321: VecRestoreArray(VHV[0],&vhv);
322: return(0);
323: }
325: static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
326: {
327: const PetscScalar *vl,*vr,*u;
328: PetscScalar *vhv;
329: PetscScalar dJdU[2][1][2]={{{0}}};
330: PetscInt i,j,k;
331: PetscErrorCode ierr;
334: VecGetArrayRead(U,&u);
335: VecGetArrayRead(Vl[0],&vl);
336: VecGetArrayRead(Vr,&vr);
337: VecGetArray(VHV[0],&vhv);
339: dJdU[1][0][0] = 1.+2.*u[1]*u[0];
340: dJdU[1][0][1] = u[0]*u[0]-1.;
341: for (j=0; j<1; j++) {
342: vhv[j] = 0;
343: for (k=0; k<2; k++)
344: for (i=0; i<2; i++)
345: vhv[j] += vl[i]*dJdU[i][j][k]*vr[k];
346: }
348: VecRestoreArrayRead(U,&u);
349: VecRestoreArrayRead(Vl[0],&vl);
350: VecRestoreArrayRead(Vr,&vr);
351: VecRestoreArray(VHV[0],&vhv);
352: return(0);
353: }
355: static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
356: {
358: return(0);
359: }
361: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
362: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
363: {
364: PetscErrorCode ierr;
365: const PetscScalar *x;
366: PetscReal tfinal, dt;
367: User user = (User)ctx;
368: Vec interpolatedX;
371: TSGetTimeStep(ts,&dt);
372: TSGetMaxTime(ts,&tfinal);
374: while (user->next_output <= t && user->next_output <= tfinal) {
375: VecDuplicate(X,&interpolatedX);
376: TSInterpolate(ts,user->next_output,interpolatedX);
377: VecGetArrayRead(interpolatedX,&x);
378: PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n",
379: (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
380: (double)PetscRealPart(x[1]));
381: VecRestoreArrayRead(interpolatedX,&x);
382: VecDestroy(&interpolatedX);
383: user->next_output += PetscRealConstant(0.1);
384: }
385: return(0);
386: }
388: int main(int argc,char **argv)
389: {
390: Vec P;
391: PetscBool monitor = PETSC_FALSE;
392: PetscScalar *x_ptr;
393: const PetscScalar *y_ptr;
394: PetscMPIInt size;
395: struct _n_User user;
396: PetscErrorCode ierr;
397: Tao tao;
398: KSP ksp;
399: PC pc;
401: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402: Initialize program
403: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
405: MPI_Comm_size(PETSC_COMM_WORLD,&size);
406: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
408: /* Create TAO solver and set desired solution method */
409: TaoCreate(PETSC_COMM_WORLD,&tao);
410: TaoSetType(tao,TAOBQNLS);
412: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
413: Set runtime options
414: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
415: user.next_output = 0.0;
416: user.mu = PetscRealConstant(1.0e3);
417: user.ftime = PetscRealConstant(0.5);
418: user.implicitform = PETSC_TRUE;
420: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
421: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
422: PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);
424: /* Create necessary matrix and vectors, solve same ODE on every process */
425: MatCreate(PETSC_COMM_WORLD,&user.A);
426: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
427: MatSetFromOptions(user.A);
428: MatSetUp(user.A);
429: MatCreateVecs(user.A,&user.U,NULL);
430: MatCreateVecs(user.A,&user.Lambda[0],NULL);
431: MatCreateVecs(user.A,&user.Lambda2[0],NULL);
432: MatCreateVecs(user.A,&user.Ihp1[0],NULL);
433: MatCreateVecs(user.A,&user.Ihp2[0],NULL);
435: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
436: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
437: MatSetFromOptions(user.Jacp);
438: MatSetUp(user.Jacp);
439: MatCreateVecs(user.Jacp,&user.Dir,NULL);
440: MatCreateVecs(user.Jacp,&user.Mup[0],NULL);
441: MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);
442: MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);
443: MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);
445: /* Create timestepping solver context */
446: TSCreate(PETSC_COMM_WORLD,&user.ts);
447: TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
448: if (user.implicitform) {
449: TSSetIFunction(user.ts,NULL,IFunction,&user);
450: TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
451: TSSetType(user.ts,TSCN);
452: } else {
453: TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
454: TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
455: TSSetType(user.ts,TSRK);
456: }
457: TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
458: TSSetMaxTime(user.ts,user.ftime);
459: TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);
461: if (monitor) {
462: TSMonitorSet(user.ts,Monitor,&user,NULL);
463: }
465: /* Set ODE initial conditions */
466: VecGetArray(user.U,&x_ptr);
467: x_ptr[0] = 2.0;
468: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
469: VecRestoreArray(user.U,&x_ptr);
470: TSSetTimeStep(user.ts,PetscRealConstant(0.001));
472: /* Set runtime options */
473: TSSetFromOptions(user.ts);
475: TSSolve(user.ts,user.U);
476: VecGetArrayRead(user.U,&y_ptr);
477: user.ob[0] = y_ptr[0];
478: user.ob[1] = y_ptr[1];
479: VecRestoreArrayRead(user.U,&y_ptr);
481: /* Save trajectory of solution so that TSAdjointSolve() may be used.
482: Skip checkpointing for the first TSSolve since no adjoint run follows it.
483: */
484: TSSetSaveTrajectory(user.ts);
486: /* Optimization starts */
487: MatCreate(PETSC_COMM_WORLD,&user.H);
488: MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);
489: MatSetUp(user.H); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */
491: /* Set initial solution guess */
492: MatCreateVecs(user.Jacp,&P,NULL);
493: VecGetArray(P,&x_ptr);
494: x_ptr[0] = PetscRealConstant(1.2);
495: VecRestoreArray(P,&x_ptr);
496: TaoSetInitialVector(tao,P);
498: /* Set routine for function and gradient evaluation */
499: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
500: TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);
502: /* Check for any TAO command line options */
503: TaoGetKSP(tao,&ksp);
504: if (ksp) {
505: KSPGetPC(ksp,&pc);
506: PCSetType(pc,PCNONE);
507: }
508: TaoSetFromOptions(tao);
510: TaoSolve(tao);
512: VecView(P,PETSC_VIEWER_STDOUT_WORLD);
513: /* Free TAO data structures */
514: TaoDestroy(&tao);
516: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
517: Free work space. All PETSc objects should be destroyed when they
518: are no longer needed.
519: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
520: MatDestroy(&user.H);
521: MatDestroy(&user.A);
522: VecDestroy(&user.U);
523: MatDestroy(&user.Jacp);
524: VecDestroy(&user.Lambda[0]);
525: VecDestroy(&user.Mup[0]);
526: VecDestroy(&user.Lambda2[0]);
527: VecDestroy(&user.Mup2[0]);
528: VecDestroy(&user.Ihp1[0]);
529: VecDestroy(&user.Ihp2[0]);
530: VecDestroy(&user.Ihp3[0]);
531: VecDestroy(&user.Ihp4[0]);
532: VecDestroy(&user.Dir);
533: TSDestroy(&user.ts);
534: VecDestroy(&P);
535: PetscFinalize();
536: return ierr;
537: }
539: /* ------------------------------------------------------------------ */
540: /*
541: FormFunctionGradient - Evaluates the function and corresponding gradient.
543: Input Parameters:
544: tao - the Tao context
545: X - the input vector
546: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
548: Output Parameters:
549: f - the newly evaluated function
550: G - the newly evaluated gradient
551: */
552: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
553: {
554: User user_ptr = (User)ctx;
555: TS ts = user_ptr->ts;
556: PetscScalar *x_ptr,*g;
557: const PetscScalar *y_ptr;
558: PetscErrorCode ierr;
561: VecGetArrayRead(P,&y_ptr);
562: user_ptr->mu = y_ptr[0];
563: VecRestoreArrayRead(P,&y_ptr);
565: TSSetTime(ts,0.0);
566: TSSetStepNumber(ts,0);
567: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
568: TSSetFromOptions(ts);
569: VecGetArray(user_ptr->U,&x_ptr);
570: x_ptr[0] = 2.0;
571: x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
572: VecRestoreArray(user_ptr->U,&x_ptr);
574: TSSolve(ts,user_ptr->U);
576: VecGetArrayRead(user_ptr->U,&y_ptr);
577: *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]);
579: /* Reset initial conditions for the adjoint integration */
580: VecGetArray(user_ptr->Lambda[0],&x_ptr);
581: x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
582: x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
583: VecRestoreArrayRead(user_ptr->U,&y_ptr);
584: VecRestoreArray(user_ptr->Lambda[0],&x_ptr);
586: VecGetArray(user_ptr->Mup[0],&x_ptr);
587: x_ptr[0] = 0.0;
588: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
589: TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);
591: TSAdjointSolve(ts);
593: VecGetArray(user_ptr->Mup[0],&x_ptr);
594: VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
595: VecGetArray(G,&g);
596: 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];
597: VecRestoreArray(user_ptr->Mup[0],&x_ptr);
598: VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
599: VecRestoreArray(G,&g);
600: return(0);
601: }
603: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
604: {
605: User user_ptr = (User)ctx;
606: PetscScalar harr[1];
607: const PetscInt rows[1] = {0};
608: PetscInt col = 0;
612: Adjoint2(P,harr,user_ptr);
613: MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);
615: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
616: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
617: if (H != Hpre) {
618: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
619: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
620: }
621: return(0);
622: }
624: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
625: {
626: TS ts = ctx->ts;
627: const PetscScalar *z_ptr;
628: PetscScalar *x_ptr,*y_ptr,dzdp,dzdp2;
629: Mat tlmsen;
630: PetscErrorCode ierr;
633: /* Reset TSAdjoint so that AdjointSetUp will be called again */
634: TSAdjointReset(ts);
636: /* The directional vector should be 1 since it is one-dimensional */
637: VecGetArray(ctx->Dir,&x_ptr);
638: x_ptr[0] = 1.;
639: VecRestoreArray(ctx->Dir,&x_ptr);
641: VecGetArrayRead(P,&z_ptr);
642: ctx->mu = z_ptr[0];
643: VecRestoreArrayRead(P,&z_ptr);
645: dzdp = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
646: 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);
648: TSSetTime(ts,0.0);
649: TSSetStepNumber(ts,0);
650: TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
651: TSSetFromOptions(ts);
652: TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);
654: MatZeroEntries(ctx->Jacp);
655: MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
656: MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
657: MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);
659: TSAdjointSetForward(ts,ctx->Jacp);
660: VecGetArray(ctx->U,&y_ptr);
661: y_ptr[0] = 2.0;
662: y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
663: VecRestoreArray(ctx->U,&y_ptr);
664: TSSolve(ts,ctx->U);
666: /* Set terminal conditions for first- and second-order adjonts */
667: VecGetArrayRead(ctx->U,&z_ptr);
668: VecGetArray(ctx->Lambda[0],&y_ptr);
669: y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
670: y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
671: VecRestoreArray(ctx->Lambda[0],&y_ptr);
672: VecRestoreArrayRead(ctx->U,&z_ptr);
673: VecGetArray(ctx->Mup[0],&y_ptr);
674: y_ptr[0] = 0.0;
675: VecRestoreArray(ctx->Mup[0],&y_ptr);
676: TSForwardGetSensitivities(ts,NULL,&tlmsen);
677: MatDenseGetColumn(tlmsen,0,&x_ptr);
678: VecGetArray(ctx->Lambda2[0],&y_ptr);
679: y_ptr[0] = 2.*x_ptr[0];
680: y_ptr[1] = 2.*x_ptr[1];
681: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
682: VecGetArray(ctx->Mup2[0],&y_ptr);
683: y_ptr[0] = 0.0;
684: VecRestoreArray(ctx->Mup2[0],&y_ptr);
685: MatDenseRestoreColumn(tlmsen,&x_ptr);
686: TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
687: if (ctx->implicitform) {
688: TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
689: } else {
690: TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
691: }
692: TSAdjointSolve(ts);
694: VecGetArray(ctx->Lambda[0],&x_ptr);
695: VecGetArray(ctx->Lambda2[0],&y_ptr);
696: VecGetArrayRead(ctx->Mup2[0],&z_ptr);
698: arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];
700: VecRestoreArray(ctx->Lambda2[0],&x_ptr);
701: VecRestoreArray(ctx->Lambda2[0],&y_ptr);
702: VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);
704: /* Disable second-order adjoint mode */
705: TSAdjointReset(ts);
706: TSAdjointResetForward(ts);
707: return(0);
708: }
710: /*TEST
711: build:
712: requires: !complex !single
713: test:
714: args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
715: output_file: output/ex20opt_p_1.out
717: test:
718: suffix: 2
719: 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
720: output_file: output/ex20opt_p_2.out
722: test:
723: suffix: 3
724: args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
725: output_file: output/ex20opt_p_3.out
727: test:
728: suffix: 4
729: 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
730: output_file: output/ex20opt_p_4.out
732: TEST*/