1 static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n";
2
3 /*
4 Notes:
5 This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS.
6 The nonlinear problem is written in an ODE equivalent form.
7 The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions.
8 The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details.
9 */
10
11 #include <petsctao.h>
12 #include <petscts.h>
13
14 typedef struct _n_User *User;
15 struct _n_User {
16 TS ts;
17 PetscReal mu;
18 PetscReal next_output;
19
20 /* Sensitivity analysis support */
21 PetscInt steps;
22 PetscReal ftime;
23 Mat A; /* Jacobian matrix for ODE */
24 Mat Jacp; /* JacobianP matrix for ODE*/
25 Mat H; /* Hessian matrix for optimization */
26 Vec U, Lambda[1], Mup[1]; /* first-order adjoint variables */
27 Vec Lambda2[2]; /* second-order adjoint variables */
28 Vec Ihp1[1]; /* working space for Hessian evaluations */
29 Vec Dir; /* direction vector */
30 PetscReal ob[2]; /* observation used by the cost function */
31 PetscBool implicitform; /* implicit ODE? */
32 };
33 PetscErrorCode Adjoint2(Vec, PetscScalar[], User);
34
35 /* ----------------------- Explicit form of the ODE -------------------- */
36
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,PetscCtx ctx)37 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
38 {
39 User user = (User)ctx;
40 PetscScalar *f;
41 const PetscScalar *u;
42
43 PetscFunctionBeginUser;
44 PetscCall(VecGetArrayRead(U, &u));
45 PetscCall(VecGetArray(F, &f));
46 f[0] = u[1];
47 f[1] = user->mu * ((1. - u[0] * u[0]) * u[1] - u[0]);
48 PetscCall(VecRestoreArrayRead(U, &u));
49 PetscCall(VecRestoreArray(F, &f));
50 PetscFunctionReturn(PETSC_SUCCESS);
51 }
52
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,PetscCtx ctx)53 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, PetscCtx ctx)
54 {
55 User user = (User)ctx;
56 PetscReal mu = user->mu;
57 PetscInt rowcol[] = {0, 1};
58 PetscScalar J[2][2];
59 const PetscScalar *u;
60
61 PetscFunctionBeginUser;
62 PetscCall(VecGetArrayRead(U, &u));
63 J[0][0] = 0;
64 J[1][0] = -mu * (2.0 * u[1] * u[0] + 1.);
65 J[0][1] = 1.0;
66 J[1][1] = mu * (1.0 - u[0] * u[0]);
67 PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
68 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
69 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
70 if (A != B) {
71 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
72 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
73 }
74 PetscCall(VecRestoreArrayRead(U, &u));
75 PetscFunctionReturn(PETSC_SUCCESS);
76 }
77
RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV,PetscCtx ctx)78 static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
79 {
80 const PetscScalar *vl, *vr, *u;
81 PetscScalar *vhv;
82 PetscScalar dJdU[2][2][2] = {{{0}}};
83 PetscInt i, j, k;
84 User user = (User)ctx;
85
86 PetscFunctionBeginUser;
87 PetscCall(VecGetArrayRead(U, &u));
88 PetscCall(VecGetArrayRead(Vl[0], &vl));
89 PetscCall(VecGetArrayRead(Vr, &vr));
90 PetscCall(VecGetArray(VHV[0], &vhv));
91
92 dJdU[1][0][0] = -2. * user->mu * u[1];
93 dJdU[1][1][0] = -2. * user->mu * u[0];
94 dJdU[1][0][1] = -2. * user->mu * u[0];
95 for (j = 0; j < 2; j++) {
96 vhv[j] = 0;
97 for (k = 0; k < 2; k++)
98 for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
99 }
100 PetscCall(VecRestoreArrayRead(U, &u));
101 PetscCall(VecRestoreArrayRead(Vl[0], &vl));
102 PetscCall(VecRestoreArrayRead(Vr, &vr));
103 PetscCall(VecRestoreArray(VHV[0], &vhv));
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
107 /* ----------------------- Implicit form of the ODE -------------------- */
108
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,PetscCtx ctx)109 static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, PetscCtx ctx)
110 {
111 User user = (User)ctx;
112 const PetscScalar *u, *udot;
113 PetscScalar *f;
114
115 PetscFunctionBeginUser;
116 PetscCall(VecGetArrayRead(U, &u));
117 PetscCall(VecGetArrayRead(Udot, &udot));
118 PetscCall(VecGetArray(F, &f));
119 f[0] = udot[0] - u[1];
120 f[1] = udot[1] - user->mu * ((1.0 - u[0] * u[0]) * u[1] - u[0]);
121 PetscCall(VecRestoreArrayRead(U, &u));
122 PetscCall(VecRestoreArrayRead(Udot, &udot));
123 PetscCall(VecRestoreArray(F, &f));
124 PetscFunctionReturn(PETSC_SUCCESS);
125 }
126
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,PetscCtx ctx)127 static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, PetscCtx ctx)
128 {
129 User user = (User)ctx;
130 PetscInt rowcol[] = {0, 1};
131 PetscScalar J[2][2];
132 const PetscScalar *u;
133
134 PetscFunctionBeginUser;
135 PetscCall(VecGetArrayRead(U, &u));
136 J[0][0] = a;
137 J[0][1] = -1.0;
138 J[1][0] = user->mu * (1.0 + 2.0 * u[0] * u[1]);
139 J[1][1] = a - user->mu * (1.0 - u[0] * u[0]);
140 PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
141 PetscCall(VecRestoreArrayRead(U, &u));
142 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
143 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
144 if (A != B) {
145 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
146 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
147 }
148 PetscFunctionReturn(PETSC_SUCCESS);
149 }
150
151 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
Monitor(TS ts,PetscInt step,PetscReal t,Vec U,PetscCtx ctx)152 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec U, PetscCtx ctx)
153 {
154 const PetscScalar *u;
155 PetscReal tfinal, dt;
156 User user = (User)ctx;
157 Vec interpolatedU;
158
159 PetscFunctionBeginUser;
160 PetscCall(TSGetTimeStep(ts, &dt));
161 PetscCall(TSGetMaxTime(ts, &tfinal));
162
163 while (user->next_output <= t && user->next_output <= tfinal) {
164 PetscCall(VecDuplicate(U, &interpolatedU));
165 PetscCall(TSInterpolate(ts, user->next_output, interpolatedU));
166 PetscCall(VecGetArrayRead(interpolatedU, &u));
167 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(u[0]), (double)PetscRealPart(u[1])));
168 PetscCall(VecRestoreArrayRead(interpolatedU, &u));
169 PetscCall(VecDestroy(&interpolatedU));
170 user->next_output += 0.1;
171 }
172 PetscFunctionReturn(PETSC_SUCCESS);
173 }
174
IHessianProductUU(TS ts,PetscReal t,Vec U,Vec * Vl,Vec Vr,Vec * VHV,PetscCtx ctx)175 static PetscErrorCode IHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
176 {
177 const PetscScalar *vl, *vr, *u;
178 PetscScalar *vhv;
179 PetscScalar dJdU[2][2][2] = {{{0}}};
180 PetscInt i, j, k;
181 User user = (User)ctx;
182
183 PetscFunctionBeginUser;
184 PetscCall(VecGetArrayRead(U, &u));
185 PetscCall(VecGetArrayRead(Vl[0], &vl));
186 PetscCall(VecGetArrayRead(Vr, &vr));
187 PetscCall(VecGetArray(VHV[0], &vhv));
188 dJdU[1][0][0] = 2. * user->mu * u[1];
189 dJdU[1][1][0] = 2. * user->mu * u[0];
190 dJdU[1][0][1] = 2. * user->mu * u[0];
191 for (j = 0; j < 2; j++) {
192 vhv[j] = 0;
193 for (k = 0; k < 2; k++)
194 for (i = 0; i < 2; i++) vhv[j] += vl[i] * dJdU[i][j][k] * vr[k];
195 }
196 PetscCall(VecRestoreArrayRead(U, &u));
197 PetscCall(VecRestoreArrayRead(Vl[0], &vl));
198 PetscCall(VecRestoreArrayRead(Vr, &vr));
199 PetscCall(VecRestoreArray(VHV[0], &vhv));
200 PetscFunctionReturn(PETSC_SUCCESS);
201 }
202
203 /* ------------------ User-defined routines for TAO -------------------------- */
204
FormFunctionGradient(Tao tao,Vec IC,PetscReal * f,Vec G,PetscCtx ctx)205 static PetscErrorCode FormFunctionGradient(Tao tao, Vec IC, PetscReal *f, Vec G, PetscCtx ctx)
206 {
207 User user_ptr = (User)ctx;
208 TS ts = user_ptr->ts;
209 const PetscScalar *x_ptr;
210 PetscScalar *y_ptr;
211
212 PetscFunctionBeginUser;
213 PetscCall(VecCopy(IC, user_ptr->U)); /* set up the initial condition */
214
215 PetscCall(TSSetTime(ts, 0.0));
216 PetscCall(TSSetStepNumber(ts, 0));
217 PetscCall(TSSetTimeStep(ts, 0.001)); /* can be overwritten by command line options */
218 PetscCall(TSSetFromOptions(ts));
219 PetscCall(TSSolve(ts, user_ptr->U));
220
221 PetscCall(VecGetArrayRead(user_ptr->U, &x_ptr));
222 PetscCall(VecGetArray(user_ptr->Lambda[0], &y_ptr));
223 *f = (x_ptr[0] - user_ptr->ob[0]) * (x_ptr[0] - user_ptr->ob[0]) + (x_ptr[1] - user_ptr->ob[1]) * (x_ptr[1] - user_ptr->ob[1]);
224 y_ptr[0] = 2. * (x_ptr[0] - user_ptr->ob[0]);
225 y_ptr[1] = 2. * (x_ptr[1] - user_ptr->ob[1]);
226 PetscCall(VecRestoreArray(user_ptr->Lambda[0], &y_ptr));
227 PetscCall(VecRestoreArrayRead(user_ptr->U, &x_ptr));
228
229 PetscCall(TSSetCostGradients(ts, 1, user_ptr->Lambda, NULL));
230 PetscCall(TSAdjointSolve(ts));
231 PetscCall(VecCopy(user_ptr->Lambda[0], G));
232 PetscFunctionReturn(PETSC_SUCCESS);
233 }
234
FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,PetscCtx ctx)235 static PetscErrorCode FormHessian(Tao tao, Vec U, Mat H, Mat Hpre, PetscCtx ctx)
236 {
237 User user_ptr = (User)ctx;
238 PetscScalar harr[2];
239 PetscScalar *x_ptr;
240 const PetscInt rows[2] = {0, 1};
241 PetscInt col;
242
243 PetscFunctionBeginUser;
244 PetscCall(VecCopy(U, user_ptr->U));
245 PetscCall(VecGetArray(user_ptr->Dir, &x_ptr));
246 x_ptr[0] = 1.;
247 x_ptr[1] = 0.;
248 PetscCall(VecRestoreArray(user_ptr->Dir, &x_ptr));
249 PetscCall(Adjoint2(user_ptr->U, harr, user_ptr));
250 col = 0;
251 PetscCall(MatSetValues(H, 2, rows, 1, &col, harr, INSERT_VALUES));
252
253 PetscCall(VecCopy(U, user_ptr->U));
254 PetscCall(VecGetArray(user_ptr->Dir, &x_ptr));
255 x_ptr[0] = 0.;
256 x_ptr[1] = 1.;
257 PetscCall(VecRestoreArray(user_ptr->Dir, &x_ptr));
258 PetscCall(Adjoint2(user_ptr->U, harr, user_ptr));
259 col = 1;
260 PetscCall(MatSetValues(H, 2, rows, 1, &col, harr, INSERT_VALUES));
261
262 PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
263 PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
264 if (H != Hpre) {
265 PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
266 PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
267 }
268 PetscFunctionReturn(PETSC_SUCCESS);
269 }
270
MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,PetscCtx ctx)271 static PetscErrorCode MatrixFreeHessian(Tao tao, Vec U, Mat H, Mat Hpre, PetscCtx ctx)
272 {
273 User user_ptr = (User)ctx;
274
275 PetscFunctionBeginUser;
276 PetscCall(VecCopy(U, user_ptr->U));
277 PetscFunctionReturn(PETSC_SUCCESS);
278 }
279
280 /* ------------ Routines calculating second-order derivatives -------------- */
281
282 /*
283 Compute the Hessian-vector product for the cost function using Second-order adjoint
284 */
Adjoint2(Vec U,PetscScalar arr[],User ctx)285 PetscErrorCode Adjoint2(Vec U, PetscScalar arr[], User ctx)
286 {
287 TS ts = ctx->ts;
288 PetscScalar *x_ptr, *y_ptr;
289 Mat tlmsen;
290
291 PetscFunctionBeginUser;
292 PetscCall(TSAdjointReset(ts));
293
294 PetscCall(TSSetTime(ts, 0.0));
295 PetscCall(TSSetStepNumber(ts, 0));
296 PetscCall(TSSetTimeStep(ts, 0.001));
297 PetscCall(TSSetFromOptions(ts));
298 PetscCall(TSSetCostHessianProducts(ts, 1, ctx->Lambda2, NULL, ctx->Dir));
299 PetscCall(TSAdjointSetForward(ts, NULL));
300 PetscCall(TSSolve(ts, U));
301
302 /* Set terminal conditions for first- and second-order adjonts */
303 PetscCall(VecGetArray(U, &x_ptr));
304 PetscCall(VecGetArray(ctx->Lambda[0], &y_ptr));
305 y_ptr[0] = 2. * (x_ptr[0] - ctx->ob[0]);
306 y_ptr[1] = 2. * (x_ptr[1] - ctx->ob[1]);
307 PetscCall(VecRestoreArray(ctx->Lambda[0], &y_ptr));
308 PetscCall(VecRestoreArray(U, &x_ptr));
309 PetscCall(TSForwardGetSensitivities(ts, NULL, &tlmsen));
310 PetscCall(MatDenseGetColumn(tlmsen, 0, &x_ptr));
311 PetscCall(VecGetArray(ctx->Lambda2[0], &y_ptr));
312 y_ptr[0] = 2. * x_ptr[0];
313 y_ptr[1] = 2. * x_ptr[1];
314 PetscCall(VecRestoreArray(ctx->Lambda2[0], &y_ptr));
315 PetscCall(MatDenseRestoreColumn(tlmsen, &x_ptr));
316
317 PetscCall(TSSetCostGradients(ts, 1, ctx->Lambda, NULL));
318 if (ctx->implicitform) {
319 PetscCall(TSSetIHessianProduct(ts, ctx->Ihp1, IHessianProductUU, NULL, NULL, NULL, NULL, NULL, NULL, ctx));
320 } else {
321 PetscCall(TSSetRHSHessianProduct(ts, ctx->Ihp1, RHSHessianProductUU, NULL, NULL, NULL, NULL, NULL, NULL, ctx));
322 }
323 PetscCall(TSAdjointSolve(ts));
324
325 PetscCall(VecGetArray(ctx->Lambda2[0], &x_ptr));
326 arr[0] = x_ptr[0];
327 arr[1] = x_ptr[1];
328 PetscCall(VecRestoreArray(ctx->Lambda2[0], &x_ptr));
329
330 PetscCall(TSAdjointReset(ts));
331 PetscCall(TSAdjointResetForward(ts));
332 PetscFunctionReturn(PETSC_SUCCESS);
333 }
334
FiniteDiff(Vec U,PetscScalar arr[],User ctx)335 PetscErrorCode FiniteDiff(Vec U, PetscScalar arr[], User ctx)
336 {
337 Vec Up, G, Gp;
338 const PetscScalar eps = PetscRealConstant(1e-7);
339 PetscScalar *u;
340 Tao tao = NULL;
341 PetscReal f;
342
343 PetscFunctionBeginUser;
344 PetscCall(VecDuplicate(U, &Up));
345 PetscCall(VecDuplicate(U, &G));
346 PetscCall(VecDuplicate(U, &Gp));
347
348 PetscCall(FormFunctionGradient(tao, U, &f, G, ctx));
349
350 PetscCall(VecCopy(U, Up));
351 PetscCall(VecGetArray(Up, &u));
352 u[0] += eps;
353 PetscCall(VecRestoreArray(Up, &u));
354 PetscCall(FormFunctionGradient(tao, Up, &f, Gp, ctx));
355 PetscCall(VecAXPY(Gp, -1, G));
356 PetscCall(VecScale(Gp, 1. / eps));
357 PetscCall(VecGetArray(Gp, &u));
358 arr[0] = u[0];
359 arr[1] = u[1];
360 PetscCall(VecRestoreArray(Gp, &u));
361
362 PetscCall(VecCopy(U, Up));
363 PetscCall(VecGetArray(Up, &u));
364 u[1] += eps;
365 PetscCall(VecRestoreArray(Up, &u));
366 PetscCall(FormFunctionGradient(tao, Up, &f, Gp, ctx));
367 PetscCall(VecAXPY(Gp, -1, G));
368 PetscCall(VecScale(Gp, 1. / eps));
369 PetscCall(VecGetArray(Gp, &u));
370 arr[2] = u[0];
371 arr[3] = u[1];
372 PetscCall(VecRestoreArray(Gp, &u));
373
374 PetscCall(VecDestroy(&G));
375 PetscCall(VecDestroy(&Gp));
376 PetscCall(VecDestroy(&Up));
377 PetscFunctionReturn(PETSC_SUCCESS);
378 }
379
HessianProductMat(Mat mat,Vec svec,Vec y)380 static PetscErrorCode HessianProductMat(Mat mat, Vec svec, Vec y)
381 {
382 User user_ptr;
383 PetscScalar *y_ptr;
384
385 PetscFunctionBeginUser;
386 PetscCall(MatShellGetContext(mat, &user_ptr));
387 PetscCall(VecCopy(svec, user_ptr->Dir));
388 PetscCall(VecGetArray(y, &y_ptr));
389 PetscCall(Adjoint2(user_ptr->U, y_ptr, user_ptr));
390 PetscCall(VecRestoreArray(y, &y_ptr));
391 PetscFunctionReturn(PETSC_SUCCESS);
392 }
393
main(int argc,char ** argv)394 int main(int argc, char **argv)
395 {
396 PetscBool monitor = PETSC_FALSE, mf = PETSC_TRUE;
397 PetscInt mode = 0;
398 PetscMPIInt size;
399 struct _n_User user;
400 Vec x; /* working vector for TAO */
401 PetscScalar *x_ptr, arr[4];
402 PetscScalar ic1 = 2.2, ic2 = -0.7; /* initial guess for TAO */
403 Tao tao;
404 KSP ksp;
405 PC pc;
406
407 /* Initialize program */
408 PetscFunctionBeginUser;
409 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
410 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
411 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
412
413 /* Set runtime options */
414 user.next_output = 0.0;
415 user.mu = 1.0e3;
416 user.steps = 0;
417 user.ftime = 0.5;
418 user.implicitform = PETSC_TRUE;
419
420 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL));
421 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL));
422 PetscCall(PetscOptionsGetInt(NULL, NULL, "-mode", &mode, NULL));
423 PetscCall(PetscOptionsGetReal(NULL, NULL, "-ic1", &ic1, NULL));
424 PetscCall(PetscOptionsGetReal(NULL, NULL, "-ic2", &ic2, NULL));
425 PetscCall(PetscOptionsGetBool(NULL, NULL, "-my_tao_mf", &mf, NULL)); /* matrix-free hessian for optimization */
426 PetscCall(PetscOptionsGetBool(NULL, NULL, "-implicitform", &user.implicitform, NULL));
427
428 /* Create necessary matrix and vectors, solve same ODE on every process */
429 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.A));
430 PetscCall(MatSetSizes(user.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
431 PetscCall(MatSetFromOptions(user.A));
432 PetscCall(MatSetUp(user.A));
433 PetscCall(MatCreateVecs(user.A, &user.U, NULL));
434 PetscCall(MatCreateVecs(user.A, &user.Dir, NULL));
435 PetscCall(MatCreateVecs(user.A, &user.Lambda[0], NULL));
436 PetscCall(MatCreateVecs(user.A, &user.Lambda2[0], NULL));
437 PetscCall(MatCreateVecs(user.A, &user.Ihp1[0], NULL));
438
439 /* Create timestepping solver context */
440 PetscCall(TSCreate(PETSC_COMM_WORLD, &user.ts));
441 PetscCall(TSSetEquationType(user.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
442 if (user.implicitform) {
443 PetscCall(TSSetIFunction(user.ts, NULL, IFunction, &user));
444 PetscCall(TSSetIJacobian(user.ts, user.A, user.A, IJacobian, &user));
445 PetscCall(TSSetType(user.ts, TSCN));
446 } else {
447 PetscCall(TSSetRHSFunction(user.ts, NULL, RHSFunction, &user));
448 PetscCall(TSSetRHSJacobian(user.ts, user.A, user.A, RHSJacobian, &user));
449 PetscCall(TSSetType(user.ts, TSRK));
450 }
451 PetscCall(TSSetMaxTime(user.ts, user.ftime));
452 PetscCall(TSSetExactFinalTime(user.ts, TS_EXACTFINALTIME_MATCHSTEP));
453
454 if (monitor) PetscCall(TSMonitorSet(user.ts, Monitor, &user, NULL));
455
456 /* Set ODE initial conditions */
457 PetscCall(VecGetArray(user.U, &x_ptr));
458 x_ptr[0] = 2.0;
459 x_ptr[1] = -2.0 / 3.0 + 10.0 / (81.0 * user.mu) - 292.0 / (2187.0 * user.mu * user.mu);
460 PetscCall(VecRestoreArray(user.U, &x_ptr));
461
462 /* Set runtime options */
463 PetscCall(TSSetFromOptions(user.ts));
464
465 /* Obtain the observation */
466 PetscCall(TSSolve(user.ts, user.U));
467 PetscCall(VecGetArray(user.U, &x_ptr));
468 user.ob[0] = x_ptr[0];
469 user.ob[1] = x_ptr[1];
470 PetscCall(VecRestoreArray(user.U, &x_ptr));
471
472 PetscCall(VecDuplicate(user.U, &x));
473 PetscCall(VecGetArray(x, &x_ptr));
474 x_ptr[0] = ic1;
475 x_ptr[1] = ic2;
476 PetscCall(VecRestoreArray(x, &x_ptr));
477
478 /* Save trajectory of solution so that TSAdjointSolve() may be used */
479 PetscCall(TSSetSaveTrajectory(user.ts));
480
481 /* Compare finite difference and second-order adjoint. */
482 switch (mode) {
483 case 2:
484 PetscCall(FiniteDiff(x, arr, &user));
485 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n Finite difference approximation of the Hessian\n"));
486 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g %g\n%g %g\n", (double)arr[0], (double)arr[1], (double)arr[2], (double)arr[3]));
487 break;
488 case 1: /* Compute the Hessian column by column */
489 PetscCall(VecCopy(x, user.U));
490 PetscCall(VecGetArray(user.Dir, &x_ptr));
491 x_ptr[0] = 1.;
492 x_ptr[1] = 0.;
493 PetscCall(VecRestoreArray(user.Dir, &x_ptr));
494 PetscCall(Adjoint2(user.U, arr, &user));
495 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nFirst column of the Hessian\n"));
496 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g\n%g\n", (double)arr[0], (double)arr[1]));
497 PetscCall(VecCopy(x, user.U));
498 PetscCall(VecGetArray(user.Dir, &x_ptr));
499 x_ptr[0] = 0.;
500 x_ptr[1] = 1.;
501 PetscCall(VecRestoreArray(user.Dir, &x_ptr));
502 PetscCall(Adjoint2(user.U, arr, &user));
503 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSecond column of the Hessian\n"));
504 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%g\n%g\n", (double)arr[0], (double)arr[1]));
505 break;
506 case 0:
507 /* Create optimization context and set up */
508 PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
509 PetscCall(TaoSetType(tao, TAOBLMVM));
510 PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormFunctionGradient, (void *)&user));
511
512 if (mf) {
513 PetscCall(MatCreateShell(PETSC_COMM_SELF, 2, 2, 2, 2, (void *)&user, &user.H));
514 PetscCall(MatShellSetOperation(user.H, MATOP_MULT, (PetscErrorCodeFn *)HessianProductMat));
515 PetscCall(MatSetOption(user.H, MAT_SYMMETRIC, PETSC_TRUE));
516 PetscCall(TaoSetHessian(tao, user.H, user.H, MatrixFreeHessian, (void *)&user));
517 } else { /* Create Hessian matrix */
518 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.H));
519 PetscCall(MatSetSizes(user.H, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
520 PetscCall(MatSetUp(user.H));
521 PetscCall(TaoSetHessian(tao, user.H, user.H, FormHessian, (void *)&user));
522 }
523
524 /* Not use any preconditioner */
525 PetscCall(TaoGetKSP(tao, &ksp));
526 if (ksp) {
527 PetscCall(KSPGetPC(ksp, &pc));
528 PetscCall(PCSetType(pc, PCNONE));
529 }
530
531 /* Set initial solution guess */
532 PetscCall(TaoSetSolution(tao, x));
533 PetscCall(TaoSetFromOptions(tao));
534 PetscCall(TaoSolve(tao));
535 PetscCall(TaoDestroy(&tao));
536 PetscCall(MatDestroy(&user.H));
537 break;
538 default:
539 break;
540 }
541
542 /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */
543 PetscCall(MatDestroy(&user.A));
544 PetscCall(VecDestroy(&user.U));
545 PetscCall(VecDestroy(&user.Lambda[0]));
546 PetscCall(VecDestroy(&user.Lambda2[0]));
547 PetscCall(VecDestroy(&user.Ihp1[0]));
548 PetscCall(VecDestroy(&user.Dir));
549 PetscCall(TSDestroy(&user.ts));
550 PetscCall(VecDestroy(&x));
551 PetscCall(PetscFinalize());
552 return 0;
553 }
554
555 /*TEST
556 build:
557 requires: !complex !single
558
559 test:
560 args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_time_step 0.03125
561 output_file: output/ex20opt_ic_1.out
562
563 test:
564 suffix: 2
565 args: -ts_type beuler -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_time_step 0.01 -tao_type bntr -tao_bnk_pc_type none
566 output_file: output/ex20opt_ic_2.out
567
568 test:
569 suffix: 3
570 args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_time_step 0.01 -tao_type bntr -tao_bnk_pc_type none
571 output_file: output/ex20opt_ic_3.out
572
573 test:
574 suffix: 4
575 args: -implicitform 0 -ts_time_step 0.01 -ts_max_steps 2 -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view -mode 1 -my_tao_mf
576 TEST*/
577