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*/