Actual source code: ex20opt_p.c


  2: static char help[] = "Solves the van der Pol equation.\n\
  3: Input parameters include:\n";

  5: /* ------------------------------------------------------------------------

  7:   Notes:
  8:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
  9:   The nonlinear problem is written in a DAE equivalent form.
 10:   The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
 11:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
 12:   ------------------------------------------------------------------------- */
 13: #include <petsctao.h>
 14: #include <petscts.h>

 16: typedef struct _n_User *User;
 17: struct _n_User {
 18:   TS        ts;
 19:   PetscReal mu;
 20:   PetscReal next_output;

 22:   /* Sensitivity analysis support */
 23:   PetscReal ftime;
 24:   Mat       A;                    /* Jacobian matrix */
 25:   Mat       Jacp;                 /* JacobianP matrix */
 26:   Mat       H;                    /* Hessian matrix for optimization */
 27:   Vec       U, Lambda[1], Mup[1]; /* adjoint variables */
 28:   Vec       Lambda2[1], Mup2[1];  /* second-order adjoint variables */
 29:   Vec       Ihp1[1];              /* working space for Hessian evaluations */
 30:   Vec       Ihp2[1];              /* working space for Hessian evaluations */
 31:   Vec       Ihp3[1];              /* working space for Hessian evaluations */
 32:   Vec       Ihp4[1];              /* working space for Hessian evaluations */
 33:   Vec       Dir;                  /* direction vector */
 34:   PetscReal ob[2];                /* observation used by the cost function */
 35:   PetscBool implicitform;         /* implicit ODE? */
 36: };

 38: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
 39: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
 40: PetscErrorCode Adjoint2(Vec, PetscScalar[], User);

 42: /* ----------------------- Explicit form of the ODE  -------------------- */

 44: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
 45: {
 46:   User               user = (User)ctx;
 47:   PetscScalar       *f;
 48:   const PetscScalar *u;

 50:   PetscFunctionBeginUser;
 51:   PetscCall(VecGetArrayRead(U, &u));
 52:   PetscCall(VecGetArray(F, &f));
 53:   f[0] = u[1];
 54:   f[1] = user->mu * ((1. - u[0] * u[0]) * u[1] - u[0]);
 55:   PetscCall(VecRestoreArrayRead(U, &u));
 56:   PetscCall(VecRestoreArray(F, &f));
 57:   PetscFunctionReturn(PETSC_SUCCESS);
 58: }

 60: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
 61: {
 62:   User               user     = (User)ctx;
 63:   PetscReal          mu       = user->mu;
 64:   PetscInt           rowcol[] = {0, 1};
 65:   PetscScalar        J[2][2];
 66:   const PetscScalar *u;

 68:   PetscFunctionBeginUser;
 69:   PetscCall(VecGetArrayRead(U, &u));

 71:   J[0][0] = 0;
 72:   J[1][0] = -mu * (2.0 * u[1] * u[0] + 1.);
 73:   J[0][1] = 1.0;
 74:   J[1][1] = mu * (1.0 - u[0] * u[0]);
 75:   PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
 76:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 77:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 78:   if (B && A != B) {
 79:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 80:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 81:   }

 83:   PetscCall(VecRestoreArrayRead(U, &u));
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

 87: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
 88: {
 89:   const PetscScalar *vl, *vr, *u;
 90:   PetscScalar       *vhv;
 91:   PetscScalar        dJdU[2][2][2] = {{{0}}};
 92:   PetscInt           i, j, k;
 93:   User               user = (User)ctx;

 95:   PetscFunctionBeginUser;
 96:   PetscCall(VecGetArrayRead(U, &u));
 97:   PetscCall(VecGetArrayRead(Vl[0], &vl));
 98:   PetscCall(VecGetArrayRead(Vr, &vr));
 99:   PetscCall(VecGetArray(VHV[0], &vhv));

101:   dJdU[1][0][0] = -2. * user->mu * u[1];
102:   dJdU[1][1][0] = -2. * user->mu * u[0];
103:   dJdU[1][0][1] = -2. * user->mu * u[0];
104:   for (j = 0; j < 2; j++) {
105:     vhv[j] = 0;
106:     for (k = 0; k < 2; k++)
107:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
108:   }

110:   PetscCall(VecRestoreArrayRead(U, &u));
111:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
112:   PetscCall(VecRestoreArrayRead(Vr, &vr));
113:   PetscCall(VecRestoreArray(VHV[0], &vhv));
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: static PetscErrorCode RHSHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
118: {
119:   const PetscScalar *vl, *vr, *u;
120:   PetscScalar       *vhv;
121:   PetscScalar        dJdP[2][2][1] = {{{0}}};
122:   PetscInt           i, j, k;

124:   PetscFunctionBeginUser;
125:   PetscCall(VecGetArrayRead(U, &u));
126:   PetscCall(VecGetArrayRead(Vl[0], &vl));
127:   PetscCall(VecGetArrayRead(Vr, &vr));
128:   PetscCall(VecGetArray(VHV[0], &vhv));

130:   dJdP[1][0][0] = -(1. + 2. * u[0] * u[1]);
131:   dJdP[1][1][0] = 1. - u[0] * u[0];
132:   for (j = 0; j < 2; j++) {
133:     vhv[j] = 0;
134:     for (k = 0; k < 1; k++)
135:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
136:   }

138:   PetscCall(VecRestoreArrayRead(U, &u));
139:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
140:   PetscCall(VecRestoreArrayRead(Vr, &vr));
141:   PetscCall(VecRestoreArray(VHV[0], &vhv));
142:   PetscFunctionReturn(PETSC_SUCCESS);
143: }

145: static PetscErrorCode RHSHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
146: {
147:   const PetscScalar *vl, *vr, *u;
148:   PetscScalar       *vhv;
149:   PetscScalar        dJdU[2][1][2] = {{{0}}};
150:   PetscInt           i, j, k;

152:   PetscFunctionBeginUser;
153:   PetscCall(VecGetArrayRead(U, &u));
154:   PetscCall(VecGetArrayRead(Vl[0], &vl));
155:   PetscCall(VecGetArrayRead(Vr, &vr));
156:   PetscCall(VecGetArray(VHV[0], &vhv));

158:   dJdU[1][0][0] = -1. - 2. * u[1] * u[0];
159:   dJdU[1][0][1] = 1. - u[0] * u[0];
160:   for (j = 0; j < 1; j++) {
161:     vhv[j] = 0;
162:     for (k = 0; k < 2; k++)
163:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
164:   }

166:   PetscCall(VecRestoreArrayRead(U, &u));
167:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
168:   PetscCall(VecRestoreArrayRead(Vr, &vr));
169:   PetscCall(VecRestoreArray(VHV[0], &vhv));
170:   PetscFunctionReturn(PETSC_SUCCESS);
171: }

173: static PetscErrorCode RHSHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
174: {
175:   PetscFunctionBeginUser;
176:   PetscFunctionReturn(PETSC_SUCCESS);
177: }

179: /* ----------------------- Implicit form of the ODE  -------------------- */

181: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
182: {
183:   User               user = (User)ctx;
184:   PetscScalar       *f;
185:   const PetscScalar *u, *udot;

187:   PetscFunctionBeginUser;
188:   PetscCall(VecGetArrayRead(U, &u));
189:   PetscCall(VecGetArrayRead(Udot, &udot));
190:   PetscCall(VecGetArray(F, &f));

192:   f[0] = udot[0] - u[1];
193:   f[1] = udot[1] - user->mu * ((1.0 - u[0] * u[0]) * u[1] - u[0]);

195:   PetscCall(VecRestoreArrayRead(U, &u));
196:   PetscCall(VecRestoreArrayRead(Udot, &udot));
197:   PetscCall(VecRestoreArray(F, &f));
198:   PetscFunctionReturn(PETSC_SUCCESS);
199: }

201: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
202: {
203:   User               user     = (User)ctx;
204:   PetscInt           rowcol[] = {0, 1};
205:   PetscScalar        J[2][2];
206:   const PetscScalar *u;

208:   PetscFunctionBeginUser;
209:   PetscCall(VecGetArrayRead(U, &u));

211:   J[0][0] = a;
212:   J[0][1] = -1.0;
213:   J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
214:   J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
215:   PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
216:   PetscCall(VecRestoreArrayRead(U, &u));
217:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
218:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
219:   if (A != B) {
220:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
221:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
222:   }

224:   PetscCall(VecRestoreArrayRead(U, &u));
225:   PetscFunctionReturn(PETSC_SUCCESS);
226: }

228: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec U, Mat A, void *ctx)
229: {
230:   PetscInt           row[] = {0, 1}, col[] = {0};
231:   PetscScalar        J[2][1];
232:   const PetscScalar *u;

234:   PetscFunctionBeginUser;
235:   PetscCall(VecGetArrayRead(U, &u));

237:   J[0][0] = 0;
238:   J[1][0] = (1. - u[0] * u[0]) * u[1] - u[0];
239:   PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
240:   PetscCall(VecRestoreArrayRead(U, &u));
241:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
242:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

244:   PetscCall(VecRestoreArrayRead(U, &u));
245:   PetscFunctionReturn(PETSC_SUCCESS);
246: }

248: static PetscErrorCode IHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
249: {
250:   const PetscScalar *vl, *vr, *u;
251:   PetscScalar       *vhv;
252:   PetscScalar        dJdU[2][2][2] = {{{0}}};
253:   PetscInt           i, j, k;
254:   User               user = (User)ctx;

256:   PetscFunctionBeginUser;
257:   PetscCall(VecGetArrayRead(U, &u));
258:   PetscCall(VecGetArrayRead(Vl[0], &vl));
259:   PetscCall(VecGetArrayRead(Vr, &vr));
260:   PetscCall(VecGetArray(VHV[0], &vhv));

262:   dJdU[1][0][0] = 2. * user->mu * u[1];
263:   dJdU[1][1][0] = 2. * user->mu * u[0];
264:   dJdU[1][0][1] = 2. * user->mu * u[0];
265:   for (j = 0; j < 2; j++) {
266:     vhv[j] = 0;
267:     for (k = 0; k < 2; k++)
268:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
269:   }

271:   PetscCall(VecRestoreArrayRead(U, &u));
272:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
273:   PetscCall(VecRestoreArrayRead(Vr, &vr));
274:   PetscCall(VecRestoreArray(VHV[0], &vhv));
275:   PetscFunctionReturn(PETSC_SUCCESS);
276: }

278: static PetscErrorCode IHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
279: {
280:   const PetscScalar *vl, *vr, *u;
281:   PetscScalar       *vhv;
282:   PetscScalar        dJdP[2][2][1] = {{{0}}};
283:   PetscInt           i, j, k;

285:   PetscFunctionBeginUser;
286:   PetscCall(VecGetArrayRead(U, &u));
287:   PetscCall(VecGetArrayRead(Vl[0], &vl));
288:   PetscCall(VecGetArrayRead(Vr, &vr));
289:   PetscCall(VecGetArray(VHV[0], &vhv));

291:   dJdP[1][0][0] = 1. + 2. * u[0] * u[1];
292:   dJdP[1][1][0] = u[0] * u[0] - 1.;
293:   for (j = 0; j < 2; j++) {
294:     vhv[j] = 0;
295:     for (k = 0; k < 1; k++)
296:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
297:   }

299:   PetscCall(VecRestoreArrayRead(U, &u));
300:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
301:   PetscCall(VecRestoreArrayRead(Vr, &vr));
302:   PetscCall(VecRestoreArray(VHV[0], &vhv));
303:   PetscFunctionReturn(PETSC_SUCCESS);
304: }

306: static PetscErrorCode IHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
307: {
308:   const PetscScalar *vl, *vr, *u;
309:   PetscScalar       *vhv;
310:   PetscScalar        dJdU[2][1][2] = {{{0}}};
311:   PetscInt           i, j, k;

313:   PetscFunctionBeginUser;
314:   PetscCall(VecGetArrayRead(U, &u));
315:   PetscCall(VecGetArrayRead(Vl[0], &vl));
316:   PetscCall(VecGetArrayRead(Vr, &vr));
317:   PetscCall(VecGetArray(VHV[0], &vhv));

319:   dJdU[1][0][0] = 1. + 2. * u[1] * u[0];
320:   dJdU[1][0][1] = u[0] * u[0] - 1.;
321:   for (j = 0; j < 1; j++) {
322:     vhv[j] = 0;
323:     for (k = 0; k < 2; k++)
324:       for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
325:   }

327:   PetscCall(VecRestoreArrayRead(U, &u));
328:   PetscCall(VecRestoreArrayRead(Vl[0], &vl));
329:   PetscCall(VecRestoreArrayRead(Vr, &vr));
330:   PetscCall(VecRestoreArray(VHV[0], &vhv));
331:   PetscFunctionReturn(PETSC_SUCCESS);
332: }

334: static PetscErrorCode IHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
335: {
336:   PetscFunctionBeginUser;
337:   PetscFunctionReturn(PETSC_SUCCESS);
338: }

340: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
341: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
342: {
343:   const PetscScalar *x;
344:   PetscReal          tfinal, dt;
345:   User               user = (User)ctx;
346:   Vec                interpolatedX;

348:   PetscFunctionBeginUser;
349:   PetscCall(TSGetTimeStep(ts, &dt));
350:   PetscCall(TSGetMaxTime(ts, &tfinal));

352:   while (user->next_output <= t && user->next_output <= tfinal) {
353:     PetscCall(VecDuplicate(X, &interpolatedX));
354:     PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
355:     PetscCall(VecGetArrayRead(interpolatedX, &x));
356:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%g] %" PetscInt_FMT " TS %g (dt = %g) X %g %g\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1])));
357:     PetscCall(VecRestoreArrayRead(interpolatedX, &x));
358:     PetscCall(VecDestroy(&interpolatedX));
359:     user->next_output += PetscRealConstant(0.1);
360:   }
361:   PetscFunctionReturn(PETSC_SUCCESS);
362: }

