Actual source code: ex20opt_p.c

petsc-3.14.6 2021-03-30
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: 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*/