Actual source code: ex20opt_p.c
1: static char help[] = "Solves the van der Pol equation.\n\
2: Input parameters include:\n";
4: /* ------------------------------------------------------------------------
6: Notes:
7: This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
8: The nonlinear problem is written in a DAE equivalent form.
9: The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
10: The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
11: ------------------------------------------------------------------------- */
12: #include <petsctao.h>
13: #include <petscts.h>
15: typedef struct _n_User *User;
16: struct _n_User {
17: TS ts;
18: PetscReal mu;
19: PetscReal next_output;
21: /* Sensitivity analysis support */
22: PetscReal ftime;
23: Mat A; /* Jacobian matrix */
24: Mat Jacp; /* JacobianP matrix */
25: Mat H; /* Hessian matrix for optimization */
26: Vec U, Lambda[1], Mup[1]; /* adjoint variables */
27: Vec Lambda2[1], Mup2[1]; /* second-order adjoint variables */
28: Vec Ihp1[1]; /* working space for Hessian evaluations */
29: Vec Ihp2[1]; /* working space for Hessian evaluations */
30: Vec Ihp3[1]; /* working space for Hessian evaluations */
31: Vec Ihp4[1]; /* working space for Hessian evaluations */
32: Vec Dir; /* direction vector */
33: PetscReal ob[2]; /* observation used by the cost function */
34: PetscBool implicitform; /* implicit ODE? */
35: };
37: PetscErrorCode FormFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
38: PetscErrorCode FormHessian(Tao, Vec, Mat, Mat, void *);
39: PetscErrorCode Adjoint2(Vec, PetscScalar[], User);
41: /* ----------------------- Explicit form of the ODE -------------------- */
43: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *ctx)
44: {
45: User user = (User)ctx;
46: PetscScalar *f;
47: const PetscScalar *u;
49: PetscFunctionBeginUser;
50: PetscCall(VecGetArrayRead(U, &u));
51: PetscCall(VecGetArray(F, &f));
52: f[0] = u[1];
53: f[1] = user->mu * ((1. - u[0] * u[0]) * u[1] - u[0]);
54: PetscCall(VecRestoreArrayRead(U, &u));
55: PetscCall(VecRestoreArray(F, &f));
56: PetscFunctionReturn(PETSC_SUCCESS);
57: }
59: static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, void *ctx)
60: {
61: User user = (User)ctx;
62: PetscReal mu = user->mu;
63: PetscInt rowcol[] = {0, 1};
64: PetscScalar J[2][2];
65: const PetscScalar *u;
67: PetscFunctionBeginUser;
68: PetscCall(VecGetArrayRead(U, &u));
70: J[0][0] = 0;
71: J[1][0] = -mu * (2.0 * u[1] * u[0] + 1.);
72: J[0][1] = 1.0;
73: J[1][1] = mu * (1.0 - u[0] * u[0]);
74: PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
75: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
76: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
77: if (B && A != B) {
78: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
79: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
80: }
82: PetscCall(VecRestoreArrayRead(U, &u));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
87: {
88: const PetscScalar *vl, *vr, *u;
89: PetscScalar *vhv;
90: PetscScalar dJdU[2][2][2] = {{{0}}};
91: PetscInt i, j, k;
92: User user = (User)ctx;
94: PetscFunctionBeginUser;
95: PetscCall(VecGetArrayRead(U, &u));
96: PetscCall(VecGetArrayRead(Vl[0], &vl));
97: PetscCall(VecGetArrayRead(Vr, &vr));
98: PetscCall(VecGetArray(VHV[0], &vhv));
100: dJdU[1][0][0] = -2. * user->mu * u[1];
101: dJdU[1][1][0] = -2. * user->mu * u[0];
102: dJdU[1][0][1] = -2. * user->mu * u[0];
103: for (j = 0; j < 2; j++) {
104: vhv[j] = 0;
105: for (k = 0; k < 2; k++)
106: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
107: }
109: PetscCall(VecRestoreArrayRead(U, &u));
110: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
111: PetscCall(VecRestoreArrayRead(Vr, &vr));
112: PetscCall(VecRestoreArray(VHV[0], &vhv));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode RHSHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
117: {
118: const PetscScalar *vl, *vr, *u;
119: PetscScalar *vhv;
120: PetscScalar dJdP[2][2][1] = {{{0}}};
121: PetscInt i, j, k;
123: PetscFunctionBeginUser;
124: PetscCall(VecGetArrayRead(U, &u));
125: PetscCall(VecGetArrayRead(Vl[0], &vl));
126: PetscCall(VecGetArrayRead(Vr, &vr));
127: PetscCall(VecGetArray(VHV[0], &vhv));
129: dJdP[1][0][0] = -(1. + 2. * u[0] * u[1]);
130: dJdP[1][1][0] = 1. - u[0] * u[0];
131: for (j = 0; j < 2; j++) {
132: vhv[j] = 0;
133: for (k = 0; k < 1; k++)
134: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
135: }
137: PetscCall(VecRestoreArrayRead(U, &u));
138: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
139: PetscCall(VecRestoreArrayRead(Vr, &vr));
140: PetscCall(VecRestoreArray(VHV[0], &vhv));
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode RHSHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
145: {
146: const PetscScalar *vl, *vr, *u;
147: PetscScalar *vhv;
148: PetscScalar dJdU[2][1][2] = {{{0}}};
149: PetscInt i, j, k;
151: PetscFunctionBeginUser;
152: PetscCall(VecGetArrayRead(U, &u));
153: PetscCall(VecGetArrayRead(Vl[0], &vl));
154: PetscCall(VecGetArrayRead(Vr, &vr));
155: PetscCall(VecGetArray(VHV[0], &vhv));
157: dJdU[1][0][0] = -1. - 2. * u[1] * u[0];
158: dJdU[1][0][1] = 1. - u[0] * u[0];
159: for (j = 0; j < 1; j++) {
160: vhv[j] = 0;
161: for (k = 0; k < 2; k++)
162: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
163: }
165: PetscCall(VecRestoreArrayRead(U, &u));
166: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
167: PetscCall(VecRestoreArrayRead(Vr, &vr));
168: PetscCall(VecRestoreArray(VHV[0], &vhv));
169: PetscFunctionReturn(PETSC_SUCCESS);
170: }
172: static PetscErrorCode RHSHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
173: {
174: PetscFunctionBeginUser;
175: PetscFunctionReturn(PETSC_SUCCESS);
176: }
178: /* ----------------------- Implicit form of the ODE -------------------- */
180: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
181: {
182: User user = (User)ctx;
183: PetscScalar *f;
184: const PetscScalar *u, *udot;
186: PetscFunctionBeginUser;
187: PetscCall(VecGetArrayRead(U, &u));
188: PetscCall(VecGetArrayRead(Udot, &udot));
189: PetscCall(VecGetArray(F, &f));
191: f[0] = udot[0] - u[1];
192: f[1] = udot[1] - user->mu * ((1.0 - u[0] * u[0]) * u[1] - u[0]);
194: PetscCall(VecRestoreArrayRead(U, &u));
195: PetscCall(VecRestoreArrayRead(Udot, &udot));
196: PetscCall(VecRestoreArray(F, &f));
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
201: {
202: User user = (User)ctx;
203: PetscInt rowcol[] = {0, 1};
204: PetscScalar J[2][2];
205: const PetscScalar *u;
207: PetscFunctionBeginUser;
208: PetscCall(VecGetArrayRead(U, &u));
210: J[0][0] = a;
211: J[0][1] = -1.0;
212: J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
213: J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
214: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
215: PetscCall(VecRestoreArrayRead(U, &u));
216: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
217: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
218: if (A != B) {
219: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
220: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
221: }
223: PetscCall(VecRestoreArrayRead(U, &u));
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec U, Mat A, void *ctx)
228: {
229: PetscInt row[] = {0, 1}, col[] = {0};
230: PetscScalar J[2][1];
231: const PetscScalar *u;
233: PetscFunctionBeginUser;
234: PetscCall(VecGetArrayRead(U, &u));
236: J[0][0] = 0;
237: J[1][0] = (1. - u[0] * u[0]) * u[1] - u[0];
238: PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
239: PetscCall(VecRestoreArrayRead(U, &u));
240: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
241: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
243: PetscCall(VecRestoreArrayRead(U, &u));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: static PetscErrorCode IHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
248: {
249: const PetscScalar *vl, *vr, *u;
250: PetscScalar *vhv;
251: PetscScalar dJdU[2][2][2] = {{{0}}};
252: PetscInt i, j, k;
253: User user = (User)ctx;
255: PetscFunctionBeginUser;
256: PetscCall(VecGetArrayRead(U, &u));
257: PetscCall(VecGetArrayRead(Vl[0], &vl));
258: PetscCall(VecGetArrayRead(Vr, &vr));
259: PetscCall(VecGetArray(VHV[0], &vhv));
261: dJdU[1][0][0] = 2. * user->mu * u[1];
262: dJdU[1][1][0] = 2. * user->mu * u[0];
263: dJdU[1][0][1] = 2. * user->mu * u[0];
264: for (j = 0; j < 2; j++) {
265: vhv[j] = 0;
266: for (k = 0; k < 2; k++)
267: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
268: }
270: PetscCall(VecRestoreArrayRead(U, &u));
271: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
272: PetscCall(VecRestoreArrayRead(Vr, &vr));
273: PetscCall(VecRestoreArray(VHV[0], &vhv));
274: PetscFunctionReturn(PETSC_SUCCESS);
275: }
277: static PetscErrorCode IHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
278: {
279: const PetscScalar *vl, *vr, *u;
280: PetscScalar *vhv;
281: PetscScalar dJdP[2][2][1] = {{{0}}};
282: PetscInt i, j, k;
284: PetscFunctionBeginUser;
285: PetscCall(VecGetArrayRead(U, &u));
286: PetscCall(VecGetArrayRead(Vl[0], &vl));
287: PetscCall(VecGetArrayRead(Vr, &vr));
288: PetscCall(VecGetArray(VHV[0], &vhv));
290: dJdP[1][0][0] = 1. + 2. * u[0] * u[1];
291: dJdP[1][1][0] = u[0] * u[0] - 1.;
292: for (j = 0; j < 2; j++) {
293: vhv[j] = 0;
294: for (k = 0; k < 1; k++)
295: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdP[i][j][k] * vr[k];
296: }
298: PetscCall(VecRestoreArrayRead(U, &u));
299: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
300: PetscCall(VecRestoreArrayRead(Vr, &vr));
301: PetscCall(VecRestoreArray(VHV[0], &vhv));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode IHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
306: {
307: const PetscScalar *vl, *vr, *u;
308: PetscScalar *vhv;
309: PetscScalar dJdU[2][1][2] = {{{0}}};
310: PetscInt i, j, k;
312: PetscFunctionBeginUser;
313: PetscCall(VecGetArrayRead(U, &u));
314: PetscCall(VecGetArrayRead(Vl[0], &vl));
315: PetscCall(VecGetArrayRead(Vr, &vr));
316: PetscCall(VecGetArray(VHV[0], &vhv));
318: dJdU[1][0][0] = 1. + 2. * u[1] * u[0];
319: dJdU[1][0][1] = u[0] * u[0] - 1.;
320: for (j = 0; j < 1; j++) {
321: vhv[j] = 0;
322: for (k = 0; k < 2; k++)
323: for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
324: }
326: PetscCall(VecRestoreArrayRead(U, &u));
327: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
328: PetscCall(VecRestoreArrayRead(Vr, &vr));
329: PetscCall(VecRestoreArray(VHV[0], &vhv));
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode IHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, void *ctx)
334: {
335: PetscFunctionBeginUser;
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
340: static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx)
341: {
342: const PetscScalar *x;
343: PetscReal tfinal, dt;
344: User user = (User)ctx;
345: Vec interpolatedX;
347: PetscFunctionBeginUser;
348: PetscCall(TSGetTimeStep(ts, &dt));
349: PetscCall(TSGetMaxTime(ts, &tfinal));
351: while (user->next_output <= t && user->next_output <= tfinal) {
352: PetscCall(VecDuplicate(X, &interpolatedX));
353: PetscCall(TSInterpolate(ts, user->next_output, interpolatedX));
354: PetscCall(VecGetArrayRead(interpolatedX, &x));
355: 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])));
356: PetscCall(VecRestoreArrayRead(interpolatedX, &x));
357: PetscCall(VecDestroy(&interpolatedX));
358: user->next_output += PetscRealConstant(0.1);
359: }
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: int main(int argc, char **argv)
364: {
365: Vec P;
366: PetscBool monitor = PETSC_FALSE;
367: PetscScalar *x_ptr;
368: const PetscScalar *y_ptr;
369: PetscMPIInt size;
370: struct _n_User user;
371: Tao tao;
372: KSP ksp;
373: PC pc;
375: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
376: Initialize program
377: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
378: PetscFunctionBeginUser;
379: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
380: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
381: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
383: /* Create TAO solver and set desired solution method */
384: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
385: PetscCall(TaoSetType(tao, TAOBQNLS));
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Set runtime options
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: user.next_output = 0.0;
391: user.mu = PetscRealConstant(1.0e3);
392: user.ftime = PetscRealConstant(0.5);
393: user.implicitform = PETSC_TRUE;
395: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
396: PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
397: PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &user.implicitform, NULL));
399: /* Create necessary matrix and vectors, solve same ODE on every process */
400: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.A));
401: PetscCall(MatSetSizes(user.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
402: PetscCall(MatSetFromOptions(user.A));
403: PetscCall(MatSetUp(user.A));
404: PetscCall(MatCreateVecs(user.A, &user.U, NULL));
405: PetscCall(MatCreateVecs(user.A, &user.Lambda[0], NULL));
406: PetscCall(MatCreateVecs(user.A, &user.Lambda2[0], NULL));
407: PetscCall(MatCreateVecs(user.A, &user.Ihp1[0], NULL));
408: PetscCall(MatCreateVecs(user.A, &user.Ihp2[0], NULL));
410: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Jacp));
411: PetscCall(MatSetSizes(user.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
412: PetscCall(MatSetFromOptions(user.Jacp));
413: PetscCall(MatSetUp(user.Jacp));
414: PetscCall(MatCreateVecs(user.Jacp, &user.Dir, NULL));
415: PetscCall(MatCreateVecs(user.Jacp, &user.Mup[0], NULL));
416: PetscCall(MatCreateVecs(user.Jacp, &user.Mup2[0], NULL));
417: PetscCall(MatCreateVecs(user.Jacp, &user.Ihp3[0], NULL));
418: PetscCall(MatCreateVecs(user.Jacp, &user.Ihp4[0], NULL));
420: /* Create timestepping solver context */
421: PetscCall(TSCreate(PETSC_COMM_WORLD, &user.ts));
422: PetscCall(TSSetEquationType(user.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
423: if (user.implicitform) {
424: PetscCall(TSSetIFunction(user.ts, NULL, IFunction, &user));
425: PetscCall(TSSetIJacobian(user.ts, user.A, user.A, IJacobian, &user));
426: PetscCall(TSSetType(user.ts, TSCN));
427: } else {
428: PetscCall(TSSetRHSFunction(user.ts, NULL, RHSFunction, &user));
429: PetscCall(TSSetRHSJacobian(user.ts, user.A, user.A, RHSJacobian, &user));
430: PetscCall(TSSetType(user.ts, TSRK));
431: }
432: PetscCall(TSSetRHSJacobianP(user.ts, user.Jacp, RHSJacobianP, &user));
433: PetscCall(TSSetMaxTime(user.ts, user.ftime));
434: PetscCall(TSSetExactFinalTime(user.ts, TS_EXACTFINALTIME_MATCHSTEP));
436: if (monitor) PetscCall(TSMonitorSet(user.ts, Monitor, &user, NULL));
438: /* Set ODE initial conditions */
439: PetscCall(VecGetArray(user.U, &x_ptr));
440: x_ptr[0] = 2.0;
441: x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
442: PetscCall(VecRestoreArray(user.U, &x_ptr));
443: PetscCall(TSSetTimeStep(user.ts, PetscRealConstant(0.001)));
445: /* Set runtime options */
446: PetscCall(TSSetFromOptions(user.ts));
448: PetscCall(TSSolve(user.ts, user.U));
449: PetscCall(VecGetArrayRead(user.U, &y_ptr));
450: user.ob[0] = y_ptr[0];
451: user.ob[1] = y_ptr[1];
452: PetscCall(VecRestoreArrayRead(user.U, &y_ptr));
454: /* Save trajectory of solution so that TSAdjointSolve() may be used.
455: Skip checkpointing for the first TSSolve since no adjoint run follows it.
456: */
457: PetscCall(TSSetSaveTrajectory(user.ts));
459: /* Optimization starts */
460: PetscCall(MatCreate(PETSC_COMM_WORLD, &user.H));
461: PetscCall(MatSetSizes(user.H, PETSC_DECIDE, PETSC_DECIDE, 1, 1));
462: PetscCall(MatSetUp(user.H)); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */
464: /* Set initial solution guess */
465: PetscCall(MatCreateVecs(user.Jacp, &P, NULL));
466: PetscCall(VecGetArray(P, &x_ptr));
467: x_ptr[0] = PetscRealConstant(1.2);
468: PetscCall(VecRestoreArray(P, &x_ptr));
469: PetscCall(TaoSetSolution(tao, P));
471: /* Set routine for function and gradient evaluation */
472: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
473: PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
475: /* Check for any TAO command line options */
476: PetscCall(TaoGetKSP(tao, &ksp));
477: if (ksp) {
478: PetscCall(KSPGetPC(ksp, &pc));
479: PetscCall(PCSetType(pc, PCNONE));
480: }
481: PetscCall(TaoSetFromOptions(tao));
483: PetscCall(TaoSolve(tao));
485: PetscCall(VecView(P, PETSC_VIEWER_STDOUT_WORLD));
486: /* Free TAO data structures */
487: PetscCall(TaoDestroy(&tao));
489: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
490: Free work space. All PETSc objects should be destroyed when they
491: are no longer needed.
492: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
493: PetscCall(MatDestroy(&user.H));
494: PetscCall(MatDestroy(&user.A));
495: PetscCall(VecDestroy(&user.U));
496: PetscCall(MatDestroy(&user.Jacp));
497: PetscCall(VecDestroy(&user.Lambda[0]));
498: PetscCall(VecDestroy(&user.Mup[0]));
499: PetscCall(VecDestroy(&user.Lambda2[0]));
500: PetscCall(VecDestroy(&user.Mup2[0]));
501: PetscCall(VecDestroy(&user.Ihp1[0]));
502: PetscCall(VecDestroy(&user.Ihp2[0]));
503: PetscCall(VecDestroy(&user.Ihp3[0]));
504: PetscCall(VecDestroy(&user.Ihp4[0]));
505: PetscCall(VecDestroy(&user.Dir));
506: PetscCall(TSDestroy(&user.ts));
507: PetscCall(VecDestroy(&P));
508: PetscCall(PetscFinalize());
509: return 0;
510: }
512: /* ------------------------------------------------------------------ */
513: /*
514: FormFunctionGradient - Evaluates the function and corresponding gradient.
516: Input Parameters:
517: tao - the Tao context
518: X - the input vector
519: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient()
521: Output Parameters:
522: f - the newly evaluated function
523: G - the newly evaluated gradient
524: */
525: PetscErrorCode FormFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx)
526: {
527: User user_ptr = (User)ctx;
528: TS ts = user_ptr->ts;
529: PetscScalar *x_ptr, *g;
530: const PetscScalar *y_ptr;
532: PetscFunctionBeginUser;
533: PetscCall(VecGetArrayRead(P, &y_ptr));
534: user_ptr->mu = y_ptr[0];
535: PetscCall(VecRestoreArrayRead(P, &y_ptr));
537: PetscCall(TSSetTime(ts, 0.0));
538: PetscCall(TSSetStepNumber(ts, 0));
539: PetscCall(TSSetTimeStep(ts, PetscRealConstant(0.001))); /* can be overwritten by command line options */
540: PetscCall(TSSetFromOptions(ts));
541: PetscCall(VecGetArray(user_ptr->U, &x_ptr));
542: x_ptr[0] = 2.0;
543: x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user_ptr->mu) - 292.0 / (2187.0 * user_ptr->mu * user_ptr->mu);
544: PetscCall(VecRestoreArray(user_ptr->U, &x_ptr));
546: PetscCall(TSSolve(ts, user_ptr->U));
548: PetscCall(VecGetArrayRead(user_ptr->U, &y_ptr));
549: *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]);
551: /* Reset initial conditions for the adjoint integration */
552: PetscCall(VecGetArray(user_ptr->Lambda[0], &x_ptr));
553: x_ptr[0] = 2. * (y_ptr[0] - user_ptr->ob[0]);
554: x_ptr[1] = 2. * (y_ptr[1] - user_ptr->ob[1]);
555: PetscCall(VecRestoreArrayRead(user_ptr->U, &y_ptr));
556: PetscCall(VecRestoreArray(user_ptr->Lambda[0], &x_ptr));
558: PetscCall(VecGetArray(user_ptr->Mup[0], &x_ptr));
559: x_ptr[0] = 0.0;
560: PetscCall(VecRestoreArray(user_ptr->Mup[0], &x_ptr));
561: PetscCall(TSSetCostGradients(ts, 1, user_ptr->Lambda, user_ptr->Mup));
563: PetscCall(TSAdjointSolve(ts));
565: PetscCall(VecGetArray(user_ptr->Mup[0], &x_ptr));
566: PetscCall(VecGetArrayRead(user_ptr->Lambda[0], &y_ptr));
567: PetscCall(VecGetArray(G, &g));
568: 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];
569: PetscCall(VecRestoreArray(user_ptr->Mup[0], &x_ptr));
570: PetscCall(VecRestoreArrayRead(user_ptr->Lambda[0], &y_ptr));
571: PetscCall(VecRestoreArray(G, &g));
572: PetscFunctionReturn(PETSC_SUCCESS);
573: }
575: PetscErrorCode FormHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
576: {
577: User user_ptr = (User)ctx;
578: PetscScalar harr[1];
579: const PetscInt rows[1] = {0};
580: PetscInt col = 0;
582: PetscFunctionBeginUser;
583: PetscCall(Adjoint2(P, harr, user_ptr));
584: PetscCall(MatSetValues(H, 1, rows, 1, &col, harr, INSERT_VALUES));
586: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
587: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
588: if (H != Hpre) {
589: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
590: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
591: }
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: PetscErrorCode Adjoint2(Vec P, PetscScalar arr[], User ctx)
596: {
597: TS ts = ctx->ts;
598: const PetscScalar *z_ptr;
599: PetscScalar *x_ptr, *y_ptr, dzdp, dzdp2;
600: Mat tlmsen;
602: PetscFunctionBeginUser;
603: /* Reset TSAdjoint so that AdjointSetUp will be called again */
604: PetscCall(TSAdjointReset(ts));
606: /* The directional vector should be 1 since it is one-dimensional */
607: PetscCall(VecGetArray(ctx->Dir, &x_ptr));
608: x_ptr[0] = 1.;
609: PetscCall(VecRestoreArray(ctx->Dir, &x_ptr));
611: PetscCall(VecGetArrayRead(P, &z_ptr));
612: ctx->mu = z_ptr[0];
613: PetscCall(VecRestoreArrayRead(P, &z_ptr));
615: dzdp = -10.0 / (81.0 * ctx->mu * ctx->mu) + 2.0 * 292.0 / (2187.0 * ctx->mu * ctx->mu * ctx->mu);
616: 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);
618: PetscCall(TSSetTime(ts, 0.0));
619: PetscCall(TSSetStepNumber(ts, 0));
620: PetscCall(TSSetTimeStep(ts, PetscRealConstant(0.001))); /* can be overwritten by command line options */
621: PetscCall(TSSetFromOptions(ts));
622: PetscCall(TSSetCostHessianProducts(ts, 1, ctx->Lambda2, ctx->Mup2, ctx->Dir));
624: PetscCall(MatZeroEntries(ctx->Jacp));
625: PetscCall(MatSetValue(ctx->Jacp, 1, 0, dzdp, INSERT_VALUES));
626: PetscCall(MatAssemblyBegin(ctx->Jacp, MAT_FINAL_ASSEMBLY));
627: PetscCall(MatAssemblyEnd(ctx->Jacp, MAT_FINAL_ASSEMBLY));
629: PetscCall(TSAdjointSetForward(ts, ctx->Jacp));
630: PetscCall(VecGetArray(ctx->U, &y_ptr));
631: y_ptr[0] = 2.0;
632: y_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * ctx->mu) - 292.0 / (2187.0 * ctx->mu * ctx->mu);
633: PetscCall(VecRestoreArray(ctx->U, &y_ptr));
634: PetscCall(TSSolve(ts, ctx->U));
636: /* Set terminal conditions for first- and second-order adjonts */
637: PetscCall(VecGetArrayRead(ctx->U, &z_ptr));
638: PetscCall(VecGetArray(ctx->Lambda[0], &y_ptr));
639: y_ptr[0] = 2. * (z_ptr[0] - ctx->ob[0]);
640: y_ptr[1] = 2. * (z_ptr[1] - ctx->ob[1]);
641: PetscCall(VecRestoreArray(ctx->Lambda[0], &y_ptr));
642: PetscCall(VecRestoreArrayRead(ctx->U, &z_ptr));
643: PetscCall(VecGetArray(ctx->Mup[0], &y_ptr));
644: y_ptr[0] = 0.0;
645: PetscCall(VecRestoreArray(ctx->Mup[0], &y_ptr));
646: PetscCall(TSForwardGetSensitivities(ts, NULL, &tlmsen));
647: PetscCall(MatDenseGetColumn(tlmsen, 0, &x_ptr));
648: PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
649: y_ptr[0] = 2. * x_ptr[0];
650: y_ptr[1] = 2. * x_ptr[1];
651: PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
652: PetscCall(VecGetArray(ctx->Mup2[0], &y_ptr));
653: y_ptr[0] = 0.0;
654: PetscCall(VecRestoreArray(ctx->Mup2[0], &y_ptr));
655: PetscCall(MatDenseRestoreColumn(tlmsen, &x_ptr));
656: PetscCall(TSSetCostGradients(ts, 1, ctx->Lambda, ctx->Mup));
657: if (ctx->implicitform) {
658: PetscCall(TSSetIHessianProduct(ts, ctx->Ihp1, IHessianProductUU, ctx->Ihp2, IHessianProductUP, ctx->Ihp3, IHessianProductPU, ctx->Ihp4, IHessianProductPP, ctx));
659: } else {
660: PetscCall(TSSetRHSHessianProduct(ts, ctx->Ihp1, RHSHessianProductUU, ctx->Ihp2, RHSHessianProductUP, ctx->Ihp3, RHSHessianProductPU, ctx->Ihp4, RHSHessianProductPP, ctx));
661: }
662: PetscCall(TSAdjointSolve(ts));
664: PetscCall(VecGetArray(ctx->Lambda[0], &x_ptr));
665: PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
666: PetscCall(VecGetArrayRead(ctx->Mup2[0], &z_ptr));
668: arr[0] = x_ptr[1] * dzdp2 + y_ptr[1] * dzdp2 + z_ptr[0];
670: PetscCall(VecRestoreArray(ctx->Lambda2[0], &x_ptr));
671: PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
672: PetscCall(VecRestoreArrayRead(ctx->Mup2[0], &z_ptr));
674: /* Disable second-order adjoint mode */
675: PetscCall(TSAdjointReset(ts));
676: PetscCall(TSAdjointResetForward(ts));
677: PetscFunctionReturn(PETSC_SUCCESS);
678: }
680: /*TEST
681: build:
682: requires: !complex !single
683: test:
684: args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view
685: output_file: output/ex20opt_p_1.out
687: test:
688: suffix: 2
689: 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
690: output_file: output/ex20opt_p_2.out
692: test:
693: suffix: 3
694: args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view
695: output_file: output/ex20opt_p_3.out
697: test:
698: suffix: 4
699: 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
700: output_file: output/ex20opt_p_4.out
702: TEST*/