364: int main(int argc, char **argv)
365: {
366:   Vec                P;
367:   PetscBool          monitor = PETSC_FALSE;
368:   PetscScalar       *x_ptr;
369:   const PetscScalar *y_ptr;
370:   PetscMPIInt        size;
371:   struct _n_User     user;
372:   Tao                tao;
373:   KSP                ksp;
374:   PC                 pc;

376:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
377:      Initialize program
378:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
379:   PetscFunctionBeginUser;
380:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
381:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
382:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");

384:   /* Create TAO solver and set desired solution method */
385:   PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
386:   PetscCall(TaoSetType(tao, TAOBQNLS));

388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389:     Set runtime options
390:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391:   user.next_output  = 0.0;
392:   user.mu           = PetscRealConstant(1.0e3);
393:   user.ftime        = PetscRealConstant(0.5);
394:   user.implicitform = PETSC_TRUE;

396:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
397:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
398:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &user.implicitform, NULL));

400:   /* Create necessary matrix and vectors, solve same ODE on every process */
401:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.A));
402:   PetscCall(MatSetSizes(user.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
403:   PetscCall(MatSetFromOptions(user.A));
404:   PetscCall(MatSetUp(user.A));
405:   PetscCall(MatCreateVecs(user.A, &user.U, NULL));
406:   PetscCall(MatCreateVecs(user.A, &user.Lambda[0], NULL));
407:   PetscCall(MatCreateVecs(user.A, &user.Lambda2[0], NULL));
408:   PetscCall(MatCreateVecs(user.A, &user.Ihp1[0], NULL));
409:   PetscCall(MatCreateVecs(user.A, &user.Ihp2[0], NULL));

411:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
412:   PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
413:   PetscCall(MatSetFromOptions(user.Jacp));
414:   PetscCall(MatSetUp(user.Jacp));
415:   PetscCall(MatCreateVecs(user.Jacp, &user.Dir, NULL));
416:   PetscCall(MatCreateVecs(user.Jacp, &user.Mup[0], NULL));
417:   PetscCall(MatCreateVecs(user.Jacp, &user.Mup2[0], NULL));
418:   PetscCall(MatCreateVecs(user.Jacp, &user.Ihp3[0], NULL));
419:   PetscCall(MatCreateVecs(user.Jacp, &user.Ihp4[0], NULL));

421:   /* Create timestepping solver context */
422:   PetscCall(TSCreate(PETSC_COMM_WORLD, &user.ts));
423:   PetscCall(TSSetEquationType(user.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
424:   if (user.implicitform) {
425:     PetscCall(TSSetIFunction(user.ts, NULL, IFunction, &user));
426:     PetscCall(TSSetIJacobian(user.ts, user.A, user.A, IJacobian, &user));
427:     PetscCall(TSSetType(user.ts, TSCN));
428:   } else {
429:     PetscCall(TSSetRHSFunction(user.ts, NULL, RHSFunction, &user));
430:     PetscCall(TSSetRHSJacobian(user.ts, user.A, user.A, RHSJacobian, &user));
431:     PetscCall(TSSetType(user.ts, TSRK));
432:   }
433:   PetscCall(TSSetRHSJacobianP(user.ts, user.Jacp, RHSJacobianP, &user));
434:   PetscCall(TSSetMaxTime(user.ts, user.ftime));
435:   PetscCall(TSSetExactFinalTime(user.ts, TS_EXACTFINALTIME_MATCHSTEP));

437:   if (monitor) PetscCall(TSMonitorSet(user.ts, Monitor, &user, NULL));

439:   /* Set ODE initial conditions */
440:   PetscCall(VecGetArray(user.U, &x_ptr));
441:   x_ptr[0] = 2.0;
442:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
443:   PetscCall(VecRestoreArray(user.U, &x_ptr));
444:   PetscCall(TSSetTimeStep(user.ts, PetscRealConstant(0.001)));

446:   /* Set runtime options */
447:   PetscCall(TSSetFromOptions(user.ts));

449:   PetscCall(TSSolve(user.ts, user.U));
450:   PetscCall(VecGetArrayRead(user.U, &y_ptr));
451:   user.ob[0] = y_ptr[0];
452:   user.ob[1] = y_ptr[1];
453:   PetscCall(VecRestoreArrayRead(user.U, &y_ptr));

455:   /* Save trajectory of solution so that TSAdjointSolve() may be used.
456:      Skip checkpointing for the first TSSolve since no adjoint run follows it.
457:    */
458:   PetscCall(TSSetSaveTrajectory(user.ts));

460:   /* Optimization starts */
461:   PetscCall(MatCreate(PETSC_COMM_WORLD, &user.H));
462:   PetscCall(MatSetSizes(user.H, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
463:   PetscCall(MatSetUp(user.H)); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */

465:   /* Set initial solution guess */
466:   PetscCall(MatCreateVecs(user.Jacp, &P, NULL));
467:   PetscCall(VecGetArray(P, &x_ptr));
468:   x_ptr[0] = PetscRealConstant(1.2);
469:   PetscCall(VecRestoreArray(P, &x_ptr));
470:   PetscCall(TaoSetSolution(tao, P));

472:   /* Set routine for function and gradient evaluation */
473:   PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
474:   PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));

476:   /* Check for any TAO command line options */
477:   PetscCall(TaoGetKSP(tao, &ksp));
478:   if (ksp) {
479:     PetscCall(KSPGetPC(ksp, &pc));
480:     PetscCall(PCSetType(pc, PCNONE));
481:   }
482:   PetscCall(TaoSetFromOptions(tao));

484:   PetscCall(TaoSolve(tao));

486:   PetscCall(VecView(P, PETSC_VIEWER_STDOUT_WORLD));
487:   /* Free TAO data structures */
488:   PetscCall(TaoDestroy(&tao));

490:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
491:      Free work space.  All PETSc objects should be destroyed when they
492:      are no longer needed.
493:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
494:   PetscCall(MatDestroy(&user.H));
495:   PetscCall(MatDestroy(&user.A));
496:   PetscCall(VecDestroy(&user.U));
497:   PetscCall(MatDestroy(&user.Jacp));
498:   PetscCall(VecDestroy(&user.Lambda[0]));
499:   PetscCall(VecDestroy(&user.Mup[0]));
500:   PetscCall(VecDestroy(&user.Lambda2[0]));
501:   PetscCall(VecDestroy(&user.Mup2[0]));
502:   PetscCall(VecDestroy(&user.Ihp1[0]));
503:   PetscCall(VecDestroy(&user.Ihp2[0]));
504:   PetscCall(VecDestroy(&user.Ihp3[0]));
505:   PetscCall(VecDestroy(&user.Ihp4[0]));
506:   PetscCall(VecDestroy(&user.Dir));
507:   PetscCall(TSDestroy(&user.ts));
508:   PetscCall(VecDestroy(&P));
509:   PetscCall(PetscFinalize());
510:   return 0;
511: }

513: /* ------------------------------------------------------------------ */
514: /*
515:    FormFunctionGradient - Evaluates the function and corresponding gradient.

517:    Input Parameters:
518:    tao - the Tao context
519:    X   - the input vector
520:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()

522:    Output Parameters:
523:    f   - the newly evaluated function
524:    G   - the newly evaluated gradient
525: */
526: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx)
527: {
528:   User               user_ptr = (User)ctx;
529:   TS                 ts       = user_ptr->ts;
530:   PetscScalar       *x_ptr, *g;
531:   const PetscScalar *y_ptr;

533:   PetscFunctionBeginUser;
534:   PetscCall(VecGetArrayRead(P, &y_ptr));
535:   user_ptr->mu = y_ptr[0];
536:   PetscCall(VecRestoreArrayRead(P, &y_ptr));

538:   PetscCall(TSSetTime(ts, 0.0));
539:   PetscCall(TSSetStepNumber(ts, 0));
540:   PetscCall(TSSetTimeStep(ts, PetscRealConstant(0.001))); /* can be overwritten by command line options */
541:   PetscCall(TSSetFromOptions(ts));
542:   PetscCall(VecGetArray(user_ptr->U, &x_ptr));
543:   x_ptr[0] = 2.0;
544:   x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user_ptr->mu) - 292.0 / (2187.0 * user_ptr->mu * user_ptr->mu);
545:   PetscCall(VecRestoreArray(user_ptr->U, &x_ptr));

547:   PetscCall(TSSolve(ts, user_ptr->U));

549:   PetscCall(VecGetArrayRead(user_ptr->U, &y_ptr));
550:   *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]);

