xref: /petsc/src/ts/tutorials/ex20opt_p.c (revision 2da392cc7c10228af19ad9843ce5155178acb644) 
1 
2 static char help[] = "Solves the van der Pol equation.\n\
3 Input parameters include:\n";
4 
5 /*
6    Concepts: TS^time-dependent nonlinear problems
7    Concepts: TS^van der Pol equation DAE equivalent
8    Concepts: Optimization using adjoint sensitivity analysis
9    Processors: 1
10 */
11 /* ------------------------------------------------------------------------
12 
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>
21 
22 typedef struct _n_User *User;
23 struct _n_User {
24   TS        ts;
25   PetscReal mu;
26   PetscReal next_output;
27 
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 };
43 
44 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
45 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*);
46 PetscErrorCode Adjoint2(Vec,PetscScalar[],User);
47 
48 /* ----------------------- Explicit form of the ODE  -------------------- */
49 
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;
56 
57   PetscFunctionBeginUser;
58   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
59   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
60   f[0] = u[1];
61   f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]);
62   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
63   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
64   PetscFunctionReturn(0);
65 }
66 
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;
75 
76   PetscFunctionBeginUser;
77   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
78 
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   ierr    = MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
84   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
85   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
86   if (B && A != B) {
87     ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
88     ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
89   }
90 
91   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
92   PetscFunctionReturn(0);
93 }
94 
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;
103 
104   PetscFunctionBeginUser;
105   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
106   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
107   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
108   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
109 
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   }
119 
120   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
121   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
122   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
123   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
124   PetscFunctionReturn(0);
125 }
126 
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;
134 
135   PetscFunctionBeginUser;
136   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
137   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
138   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
139   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
140 
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   }
149 
150   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
151   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
152   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
153   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
154   PetscFunctionReturn(0);
155 }
156 
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;
164 
165   PetscFunctionBeginUser;
166   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
167   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
168   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
169   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
170 
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   }
179 
180   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
181   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
182   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
183   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
184   PetscFunctionReturn(0);
185 }
186 
187 static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
188 {
189   PetscFunctionBeginUser;
190   PetscFunctionReturn(0);
191 }
192 
193 /* ----------------------- Implicit form of the ODE  -------------------- */
194 
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;
201 
202   PetscFunctionBeginUser;
203   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
204   ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr);
205   ierr = VecGetArray(F,&f);CHKERRQ(ierr);
206 
207   f[0] = udot[0] - u[1];
208   f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ;
209 
210   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
211   ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr);
212   ierr = VecRestoreArray(F,&f);CHKERRQ(ierr);
213   PetscFunctionReturn(0);
214 }
215 
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;
223 
224   PetscFunctionBeginUser;
225   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
226 
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   ierr    = MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
230   ierr    = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
231   ierr    = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
232   ierr    = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
233   if (A != B) {
234     ierr  = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
235     ierr  = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
236   }
237 
238   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
239   PetscFunctionReturn(0);
240 }
241 
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;
248 
249   PetscFunctionBeginUser;
250   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
251 
252   J[0][0] = 0;
253   J[1][0] = (1.-u[0]*u[0])*u[1]-u[0];
254   ierr    = MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);CHKERRQ(ierr);
255   ierr    = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
256   ierr    = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
257   ierr    = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
258 
259   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
260   PetscFunctionReturn(0);
261 }
262 
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;
271 
272   PetscFunctionBeginUser;
273   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
274   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
275   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
276   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
277 
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   }
287 
288   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
289   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
290   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
291   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
292   PetscFunctionReturn(0);
293 }
294 
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;
302 
303   PetscFunctionBeginUser;
304   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
305   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
306   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
307   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
308 
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   }
317 
318   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
319   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
320   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
321   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
322   PetscFunctionReturn(0);
323 }
324 
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;
332 
333   PetscFunctionBeginUser;
334   ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr);
335   ierr = VecGetArrayRead(Vl[0],&vl);CHKERRQ(ierr);
336   ierr = VecGetArrayRead(Vr,&vr);CHKERRQ(ierr);
337   ierr = VecGetArray(VHV[0],&vhv);CHKERRQ(ierr);
338 
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   }
347 
348   ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr);
349   ierr = VecRestoreArrayRead(Vl[0],&vl);CHKERRQ(ierr);
350   ierr = VecRestoreArrayRead(Vr,&vr);CHKERRQ(ierr);
351   ierr = VecRestoreArray(VHV[0],&vhv);CHKERRQ(ierr);
352   PetscFunctionReturn(0);
353 }
354 
355 static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
356 {
357   PetscFunctionBeginUser;
358   PetscFunctionReturn(0);
359 }
360 
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;
369 
370   PetscFunctionBeginUser;
371   ierr = TSGetTimeStep(ts,&dt);CHKERRQ(ierr);
372   ierr = TSGetMaxTime(ts,&tfinal);CHKERRQ(ierr);
373 
374   while (user->next_output <= t && user->next_output <= tfinal) {
375     ierr = VecDuplicate(X,&interpolatedX);CHKERRQ(ierr);
376     ierr = TSInterpolate(ts,user->next_output,interpolatedX);CHKERRQ(ierr);
377     ierr = VecGetArrayRead(interpolatedX,&x);CHKERRQ(ierr);
378     ierr = 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]));CHKERRQ(ierr);
381     ierr = VecRestoreArrayRead(interpolatedX,&x);CHKERRQ(ierr);
382     ierr = VecDestroy(&interpolatedX);CHKERRQ(ierr);
383     user->next_output += PetscRealConstant(0.1);
384   }
385   PetscFunctionReturn(0);
386 }
387 
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;
400 
401   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
402      Initialize program
403      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
404   ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
405   ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);
406   if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
407 
408   /* Create TAO solver and set desired solution method */
409   ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
410   ierr = TaoSetType(tao,TAOBQNLS);CHKERRQ(ierr);
411 
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;
419 
420   ierr = PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);CHKERRQ(ierr);
421   ierr = PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);CHKERRQ(ierr);
422   ierr = PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL);CHKERRQ(ierr);
423 
424   /* Create necessary matrix and vectors, solve same ODE on every process */
425   ierr = MatCreate(PETSC_COMM_WORLD,&user.A);CHKERRQ(ierr);
426   ierr = MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
427   ierr = MatSetFromOptions(user.A);CHKERRQ(ierr);
428   ierr = MatSetUp(user.A);CHKERRQ(ierr);
429   ierr = MatCreateVecs(user.A,&user.U,NULL);CHKERRQ(ierr);
430   ierr = MatCreateVecs(user.A,&user.Lambda[0],NULL);CHKERRQ(ierr);
431   ierr = MatCreateVecs(user.A,&user.Lambda2[0],NULL);CHKERRQ(ierr);
432   ierr = MatCreateVecs(user.A,&user.Ihp1[0],NULL);CHKERRQ(ierr);
433   ierr = MatCreateVecs(user.A,&user.Ihp2[0],NULL);CHKERRQ(ierr);
434 
435   ierr = MatCreate(PETSC_COMM_WORLD,&user.Jacp);CHKERRQ(ierr);
436   ierr = MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);CHKERRQ(ierr);
437   ierr = MatSetFromOptions(user.Jacp);CHKERRQ(ierr);
438   ierr = MatSetUp(user.Jacp);CHKERRQ(ierr);
439   ierr = MatCreateVecs(user.Jacp,&user.Dir,NULL);CHKERRQ(ierr);
440   ierr = MatCreateVecs(user.Jacp,&user.Mup[0],NULL);CHKERRQ(ierr);
441   ierr = MatCreateVecs(user.Jacp,&user.Mup2[0],NULL);CHKERRQ(ierr);
442   ierr = MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL);CHKERRQ(ierr);
443   ierr = MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL);CHKERRQ(ierr);
444 
445   /* Create timestepping solver context */
446   ierr = TSCreate(PETSC_COMM_WORLD,&user.ts);CHKERRQ(ierr);
447   ierr = TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
448   if (user.implicitform) {
449     ierr = TSSetIFunction(user.ts,NULL,IFunction,&user);CHKERRQ(ierr);
450     ierr = TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);CHKERRQ(ierr);
451     ierr = TSSetType(user.ts,TSCN);CHKERRQ(ierr);
452   } else {
453     ierr = TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);CHKERRQ(ierr);
454     ierr = TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);CHKERRQ(ierr);
455     ierr = TSSetType(user.ts,TSRK);CHKERRQ(ierr);
456   }
457   ierr = TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);CHKERRQ(ierr);
458   ierr = TSSetMaxTime(user.ts,user.ftime);CHKERRQ(ierr);
459   ierr = TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
460 
461   if (monitor) {
462     ierr = TSMonitorSet(user.ts,Monitor,&user,NULL);CHKERRQ(ierr);
463   }
464 
465   /* Set ODE initial conditions */
466   ierr = VecGetArray(user.U,&x_ptr);CHKERRQ(ierr);
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   ierr = VecRestoreArray(user.U,&x_ptr);CHKERRQ(ierr);
470   ierr = TSSetTimeStep(user.ts,PetscRealConstant(0.001));CHKERRQ(ierr);
471 
472   /* Set runtime options */
473   ierr = TSSetFromOptions(user.ts);CHKERRQ(ierr);
474 
475   ierr       = TSSolve(user.ts,user.U);CHKERRQ(ierr);
476   ierr       = VecGetArrayRead(user.U,&y_ptr);CHKERRQ(ierr);
477   user.ob[0] = y_ptr[0];
478   user.ob[1] = y_ptr[1];
479   ierr       = VecRestoreArrayRead(user.U,&y_ptr);CHKERRQ(ierr);
480 
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   ierr = TSSetSaveTrajectory(user.ts);CHKERRQ(ierr);
485 
486   /* Optimization starts */
487   ierr = MatCreate(PETSC_COMM_WORLD,&user.H);CHKERRQ(ierr);
488   ierr = MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);CHKERRQ(ierr);
489   ierr = MatSetUp(user.H);CHKERRQ(ierr); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */
490 
491   /* Set initial solution guess */
492   ierr = MatCreateVecs(user.Jacp,&P,NULL);CHKERRQ(ierr);
493   ierr = VecGetArray(P,&x_ptr);CHKERRQ(ierr);
494   x_ptr[0] = PetscRealConstant(1.2);
495   ierr = VecRestoreArray(P,&x_ptr);CHKERRQ(ierr);
496   ierr = TaoSetInitialVector(tao,P);CHKERRQ(ierr);
497 
498   /* Set routine for function and gradient evaluation */
499   ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);CHKERRQ(ierr);
500   ierr = TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);CHKERRQ(ierr);
501 
502   /* Check for any TAO command line options */
503   ierr = TaoGetKSP(tao,&ksp);CHKERRQ(ierr);
504   if (ksp) {
505     ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
506     ierr = PCSetType(pc,PCNONE);CHKERRQ(ierr);
507   }
508   ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
509 
510   ierr = TaoSolve(tao);CHKERRQ(ierr);
511 
512   ierr = VecView(P,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
513   /* Free TAO data structures */
514   ierr = TaoDestroy(&tao);CHKERRQ(ierr);
515 
516   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
517      Free work space.  All PETSc objects should be destroyed when they
518      are no longer needed.
519    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
520   ierr = MatDestroy(&user.H);CHKERRQ(ierr);
521   ierr = MatDestroy(&user.A);CHKERRQ(ierr);
522   ierr = VecDestroy(&user.U);CHKERRQ(ierr);
523   ierr = MatDestroy(&user.Jacp);CHKERRQ(ierr);
524   ierr = VecDestroy(&user.Lambda[0]);CHKERRQ(ierr);
525   ierr = VecDestroy(&user.Mup[0]);CHKERRQ(ierr);
526   ierr = VecDestroy(&user.Lambda2[0]);CHKERRQ(ierr);
527   ierr = VecDestroy(&user.Mup2[0]);CHKERRQ(ierr);
528   ierr = VecDestroy(&user.Ihp1[0]);CHKERRQ(ierr);
529   ierr = VecDestroy(&user.Ihp2[0]);CHKERRQ(ierr);
530   ierr = VecDestroy(&user.Ihp3[0]);CHKERRQ(ierr);
531   ierr = VecDestroy(&user.Ihp4[0]);CHKERRQ(ierr);
532   ierr = VecDestroy(&user.Dir);CHKERRQ(ierr);
533   ierr = TSDestroy(&user.ts);CHKERRQ(ierr);
534   ierr = VecDestroy(&P);CHKERRQ(ierr);
535   ierr = PetscFinalize();
536   return ierr;
537 }
538 
539 /* ------------------------------------------------------------------ */
540 /*
541    FormFunctionGradient - Evaluates the function and corresponding gradient.
542 
543    Input Parameters:
544    tao - the Tao context
545    X   - the input vector
546    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
547 
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;
559 
560   PetscFunctionBeginUser;
561   ierr = VecGetArrayRead(P,&y_ptr);CHKERRQ(ierr);
562   user_ptr->mu = y_ptr[0];
563   ierr = VecRestoreArrayRead(P,&y_ptr);CHKERRQ(ierr);
564 
565   ierr = TSSetTime(ts,0.0);CHKERRQ(ierr);
566   ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr);
567   ierr = TSSetTimeStep(ts,PetscRealConstant(0.001));CHKERRQ(ierr); /* can be overwritten by command line options */
568   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
569   ierr = VecGetArray(user_ptr->U,&x_ptr);CHKERRQ(ierr);
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   ierr = VecRestoreArray(user_ptr->U,&x_ptr);CHKERRQ(ierr);
573 
574   ierr = TSSolve(ts,user_ptr->U);CHKERRQ(ierr);
575 
576   ierr = VecGetArrayRead(user_ptr->U,&y_ptr);CHKERRQ(ierr);
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]);
578 
579   /*   Reset initial conditions for the adjoint integration */
580   ierr = VecGetArray(user_ptr->Lambda[0],&x_ptr);CHKERRQ(ierr);
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   ierr = VecRestoreArrayRead(user_ptr->U,&y_ptr);CHKERRQ(ierr);
584   ierr = VecRestoreArray(user_ptr->Lambda[0],&x_ptr);CHKERRQ(ierr);
585 
586   ierr = VecGetArray(user_ptr->Mup[0],&x_ptr);CHKERRQ(ierr);
587   x_ptr[0] = 0.0;
588   ierr = VecRestoreArray(user_ptr->Mup[0],&x_ptr);CHKERRQ(ierr);
589   ierr = TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);CHKERRQ(ierr);
590 
591   ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
592 
593   ierr = VecGetArray(user_ptr->Mup[0],&x_ptr);CHKERRQ(ierr);
594   ierr = VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);CHKERRQ(ierr);
595   ierr = VecGetArray(G,&g);CHKERRQ(ierr);
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   ierr = VecRestoreArray(user_ptr->Mup[0],&x_ptr);CHKERRQ(ierr);
598   ierr = VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);CHKERRQ(ierr);
599   ierr = VecRestoreArray(G,&g);CHKERRQ(ierr);
600   PetscFunctionReturn(0);
601 }
602 
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;
609   PetscErrorCode ierr;
610 
611   PetscFunctionBeginUser;
612   ierr = Adjoint2(P,harr,user_ptr);CHKERRQ(ierr);
613   ierr = MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);CHKERRQ(ierr);
614 
615   ierr     = MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
616   ierr     = MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
617   if (H != Hpre) {
618     ierr   = MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
619     ierr   = MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
620   }
621   PetscFunctionReturn(0);
622 }
623 
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;
631 
632   PetscFunctionBeginUser;
633   /* Reset TSAdjoint so that AdjointSetUp will be called again */
634   ierr = TSAdjointReset(ts);CHKERRQ(ierr);
635 
636   /* The directional vector should be 1 since it is one-dimensional */
637   ierr     = VecGetArray(ctx->Dir,&x_ptr);CHKERRQ(ierr);
638   x_ptr[0] = 1.;
639   ierr     = VecRestoreArray(ctx->Dir,&x_ptr);CHKERRQ(ierr);
640 
641   ierr = VecGetArrayRead(P,&z_ptr);CHKERRQ(ierr);
642   ctx->mu = z_ptr[0];
643   ierr = VecRestoreArrayRead(P,&z_ptr);CHKERRQ(ierr);
644 
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);
647 
648   ierr = TSSetTime(ts,0.0);CHKERRQ(ierr);
649   ierr = TSSetStepNumber(ts,0);CHKERRQ(ierr);
650   ierr = TSSetTimeStep(ts,PetscRealConstant(0.001));CHKERRQ(ierr); /* can be overwritten by command line options */
651   ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
652   ierr = TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);CHKERRQ(ierr);
653 
654   ierr = MatZeroEntries(ctx->Jacp);CHKERRQ(ierr);
655   ierr = MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);CHKERRQ(ierr);
656   ierr = MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
657   ierr = MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
658 
659   ierr = TSAdjointSetForward(ts,ctx->Jacp);CHKERRQ(ierr);
660   ierr = VecGetArray(ctx->U,&y_ptr);CHKERRQ(ierr);
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   ierr = VecRestoreArray(ctx->U,&y_ptr);CHKERRQ(ierr);
664   ierr = TSSolve(ts,ctx->U);CHKERRQ(ierr);
665 
666   /* Set terminal conditions for first- and second-order adjonts */
667   ierr = VecGetArrayRead(ctx->U,&z_ptr);CHKERRQ(ierr);
668   ierr = VecGetArray(ctx->Lambda[0],&y_ptr);CHKERRQ(ierr);
669   y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
670   y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
671   ierr = VecRestoreArray(ctx->Lambda[0],&y_ptr);CHKERRQ(ierr);
672   ierr = VecRestoreArrayRead(ctx->U,&z_ptr);CHKERRQ(ierr);
673   ierr = VecGetArray(ctx->Mup[0],&y_ptr);CHKERRQ(ierr);
674   y_ptr[0] = 0.0;
675   ierr = VecRestoreArray(ctx->Mup[0],&y_ptr);CHKERRQ(ierr);
676   ierr = TSForwardGetSensitivities(ts,NULL,&tlmsen);CHKERRQ(ierr);
677   ierr = MatDenseGetColumn(tlmsen,0,&x_ptr);CHKERRQ(ierr);
678   ierr = VecGetArray(ctx->Lambda2[0],&y_ptr);CHKERRQ(ierr);
679   y_ptr[0] = 2.*x_ptr[0];
680   y_ptr[1] = 2.*x_ptr[1];
681   ierr = VecRestoreArray(ctx->Lambda2[0],&y_ptr);CHKERRQ(ierr);
682   ierr = VecGetArray(ctx->Mup2[0],&y_ptr);CHKERRQ(ierr);
683   y_ptr[0] = 0.0;
684   ierr = VecRestoreArray(ctx->Mup2[0],&y_ptr);CHKERRQ(ierr);
685   ierr = MatDenseRestoreColumn(tlmsen,&x_ptr);CHKERRQ(ierr);
686   ierr = TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);CHKERRQ(ierr);
687   if (ctx->implicitform) {
688     ierr = TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);CHKERRQ(ierr);
689   } else {
690     ierr = TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);CHKERRQ(ierr);
691   }
692   ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
693 
694   ierr = VecGetArray(ctx->Lambda[0],&x_ptr);CHKERRQ(ierr);
695   ierr = VecGetArray(ctx->Lambda2[0],&y_ptr);CHKERRQ(ierr);
696   ierr = VecGetArrayRead(ctx->Mup2[0],&z_ptr);CHKERRQ(ierr);
697 
698   arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];
699 
700   ierr = VecRestoreArray(ctx->Lambda2[0],&x_ptr);CHKERRQ(ierr);
701   ierr = VecRestoreArray(ctx->Lambda2[0],&y_ptr);CHKERRQ(ierr);
702   ierr = VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);CHKERRQ(ierr);
703 
704   /* Disable second-order adjoint mode */
705   ierr = TSAdjointReset(ts);CHKERRQ(ierr);
706   ierr = TSAdjointResetForward(ts);CHKERRQ(ierr);
707   PetscFunctionReturn(0);
708 }
709 
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
716 
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
721 
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
726 
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
731 
732 TEST*/
733