Actual source code: ex20opt_p.c

petsc-master 2019-05-18
Report Typos and Errors

  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: 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<2; k++)
146:       for (i=0; i<1; 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<2; k++)
314:       for (i=0; i<1; 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_SELF,1,"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:     if (user.implicitform) {
448:     TSSetIFunction(user.ts,NULL,IFunction,&user);
449:     TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user);
450:     TSSetType(user.ts,TSCN);
451:   } else {
452:     TSSetRHSFunction(user.ts,NULL,RHSFunction,&user);
453:     TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user);
454:     TSSetType(user.ts,TSRK);
455:   }
456:   TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user);
457:   TSSetMaxTime(user.ts,user.ftime);
458:   TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP);

460:   if (monitor) {
461:     TSMonitorSet(user.ts,Monitor,&user,NULL);
462:   }

464:   /* Set ODE initial conditions */
465:   VecGetArray(user.U,&x_ptr);
466:   x_ptr[0] = 2.0;
467:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu);
468:   VecRestoreArray(user.U,&x_ptr);
469:   TSSetTimeStep(user.ts,PetscRealConstant(0.001));

471:   /* Set runtime options */
472:   TSSetFromOptions(user.ts);

474:   TSSolve(user.ts,user.U);
475:   VecGetArrayRead(user.U,&y_ptr);
476:   user.ob[0] = y_ptr[0];
477:   user.ob[1] = y_ptr[1];
478:   VecRestoreArrayRead(user.U,&y_ptr);

480:   /* Save trajectory of solution so that TSAdjointSolve() may be used.
481:      Skip checkpointing for the first TSSolve since no adjoint run follows it.
482:    */
483:   TSSetSaveTrajectory(user.ts);

485:   /* Optimization starts */
486:   MatCreate(PETSC_COMM_WORLD,&user.H);
487:   MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1);
488:   MatSetUp(user.H); /* Hessian should be symmetic. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */

490:   /* Set initial solution guess */
491:   MatCreateVecs(user.Jacp,&P,NULL);
492:   VecGetArray(P,&x_ptr);
493:   x_ptr[0] = PetscRealConstant(1.2);
494:   VecRestoreArray(P,&x_ptr);
495:   TaoSetInitialVector(tao,P);

497:   /* Set routine for function and gradient evaluation */
498:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
499:   TaoSetHessianRoutine(tao,user.H,user.H,FormHessian,(void *)&user);

501:   /* Check for any TAO command line options */
502:   TaoGetKSP(tao,&ksp);
503:   if (ksp) {
504:     KSPGetPC(ksp,&pc);
505:     PCSetType(pc,PCNONE);
506:   }
507:   TaoSetFromOptions(tao);

509:   TaoSolve(tao);

511:   VecView(P,PETSC_VIEWER_STDOUT_WORLD);
512:   /* Free TAO data structures */
513:   TaoDestroy(&tao);

515:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
516:      Free work space.  All PETSc objects should be destroyed when they
517:      are no longer needed.
518:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
519:   MatDestroy(&user.H);
520:   MatDestroy(&user.A);
521:   VecDestroy(&user.U);
522:   MatDestroy(&user.Jacp);
523:   VecDestroy(&user.Lambda[0]);
524:   VecDestroy(&user.Mup[0]);
525:   VecDestroy(&user.Lambda2[0]);
526:   VecDestroy(&user.Mup2[0]);
527:   VecDestroy(&user.Ihp1[0]);
528:   VecDestroy(&user.Ihp2[0]);
529:   VecDestroy(&user.Ihp3[0]);
530:   VecDestroy(&user.Ihp4[0]);
531:   VecDestroy(&user.Dir);
532:   TSDestroy(&user.ts);
533:   VecDestroy(&P);
534:   PetscFinalize();
535:   return ierr;
536: }

538: /* ------------------------------------------------------------------ */
539: /*
540:    FormFunctionGradient - Evaluates the function and corresponding gradient.

542:    Input Parameters:
543:    tao - the Tao context
544:    X   - the input vector
545:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

547:    Output Parameters:
548:    f   - the newly evaluated function
549:    G   - the newly evaluated gradient
550: */
551: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
552: {
553:   User              user_ptr = (User)ctx;
554:   TS                ts = user_ptr->ts;
555:   PetscScalar       *x_ptr,*g;
556:   const PetscScalar *y_ptr;
557:   PetscErrorCode    ierr;

560:   VecGetArrayRead(P,&y_ptr);
561:   user_ptr->mu = y_ptr[0];
562:   VecRestoreArrayRead(P,&y_ptr);

564:   TSSetTime(ts,0.0);
565:   TSSetStepNumber(ts,0);
566:   TSSetTimeStep(ts,PetscRealConstant(0.001)); /* can be overwritten by command line options */
567:   TSSetFromOptions(ts);
568:   VecGetArray(user_ptr->U,&x_ptr);
569:   x_ptr[0] = 2.0;
570:   x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu);
571:   VecRestoreArray(user_ptr->U,&x_ptr);

573:   TSSolve(ts,user_ptr->U);