552:   /*   Reset initial conditions for the adjoint integration */
553:   PetscCall(VecGetArray(user_ptr->Lambda[0], &x_ptr));
554:   x_ptr[0] = 2. * (y_ptr[0] - user_ptr->ob[0]);
555:   x_ptr[1] = 2. * (y_ptr[1] - user_ptr->ob[1]);
556:   PetscCall(VecRestoreArrayRead(user_ptr->U, &y_ptr));
557:   PetscCall(VecRestoreArray(user_ptr->Lambda[0], &x_ptr));

559:   PetscCall(VecGetArray(user_ptr->Mup[0], &x_ptr));
560:   x_ptr[0] = 0.0;
561:   PetscCall(VecRestoreArray(user_ptr->Mup[0], &x_ptr));
562:   PetscCall(TSSetCostGradients(ts, 1, user_ptr->Lambda, user_ptr->Mup));

564:   PetscCall(TSAdjointSolve(ts));

566:   PetscCall(VecGetArray(user_ptr->Mup[0], &x_ptr));
567:   PetscCall(VecGetArrayRead(user_ptr->Lambda[0], &y_ptr));
568:   PetscCall(VecGetArray(G, &g));
569:   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];
570:   PetscCall(VecRestoreArray(user_ptr->Mup[0], &x_ptr));
571:   PetscCall(VecRestoreArrayRead(user_ptr->Lambda[0], &y_ptr));
572:   PetscCall(VecRestoreArray(G, &g));
573:   PetscFunctionReturn(PETSC_SUCCESS);
574: }

