Actual source code: ex1.c
1: #include <petsctao.h>
2: #include <petscts.h>
4: typedef struct _n_aircraft *Aircraft;
5: struct _n_aircraft {
6: TS ts, quadts;
7: Vec V, W; /* control variables V and W */
8: PetscInt nsteps; /* number of time steps */
9: PetscReal ftime;
10: Mat A, H;
11: Mat Jacp, DRDU, DRDP;
12: Vec U, Lambda[1], Mup[1], Lambda2[1], Mup2[1], Dir;
13: Vec rhshp1[1], rhshp2[1], rhshp3[1], rhshp4[1], inthp1[1], inthp2[1], inthp3[1], inthp4[1];
14: PetscReal lv, lw;
15: PetscBool mf, eh;
16: };
18: PetscErrorCode FormObjFunctionGradient(Tao, Vec, PetscReal *, Vec, void *);
19: PetscErrorCode FormObjHessian(Tao, Vec, Mat, Mat, void *);
20: PetscErrorCode ComputeObjHessianWithSOA(Vec, PetscScalar[], Aircraft);
21: PetscErrorCode MatrixFreeObjHessian(Tao, Vec, Mat, Mat, void *);
22: PetscErrorCode MyMatMult(Mat, Vec, Vec);
24: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, PetscCtx ctx)
25: {
26: Aircraft actx = (Aircraft)ctx;
27: const PetscScalar *u, *v, *w;
28: PetscScalar *f;
29: PetscInt step;
31: PetscFunctionBeginUser;
32: PetscCall(TSGetStepNumber(ts, &step));
33: PetscCall(VecGetArrayRead(U, &u));
34: PetscCall(VecGetArrayRead(actx->V, &v));
35: PetscCall(VecGetArrayRead(actx->W, &w));
36: PetscCall(VecGetArray(F, &f));
37: f[0] = v[step] * PetscCosReal(w[step]);
38: f[1] = v[step] * PetscSinReal(w[step]);
39: PetscCall(VecRestoreArrayRead(U, &u));
40: PetscCall(VecRestoreArrayRead(actx->V, &v));
41: PetscCall(VecRestoreArrayRead(actx->W, &w));
42: PetscCall(VecRestoreArray(F, &f));
43: PetscFunctionReturn(PETSC_SUCCESS);
44: }
46: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec U, Mat A, PetscCtx ctx)
47: {
48: Aircraft actx = (Aircraft)ctx;
49: const PetscScalar *u, *v, *w;
50: PetscInt step, rows[2] = {0, 1}, rowcol[2];
51: PetscScalar Jp[2][2];
53: PetscFunctionBeginUser;
54: PetscCall(MatZeroEntries(A));
55: PetscCall(TSGetStepNumber(ts, &step));
56: PetscCall(VecGetArrayRead(U, &u));
57: PetscCall(VecGetArrayRead(actx->V, &v));
58: PetscCall(VecGetArrayRead(actx->W, &w));
60: Jp[0][0] = PetscCosReal(w[step]);
61: Jp[0][1] = -v[step] * PetscSinReal(w[step]);
62: Jp[1][0] = PetscSinReal(w[step]);
63: Jp[1][1] = v[step] * PetscCosReal(w[step]);
65: PetscCall(VecRestoreArrayRead(U, &u));
66: PetscCall(VecRestoreArrayRead(actx->V, &v));
67: PetscCall(VecRestoreArrayRead(actx->W, &w));
69: rowcol[0] = 2 * step;
70: rowcol[1] = 2 * step + 1;
71: PetscCall(MatSetValues(A, 2, rows, 2, rowcol, &Jp[0][0], INSERT_VALUES));
73: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
74: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode RHSHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
79: {
80: PetscFunctionBeginUser;
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: static PetscErrorCode RHSHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
85: {
86: PetscFunctionBeginUser;
87: PetscFunctionReturn(PETSC_SUCCESS);
88: }
90: static PetscErrorCode RHSHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
91: {
92: PetscFunctionBeginUser;
93: PetscFunctionReturn(PETSC_SUCCESS);
94: }
96: static PetscErrorCode RHSHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
97: {
98: Aircraft actx = (Aircraft)ctx;
99: const PetscScalar *v, *w, *vl, *vr, *u;
100: PetscScalar *vhv;
101: PetscScalar dJpdP[2][2][2] = {{{0}}};
102: PetscInt step, i, j, k;
104: PetscFunctionBeginUser;
105: PetscCall(TSGetStepNumber(ts, &step));
106: PetscCall(VecGetArrayRead(U, &u));
107: PetscCall(VecGetArrayRead(actx->V, &v));
108: PetscCall(VecGetArrayRead(actx->W, &w));
109: PetscCall(VecGetArrayRead(Vl[0], &vl));
110: PetscCall(VecGetArrayRead(Vr, &vr));
111: PetscCall(VecSet(VHV[0], 0.0));
112: PetscCall(VecGetArray(VHV[0], &vhv));
114: dJpdP[0][0][1] = -PetscSinReal(w[step]);
115: dJpdP[0][1][0] = -PetscSinReal(w[step]);
116: dJpdP[0][1][1] = -v[step] * PetscCosReal(w[step]);
117: dJpdP[1][0][1] = PetscCosReal(w[step]);
118: dJpdP[1][1][0] = PetscCosReal(w[step]);
119: dJpdP[1][1][1] = -v[step] * PetscSinReal(w[step]);
121: for (j = 0; j < 2; j++) {
122: vhv[2 * step + j] = 0;
123: for (k = 0; k < 2; k++)
124: for (i = 0; i < 2; i++) vhv[2 * step + j] += vl[i] * dJpdP[i][j][k] * vr[2 * step + k];
125: }
126: PetscCall(VecRestoreArrayRead(U, &u));
127: PetscCall(VecRestoreArrayRead(Vl[0], &vl));
128: PetscCall(VecRestoreArrayRead(Vr, &vr));
129: PetscCall(VecRestoreArray(VHV[0], &vhv));
130: PetscFunctionReturn(PETSC_SUCCESS);
131: }
133: /* Vl in NULL,updates to VHV must be added */
134: static PetscErrorCode IntegrandHessianProductUU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
135: {
136: Aircraft actx = (Aircraft)ctx;
137: const PetscScalar *v, *w, *vr, *u;
138: PetscScalar *vhv;
139: PetscScalar dRudU[2][2] = {{0}};
140: PetscInt step, j, k;
142: PetscFunctionBeginUser;
143: PetscCall(TSGetStepNumber(ts, &step));
144: PetscCall(VecGetArrayRead(U, &u));
145: PetscCall(VecGetArrayRead(actx->V, &v));
146: PetscCall(VecGetArrayRead(actx->W, &w));
147: PetscCall(VecGetArrayRead(Vr, &vr));
148: PetscCall(VecGetArray(VHV[0], &vhv));
150: dRudU[0][0] = 2.0;
151: dRudU[1][1] = 2.0;
153: for (j = 0; j < 2; j++) {
154: vhv[j] = 0;
155: for (k = 0; k < 2; k++) vhv[j] += dRudU[j][k] * vr[k];
156: }
157: PetscCall(VecRestoreArrayRead(U, &u));
158: PetscCall(VecRestoreArrayRead(Vr, &vr));
159: PetscCall(VecRestoreArray(VHV[0], &vhv));
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode IntegrandHessianProductUP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
164: {
165: PetscFunctionBeginUser;
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode IntegrandHessianProductPU(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
170: {
171: PetscFunctionBeginUser;
172: PetscFunctionReturn(PETSC_SUCCESS);
173: }
175: static PetscErrorCode IntegrandHessianProductPP(TS ts, PetscReal t, Vec U, Vec *Vl, Vec Vr, Vec *VHV, PetscCtx ctx)
176: {
177: PetscFunctionBeginUser;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, PetscCtx ctx)
182: {
183: Aircraft actx = (Aircraft)ctx;
184: PetscScalar *r;
185: PetscReal dx, dy;
186: const PetscScalar *u;
188: PetscFunctionBegin;
189: PetscCall(VecGetArrayRead(U, &u));
190: PetscCall(VecGetArray(R, &r));
191: dx = u[0] - actx->lv * t * PetscCosReal(actx->lw);
192: dy = u[1] - actx->lv * t * PetscSinReal(actx->lw);
193: r[0] = dx * dx + dy * dy;
194: PetscCall(VecRestoreArray(R, &r));
195: PetscCall(VecRestoreArrayRead(U, &u));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, PetscCtx ctx)
200: {
201: Aircraft actx = (Aircraft)ctx;
202: PetscScalar drdu[2][1];
203: const PetscScalar *u;
204: PetscReal dx, dy;
205: PetscInt row[] = {0, 1}, col[] = {0};
207: PetscFunctionBegin;
208: PetscCall(VecGetArrayRead(U, &u));
209: dx = u[0] - actx->lv * t * PetscCosReal(actx->lw);
210: dy = u[1] - actx->lv * t * PetscSinReal(actx->lw);
211: drdu[0][0] = 2. * dx;
212: drdu[1][0] = 2. * dy;
213: PetscCall(MatSetValues(DRDU, 2, row, 1, col, &drdu[0][0], INSERT_VALUES));
214: PetscCall(VecRestoreArrayRead(U, &u));
215: PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
216: PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, PetscCtx ctx)
221: {
222: PetscFunctionBegin;
223: PetscCall(MatZeroEntries(DRDP));
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: int main(int argc, char **argv)
228: {
229: Vec P, PL, PU;
230: struct _n_aircraft aircraft;
231: PetscMPIInt size;
232: Tao tao;
233: KSP ksp;
234: PC pc;
235: PetscScalar *u, *p;
237: PetscFunctionBeginUser;
238: /* Initialize program */
239: PetscCall(PetscInitialize(&argc, &argv, NULL, NULL));
240: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
241: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!");
243: /* Parameter settings */
244: aircraft.ftime = 1.; /* time interval in hour */
245: aircraft.nsteps = 10; /* number of steps */
246: aircraft.lv = 2.0; /* leader speed in kmph */
247: aircraft.lw = PETSC_PI / 4.; /* leader heading angle */
249: PetscCall(PetscOptionsGetReal(NULL, NULL, "-ftime", &aircraft.ftime, NULL));
250: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nsteps", &aircraft.nsteps, NULL));
251: PetscCall(PetscOptionsHasName(NULL, NULL, "-matrixfree", &aircraft.mf));
252: PetscCall(PetscOptionsHasName(NULL, NULL, "-exacthessian", &aircraft.eh));
254: /* Create TAO solver and set desired solution method */
255: PetscCall(TaoCreate(PETSC_COMM_WORLD, &tao));
256: PetscCall(TaoSetType(tao, TAOBQNLS));
258: /* Create necessary matrix and vectors, solve same ODE on every process */
259: PetscCall(MatCreate(PETSC_COMM_WORLD, &aircraft.A));
260: PetscCall(MatSetSizes(aircraft.A, PETSC_DECIDE, PETSC_DECIDE, 2, 2));
261: PetscCall(MatSetFromOptions(aircraft.A));
262: PetscCall(MatSetUp(aircraft.A));
263: /* this is to set explicit zeros along the diagonal of the matrix */
264: PetscCall(MatAssemblyBegin(aircraft.A, MAT_FINAL_ASSEMBLY));
265: PetscCall(MatAssemblyEnd(aircraft.A, MAT_FINAL_ASSEMBLY));
266: PetscCall(MatShift(aircraft.A, 1));
267: PetscCall(MatShift(aircraft.A, -1));
269: PetscCall(MatCreate(PETSC_COMM_WORLD, &aircraft.Jacp));
270: PetscCall(MatSetSizes(aircraft.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 2 * aircraft.nsteps));
271: PetscCall(MatSetFromOptions(aircraft.Jacp));
272: PetscCall(MatSetUp(aircraft.Jacp));
273: PetscCall(MatSetOption(aircraft.Jacp, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
274: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2 * aircraft.nsteps, 1, NULL, &aircraft.DRDP));
275: PetscCall(MatSetUp(aircraft.DRDP));
276: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &aircraft.DRDU));
277: PetscCall(MatSetUp(aircraft.DRDU));
279: /* Create timestepping solver context */
280: PetscCall(TSCreate(PETSC_COMM_WORLD, &aircraft.ts));
281: PetscCall(TSSetType(aircraft.ts, TSRK));
282: PetscCall(TSSetRHSFunction(aircraft.ts, NULL, RHSFunction, &aircraft));
283: PetscCall(TSSetRHSJacobian(aircraft.ts, aircraft.A, aircraft.A, TSComputeRHSJacobianConstant, &aircraft));
284: PetscCall(TSSetRHSJacobianP(aircraft.ts, aircraft.Jacp, RHSJacobianP, &aircraft));
285: PetscCall(TSSetExactFinalTime(aircraft.ts, TS_EXACTFINALTIME_MATCHSTEP));
286: PetscCall(TSSetEquationType(aircraft.ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
288: /* Set initial conditions */
289: PetscCall(MatCreateVecs(aircraft.A, &aircraft.U, NULL));
290: PetscCall(TSSetSolution(aircraft.ts, aircraft.U));
291: PetscCall(VecGetArray(aircraft.U, &u));
292: u[0] = 1.5;
293: u[1] = 0;
294: PetscCall(VecRestoreArray(aircraft.U, &u));
295: PetscCall(VecCreate(PETSC_COMM_WORLD, &aircraft.V));
296: PetscCall(VecSetSizes(aircraft.V, PETSC_DECIDE, aircraft.nsteps));
297: PetscCall(VecSetUp(aircraft.V));
298: PetscCall(VecDuplicate(aircraft.V, &aircraft.W));
299: PetscCall(VecSet(aircraft.V, 1.));
300: PetscCall(VecSet(aircraft.W, PETSC_PI / 4.));
302: /* Save trajectory of solution so that TSAdjointSolve() may be used */
303: PetscCall(TSSetSaveTrajectory(aircraft.ts));
305: /* Set sensitivity context */
306: PetscCall(TSCreateQuadratureTS(aircraft.ts, PETSC_FALSE, &aircraft.quadts));
307: PetscCall(TSSetRHSFunction(aircraft.quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &aircraft));
308: PetscCall(TSSetRHSJacobian(aircraft.quadts, aircraft.DRDU, aircraft.DRDU, (TSRHSJacobian)DRDUJacobianTranspose, &aircraft));
309: PetscCall(TSSetRHSJacobianP(aircraft.quadts, aircraft.DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, &aircraft));
310: PetscCall(MatCreateVecs(aircraft.A, &aircraft.Lambda[0], NULL));
311: PetscCall(MatCreateVecs(aircraft.Jacp, &aircraft.Mup[0], NULL));
312: if (aircraft.eh) {
313: PetscCall(MatCreateVecs(aircraft.A, &aircraft.rhshp1[0], NULL));
314: PetscCall(MatCreateVecs(aircraft.A, &aircraft.rhshp2[0], NULL));
315: PetscCall(MatCreateVecs(aircraft.Jacp, &aircraft.rhshp3[0], NULL));
316: PetscCall(MatCreateVecs(aircraft.Jacp, &aircraft.rhshp4[0], NULL));
317: PetscCall(MatCreateVecs(aircraft.DRDU, &aircraft.inthp1[0], NULL));
318: PetscCall(MatCreateVecs(aircraft.DRDU, &aircraft.inthp2[0], NULL));
319: PetscCall(MatCreateVecs(aircraft.DRDP, &aircraft.inthp3[0], NULL));
320: PetscCall(MatCreateVecs(aircraft.DRDP, &aircraft.inthp4[0], NULL));
321: PetscCall(MatCreateVecs(aircraft.Jacp, &aircraft.Dir, NULL));
322: PetscCall(TSSetRHSHessianProduct(aircraft.ts, aircraft.rhshp1, RHSHessianProductUU, aircraft.rhshp2, RHSHessianProductUP, aircraft.rhshp3, RHSHessianProductPU, aircraft.rhshp4, RHSHessianProductPP, &aircraft));
323: PetscCall(TSSetRHSHessianProduct(aircraft.quadts, aircraft.inthp1, IntegrandHessianProductUU, aircraft.inthp2, IntegrandHessianProductUP, aircraft.inthp3, IntegrandHessianProductPU, aircraft.inthp4, IntegrandHessianProductPP, &aircraft));
324: PetscCall(MatCreateVecs(aircraft.A, &aircraft.Lambda2[0], NULL));
325: PetscCall(MatCreateVecs(aircraft.Jacp, &aircraft.Mup2[0], NULL));
326: }
327: PetscCall(TSSetFromOptions(aircraft.ts));
328: PetscCall(TSSetMaxTime(aircraft.ts, aircraft.ftime));
329: PetscCall(TSSetTimeStep(aircraft.ts, aircraft.ftime / aircraft.nsteps));
331: /* Set initial solution guess */
332: PetscCall(MatCreateVecs(aircraft.Jacp, &P, NULL));
333: PetscCall(VecGetArray(P, &p));
334: for (PetscInt i = 0; i < aircraft.nsteps; i++) {
335: p[2 * i] = 2.0;
336: p[2 * i + 1] = PETSC_PI / 2.0;
337: }
338: PetscCall(VecRestoreArray(P, &p));
339: PetscCall(VecDuplicate(P, &PU));
340: PetscCall(VecDuplicate(P, &PL));
341: PetscCall(VecGetArray(PU, &p));
342: for (PetscInt i = 0; i < aircraft.nsteps; i++) {
343: p[2 * i] = 2.0;
344: p[2 * i + 1] = PETSC_PI;
345: }
346: PetscCall(VecRestoreArray(PU, &p));
347: PetscCall(VecGetArray(PL, &p));
348: for (PetscInt i = 0; i < aircraft.nsteps; i++) {
349: p[2 * i] = 0.0;
350: p[2 * i + 1] = -PETSC_PI;
351: }
352: PetscCall(VecRestoreArray(PL, &p));
354: PetscCall(TaoSetSolution(tao, P));
355: PetscCall(TaoSetVariableBounds(tao, PL, PU));
356: /* Set routine for function and gradient evaluation */
357: PetscCall(TaoSetObjectiveAndGradient(tao, NULL, FormObjFunctionGradient, (void *)&aircraft));
359: if (aircraft.eh) {
360: if (aircraft.mf) {
361: PetscCall(MatCreateShell(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2 * aircraft.nsteps, 2 * aircraft.nsteps, (void *)&aircraft, &aircraft.H));
362: PetscCall(MatShellSetOperation(aircraft.H, MATOP_MULT, (PetscErrorCodeFn *)MyMatMult));
363: PetscCall(MatSetOption(aircraft.H, MAT_SYMMETRIC, PETSC_TRUE));
364: PetscCall(TaoSetHessian(tao, aircraft.H, aircraft.H, MatrixFreeObjHessian, (void *)&aircraft));
365: } else {
366: PetscCall(MatCreateDense(MPI_COMM_WORLD, PETSC_DETERMINE, PETSC_DETERMINE, 2 * aircraft.nsteps, 2 * aircraft.nsteps, NULL, &aircraft.H));
367: PetscCall(MatSetOption(aircraft.H, MAT_SYMMETRIC, PETSC_TRUE));
368: PetscCall(TaoSetHessian(tao, aircraft.H, aircraft.H, FormObjHessian, (void *)&aircraft));
369: }
370: }
372: /* Check for any TAO command line options */
373: PetscCall(TaoGetKSP(tao, &ksp));
374: if (ksp) {
375: PetscCall(KSPGetPC(ksp, &pc));
376: PetscCall(PCSetType(pc, PCNONE));
377: }
378: PetscCall(TaoSetFromOptions(tao));
380: PetscCall(TaoSolve(tao));
381: PetscCall(VecView(P, PETSC_VIEWER_STDOUT_WORLD));
383: /* Free TAO data structures */
384: PetscCall(TaoDestroy(&tao));
386: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
387: Free work space. All PETSc objects should be destroyed when they
388: are no longer needed.
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: PetscCall(TSDestroy(&aircraft.ts));
391: PetscCall(MatDestroy(&aircraft.A));
392: PetscCall(VecDestroy(&aircraft.U));
393: PetscCall(VecDestroy(&aircraft.V));
394: PetscCall(VecDestroy(&aircraft.W));
395: PetscCall(VecDestroy(&P));
396: PetscCall(VecDestroy(&PU));
397: PetscCall(VecDestroy(&PL));
398: PetscCall(MatDestroy(&aircraft.Jacp));
399: PetscCall(MatDestroy(&aircraft.DRDU));
400: PetscCall(MatDestroy(&aircraft.DRDP));
401: PetscCall(VecDestroy(&aircraft.Lambda[0]));
402: PetscCall(VecDestroy(&aircraft.Mup[0]));
403: PetscCall(VecDestroy(&P));
404: if (aircraft.eh) {
405: PetscCall(VecDestroy(&aircraft.Lambda2[0]));
406: PetscCall(VecDestroy(&aircraft.Mup2[0]));
407: PetscCall(VecDestroy(&aircraft.Dir));
408: PetscCall(VecDestroy(&aircraft.rhshp1[0]));
409: PetscCall(VecDestroy(&aircraft.rhshp2[0]));
410: PetscCall(VecDestroy(&aircraft.rhshp3[0]));
411: PetscCall(VecDestroy(&aircraft.rhshp4[0]));
412: PetscCall(VecDestroy(&aircraft.inthp1[0]));
413: PetscCall(VecDestroy(&aircraft.inthp2[0]));
414: PetscCall(VecDestroy(&aircraft.inthp3[0]));
415: PetscCall(VecDestroy(&aircraft.inthp4[0]));
416: PetscCall(MatDestroy(&aircraft.H));
417: }
418: PetscCall(PetscFinalize());
419: return 0;
420: }
422: /*
423: FormObjFunctionGradient - Evaluates the function and corresponding gradient.
425: Input Parameters:
426: tao - the Tao context
427: P - the input vector
428: ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradient()
430: Output Parameters:
431: f - the newly evaluated function
432: G - the newly evaluated gradient
433: */
434: PetscErrorCode FormObjFunctionGradient(Tao tao, Vec P, PetscReal *f, Vec G, PetscCtx ctx)
435: {
436: Aircraft actx = (Aircraft)ctx;
437: TS ts = actx->ts;
438: Vec Q;
439: const PetscScalar *p, *q;
440: PetscScalar *u, *v, *w;
441: PetscInt i;
443: PetscFunctionBeginUser;
444: PetscCall(VecGetArrayRead(P, &p));
445: PetscCall(VecGetArray(actx->V, &v));
446: PetscCall(VecGetArray(actx->W, &w));
447: for (i = 0; i < actx->nsteps; i++) {
448: v[i] = p[2 * i];
449: w[i] = p[2 * i + 1];
450: }
451: PetscCall(VecRestoreArrayRead(P, &p));
452: PetscCall(VecRestoreArray(actx->V, &v));
453: PetscCall(VecRestoreArray(actx->W, &w));
455: PetscCall(TSSetTime(ts, 0.0));
456: PetscCall(TSSetStepNumber(ts, 0));
457: PetscCall(TSSetFromOptions(ts));
458: PetscCall(TSSetTimeStep(ts, actx->ftime / actx->nsteps));
460: /* reinitialize system state */
461: PetscCall(VecGetArray(actx->U, &u));
462: u[0] = 2.0;
463: u[1] = 0;
464: PetscCall(VecRestoreArray(actx->U, &u));
466: /* reinitialize the integral value */
467: PetscCall(TSGetCostIntegral(ts, &Q));
468: PetscCall(VecSet(Q, 0.0));
470: PetscCall(TSSolve(ts, actx->U));
472: /* Reset initial conditions for the adjoint integration */
473: PetscCall(VecSet(actx->Lambda[0], 0.0));
474: PetscCall(VecSet(actx->Mup[0], 0.0));
475: PetscCall(TSSetCostGradients(ts, 1, actx->Lambda, actx->Mup));
477: PetscCall(TSAdjointSolve(ts));
478: PetscCall(VecCopy(actx->Mup[0], G));
479: PetscCall(TSGetCostIntegral(ts, &Q));
480: PetscCall(VecGetArrayRead(Q, &q));
481: *f = q[0];
482: PetscCall(VecRestoreArrayRead(Q, &q));
483: PetscFunctionReturn(PETSC_SUCCESS);
484: }
486: PetscErrorCode FormObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, PetscCtx ctx)
487: {
488: Aircraft actx = (Aircraft)ctx;
489: const PetscScalar *p;
490: PetscScalar *harr, *v, *w, one = 1.0;
491: PetscInt ind[1];
492: PetscInt *cols, i;
493: Vec Dir;
495: PetscFunctionBeginUser;
496: /* set up control parameters */
497: PetscCall(VecGetArrayRead(P, &p));
498: PetscCall(VecGetArray(actx->V, &v));
499: PetscCall(VecGetArray(actx->W, &w));
500: for (i = 0; i < actx->nsteps; i++) {
501: v[i] = p[2 * i];
502: w[i] = p[2 * i + 1];
503: }
504: PetscCall(VecRestoreArrayRead(P, &p));
505: PetscCall(VecRestoreArray(actx->V, &v));
506: PetscCall(VecRestoreArray(actx->W, &w));
508: PetscCall(PetscMalloc1(2 * actx->nsteps, &harr));
509: PetscCall(PetscMalloc1(2 * actx->nsteps, &cols));
510: for (i = 0; i < 2 * actx->nsteps; i++) cols[i] = i;
511: PetscCall(VecDuplicate(P, &Dir));
512: for (i = 0; i < 2 * actx->nsteps; i++) {
513: ind[0] = i;
514: PetscCall(VecSet(Dir, 0.0));
515: PetscCall(VecSetValues(Dir, 1, ind, &one, INSERT_VALUES));
516: PetscCall(VecAssemblyBegin(Dir));
517: PetscCall(VecAssemblyEnd(Dir));
518: PetscCall(ComputeObjHessianWithSOA(Dir, harr, actx));
519: PetscCall(MatSetValues(H, 1, ind, 2 * actx->nsteps, cols, harr, INSERT_VALUES));
520: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
521: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
522: if (H != Hpre) {
523: PetscCall(MatAssemblyBegin(Hpre, MAT_FINAL_ASSEMBLY));
524: PetscCall(MatAssemblyEnd(Hpre, MAT_FINAL_ASSEMBLY));
525: }
526: }
527: PetscCall(PetscFree(cols));
528: PetscCall(PetscFree(harr));
529: PetscCall(VecDestroy(&Dir));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, PetscCtx ctx)
534: {
535: Aircraft actx = (Aircraft)ctx;
536: PetscScalar *v, *w;
537: const PetscScalar *p;
539: PetscFunctionBegin;
540: PetscCall(VecGetArrayRead(P, &p));
541: PetscCall(VecGetArray(actx->V, &v));
542: PetscCall(VecGetArray(actx->W, &w));
543: for (PetscInt i = 0; i < actx->nsteps; i++) {
544: v[i] = p[2 * i];
545: w[i] = p[2 * i + 1];
546: }
547: PetscCall(VecRestoreArrayRead(P, &p));
548: PetscCall(VecRestoreArray(actx->V, &v));
549: PetscCall(VecRestoreArray(actx->W, &w));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
554: {
555: PetscScalar *y;
556: void *ptr;
558: PetscFunctionBegin;
559: PetscCall(MatShellGetContext(H_shell, &ptr));
560: PetscCall(VecGetArray(Y, &y));
561: PetscCall(ComputeObjHessianWithSOA(X, y, (Aircraft)ptr));
562: PetscCall(VecRestoreArray(Y, &y));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: PetscErrorCode ComputeObjHessianWithSOA(Vec Dir, PetscScalar arr[], Aircraft actx)
567: {
568: TS ts = actx->ts;
569: const PetscScalar *z_ptr;
570: PetscScalar *u;
571: Vec Q;
573: PetscFunctionBeginUser;
574: /* Reset TSAdjoint so that AdjointSetUp will be called again */
575: PetscCall(TSAdjointReset(ts));
577: PetscCall(TSSetTime(ts, 0.0));
578: PetscCall(TSSetStepNumber(ts, 0));
579: PetscCall(TSSetFromOptions(ts));
580: PetscCall(TSSetTimeStep(ts, actx->ftime / actx->nsteps));
581: PetscCall(TSSetCostHessianProducts(actx->ts, 1, actx->Lambda2, actx->Mup2, Dir));
583: /* reinitialize system state */
584: PetscCall(VecGetArray(actx->U, &u));
585: u[0] = 2.0;
586: u[1] = 0;
587: PetscCall(VecRestoreArray(actx->U, &u));
589: /* reinitialize the integral value */
590: PetscCall(TSGetCostIntegral(ts, &Q));
591: PetscCall(VecSet(Q, 0.0));
593: /* initialize tlm variable */
594: PetscCall(MatZeroEntries(actx->Jacp));
595: PetscCall(TSAdjointSetForward(ts, actx->Jacp));
597: PetscCall(TSSolve(ts, actx->U));
599: /* Set terminal conditions for first- and second-order adjonts */
600: PetscCall(VecSet(actx->Lambda[0], 0.0));
601: PetscCall(VecSet(actx->Mup[0], 0.0));
602: PetscCall(VecSet(actx->Lambda2[0], 0.0));
603: PetscCall(VecSet(actx->Mup2[0], 0.0));
604: PetscCall(TSSetCostGradients(ts, 1, actx->Lambda, actx->Mup));
606: PetscCall(TSGetCostIntegral(ts, &Q));
608: /* Reset initial conditions for the adjoint integration */
609: PetscCall(TSAdjointSolve(ts));
611: /* initial condition does not depend on p, so that lambda is not needed to assemble G */
612: PetscCall(VecGetArrayRead(actx->Mup2[0], &z_ptr));
613: for (PetscInt i = 0; i < 2 * actx->nsteps; i++) arr[i] = z_ptr[i];
614: PetscCall(VecRestoreArrayRead(actx->Mup2[0], &z_ptr));
616: /* Disable second-order adjoint mode */
617: PetscCall(TSAdjointReset(ts));
618: PetscCall(TSAdjointResetForward(ts));
619: PetscFunctionReturn(PETSC_SUCCESS);
620: }
622: /*TEST
623: build:
624: requires: !complex !single
626: test:
627: args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_gatol 1e-7
629: test:
630: suffix: 2
631: args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian
633: test:
634: suffix: 3
635: args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_view -tao_type bntr -tao_bnk_pc_type none -exacthessian -matrixfree
636: TEST*/