575:   VecGetArrayRead(user_ptr->U,&y_ptr);
576:   *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:   /*   Reset initial conditions for the adjoint integration */
579:   VecGetArray(user_ptr->Lambda[0],&x_ptr);
580:   x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]);
581:   x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]);
582:   VecRestoreArrayRead(user_ptr->U,&y_ptr);
583:   VecRestoreArray(user_ptr->Lambda[0],&x_ptr);

585:   VecGetArray(user_ptr->Mup[0],&x_ptr);
586:   x_ptr[0] = 0.0;
587:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
588:   TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup);

590:   TSAdjointSolve(ts);

592:   VecGetArray(user_ptr->Mup[0],&x_ptr);
593:   VecGetArrayRead(user_ptr->Lambda[0],&y_ptr);
594:   VecGetArray(G,&g);
595:   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];
596:   VecRestoreArray(user_ptr->Mup[0],&x_ptr);
597:   VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr);
598:   VecRestoreArray(G,&g);
599:   return(0);
600: }

602: PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
603: {
604:   User           user_ptr = (User)ctx;
605:   PetscScalar    harr[1];
606:   const PetscInt rows[1] = {0};
607:   PetscInt       col = 0;

611:   Adjoint2(P,harr,user_ptr);
612:   MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES);

614:   MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
615:   MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
616:   if (H != Hpre) {
617:     MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
618:     MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
619:   }
620:   return(0);
621: }

623: PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx)
624: {
625:   TS                ts = ctx->ts;
626:   const PetscScalar *z_ptr;
627:   PetscScalar       *x_ptr,*y_ptr,dzdp,dzdp2;
628:   Mat               tlmsen;
629:   PetscErrorCode    ierr;

632:   /* Reset TSAdjoint so that AdjointSetUp will be called again */
633:   TSAdjointReset(ts);

635:   /* The directional vector should be 1 since it is one-dimensional */
636:   VecGetArray(ctx->Dir,&x_ptr);
637:   x_ptr[0] = 1.;
638:   VecRestoreArray(ctx->Dir,&x_ptr);

640:   VecGetArrayRead(P,&z_ptr);
641:   ctx->mu = z_ptr[0];
642:   VecRestoreArrayRead(P,&z_ptr);

644:   dzdp  = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu);
645:   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:   TSSetTime(ts,0.0);
648:   TSSetStepNumber(ts,0);
649:   TSSetTimeStep(ts,0.001);
650:   TSSetFromOptions(ts);
651:   TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir);

653:   MatZeroEntries(ctx->Jacp);
654:   MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES);
655:   MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY);
656:   MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY);

658:   TSAdjointSetForward(ts,ctx->Jacp);
659:   VecGetArray(ctx->U,&y_ptr);
660:   y_ptr[0] = 2.0;
661:   y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu);
662:   VecRestoreArray(ctx->U,&y_ptr);
663:   TSSolve(ts,ctx->U);

665:   /* Set terminal conditions for first- and second-order adjonts */
666:   VecGetArrayRead(ctx->U,&z_ptr);
667:   VecGetArray(ctx->Lambda[0],&y_ptr);
668:   y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]);
669:   y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]);
670:   VecRestoreArray(ctx->Lambda[0],&y_ptr);
671:   VecRestoreArrayRead(ctx->U,&z_ptr);
672:   VecGetArray(ctx->Mup[0],&y_ptr);
673:   y_ptr[0] = 0.0;
674:   VecRestoreArray(ctx->Mup[0],&y_ptr);
675:   TSForwardGetSensitivities(ts,NULL,&tlmsen);
676:   MatDenseGetColumn(tlmsen,0,&x_ptr);
677:   VecGetArray(ctx->Lambda2[0],&y_ptr);
678:   y_ptr[0] = 2.*x_ptr[0];
679:   y_ptr[1] = 2.*x_ptr[1];
680:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
681:   VecGetArray(ctx->Mup2[0],&y_ptr);
682:   y_ptr[0] = 0.0;
683:   VecRestoreArray(ctx->Mup2[0],&y_ptr);
684:   MatDenseRestoreColumn(tlmsen,&x_ptr);
685:   TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup);
686:   if (ctx->implicitform) {
687:     TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx);
688:   } else {
689:     TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx);
690:   }
691:   TSAdjointSolve(ts);

693:   VecGetArray(ctx->Lambda[0],&x_ptr);
694:   VecGetArray(ctx->Lambda2[0],&y_ptr);
695:   VecGetArrayRead(ctx->Mup2[0],&z_ptr);

697:   arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0];

699:   VecRestoreArray(ctx->Lambda2[0],&x_ptr);
700:   VecRestoreArray(ctx->Lambda2[0],&y_ptr);
701:   VecRestoreArrayRead(ctx->Mup2[0],&z_ptr);

703:   /* Disable second-order adjoint mode */
704:   TSAdjointReset(ts);
705:   TSAdjointResetForward(ts);
706:   return(0);
707: }

709: /*TEST
710:     build:
711:       requires: !complex !single
712:     test:
713:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
714:       output_file: output/ex20opt_p_1.out

716:     test:
717:       suffix: 2
718:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none
719:       output_file: output/ex20opt_p_2.out

721:     test:
722:       suffix: 3
723:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01
724:       output_file: output/ex20opt_p_3.out

726:     test:
727:       suffix: 4
728:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none
729:       output_file: output/ex20opt_p_4.out

731: TEST*/