576: PetscErrorCode FormHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
577: {
578:   User           user_ptr = (User)ctx;
579:   PetscScalar    harr[1];
580:   const PetscInt rows[1] = {0};
581:   PetscInt       col     = 0;

583:   PetscFunctionBeginUser;
584:   PetscCall(Adjoint2(P, harr, user_ptr));
585:   PetscCall(MatSetValues(H, 1, rows, 1, &col, harr, INSERT_VALUES));

587:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
588:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
589:   if (H != Hpre) {
590:     PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
591:     PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
592:   }
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: PetscErrorCode Adjoint2(Vec P, PetscScalar arr[], User ctx)
597: {
598:   TS                 ts = ctx->ts;
599:   const PetscScalar *z_ptr;
600:   PetscScalar       *x_ptr, *y_ptr, dzdp, dzdp2;
601:   Mat                tlmsen;

603:   PetscFunctionBeginUser;
604:   /* Reset TSAdjoint so that AdjointSetUp will be called again */
605:   PetscCall(TSAdjointReset(ts));

607:   /* The directional vector should be 1 since it is one-dimensional */
608:   PetscCall(VecGetArray(ctx->Dir, &x_ptr));
609:   x_ptr[0] = 1.;
610:   PetscCall(VecRestoreArray(ctx->Dir, &x_ptr));

612:   PetscCall(VecGetArrayRead(P, &z_ptr));
613:   ctx->mu = z_ptr[0];
614:   PetscCall(VecRestoreArrayRead(P, &z_ptr));

616:   dzdp  = -10.0 / (81.0 * ctx->mu * ctx->mu) + 2.0 * 292.0 / (2187.0 * ctx->mu * ctx->mu * ctx->mu);
617:   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);

619:   PetscCall(TSSetTime(ts, 0.0));
620:   PetscCall(TSSetStepNumber(ts, 0));
621:   PetscCall(TSSetTimeStep(ts, PetscRealConstant(0.001))); /* can be overwritten by command line options */
622:   PetscCall(TSSetFromOptions(ts));
623:   PetscCall(TSSetCostHessianProducts(ts, 1, ctx->Lambda2, ctx->Mup2, ctx->Dir));

625:   PetscCall(MatZeroEntries(ctx->Jacp));
626:   PetscCall(MatSetValue(ctx->Jacp, 1, 0, dzdp, INSERT_VALUES));
627:   PetscCall(MatAssemblyBegin(ctx->Jacp, MAT_FINAL_ASSEMBLY));
628:   PetscCall(MatAssemblyEnd(ctx->Jacp, MAT_FINAL_ASSEMBLY));

630:   PetscCall(TSAdjointSetForward(ts, ctx->Jacp));
631:   PetscCall(VecGetArray(ctx->U, &y_ptr));
632:   y_ptr[0] = 2.0;
633:   y_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * ctx->mu) - 292.0 / (2187.0 * ctx->mu * ctx->mu);
634:   PetscCall(VecRestoreArray(ctx->U, &y_ptr));
635:   PetscCall(TSSolve(ts, ctx->U));

637:   /* Set terminal conditions for first- and second-order adjonts */
638:   PetscCall(VecGetArrayRead(ctx->U, &z_ptr));
639:   PetscCall(VecGetArray(ctx->Lambda[0], &y_ptr));
640:   y_ptr[0] = 2. * (z_ptr[0] - ctx->ob[0]);
641:   y_ptr[1] = 2. * (z_ptr[1] - ctx->ob[1]);
642:   PetscCall(VecRestoreArray(ctx->Lambda[0], &y_ptr));
643:   PetscCall(VecRestoreArrayRead(ctx->U, &z_ptr));
644:   PetscCall(VecGetArray(ctx->Mup[0], &y_ptr));
645:   y_ptr[0] = 0.0;
646:   PetscCall(VecRestoreArray(ctx->Mup[0], &y_ptr));
647:   PetscCall(TSForwardGetSensitivities(ts, NULL, &tlmsen));
648:   PetscCall(MatDenseGetColumn(tlmsen, 0, &x_ptr));
649:   PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
650:   y_ptr[0] = 2. * x_ptr[0];
651:   y_ptr[1] = 2. * x_ptr[1];
652:   PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
653:   PetscCall(VecGetArray(ctx->Mup2[0], &y_ptr));
654:   y_ptr[0] = 0.0;
655:   PetscCall(VecRestoreArray(ctx->Mup2[0], &y_ptr));
656:   PetscCall(MatDenseRestoreColumn(tlmsen, &x_ptr));
657:   PetscCall(TSSetCostGradients(ts, 1, ctx->Lambda, ctx->Mup));
658:   if (ctx->implicitform) {
659:     PetscCall(TSSetIHessianProduct(ts, ctx->Ihp1, IHessianProductUU, ctx->Ihp2, IHessianProductUP, ctx->Ihp3, IHessianProductPU, ctx->Ihp4, IHessianProductPP, ctx));
660:   } else {
661:     PetscCall(TSSetRHSHessianProduct(ts, ctx->Ihp1, RHSHessianProductUU, ctx->Ihp2, RHSHessianProductUP, ctx->Ihp3, RHSHessianProductPU, ctx->Ihp4, RHSHessianProductPP, ctx));
662:   }
663:   PetscCall(TSAdjointSolve(ts));

665:   PetscCall(VecGetArray(ctx->Lambda[0], &x_ptr));
666:   PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
667:   PetscCall(VecGetArrayRead(ctx->Mup2[0], &z_ptr));

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

671:   PetscCall(VecRestoreArray(ctx->Lambda2[0], &x_ptr));
672:   PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
673:   PetscCall(VecRestoreArrayRead(ctx->Mup2[0], &z_ptr));

675:   /* Disable second-order adjoint mode */
676:   PetscCall(TSAdjointReset(ts));
677:   PetscCall(TSAdjointResetForward(ts));
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: /*TEST
682:     build:
683:       requires: !complex !single
684:     test:
685:       args:  -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
686:       output_file: output/ex20opt_p_1.out

688:     test:
689:       suffix: 2
690:       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
691:       output_file: output/ex20opt_p_2.out

693:     test:
694:       suffix: 3
695:       args:  -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
696:       output_file: output/ex20opt_p_3.out

698:     test:
699:       suffix: 4
700:       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
701:       output_file: output/ex20opt_p_4.out

703: TEST*/