Actual source code: ex1.c
petsc-3.13.4 2020-08-01
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,void *ctx)
25: {
26: Aircraft actx = (Aircraft)ctx;
27: const PetscScalar *u,*v,*w;
28: PetscScalar *f;
29: PetscInt step;
30: PetscErrorCode ierr;
33: TSGetStepNumber(ts,&step);
34: VecGetArrayRead(U,&u);
35: VecGetArrayRead(actx->V,&v);
36: VecGetArrayRead(actx->W,&w);
37: VecGetArray(F,&f);
38: f[0] = v[step]*PetscCosReal(w[step]);
39: f[1] = v[step]*PetscSinReal(w[step]);
40: VecRestoreArrayRead(U,&u);
41: VecRestoreArrayRead(actx->V,&v);
42: VecRestoreArrayRead(actx->W,&w);
43: VecRestoreArray(F,&f);
44: return(0);
45: }
47: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx)
48: {
49: Aircraft actx = (Aircraft)ctx;
50: const PetscScalar *u,*v,*w;
51: PetscInt step,rows[2] = {0,1},rowcol[2];
52: PetscScalar Jp[2][2];
53: PetscErrorCode ierr;
56: MatZeroEntries(A);
57: TSGetStepNumber(ts,&step);
58: VecGetArrayRead(U,&u);
59: VecGetArrayRead(actx->V,&v);
60: VecGetArrayRead(actx->W,&w);
62: Jp[0][0] = PetscCosReal(w[step]);
63: Jp[0][1] = -v[step]*PetscSinReal(w[step]);
64: Jp[1][0] = PetscSinReal(w[step]);
65: Jp[1][1] = v[step]*PetscCosReal(w[step]);
67: VecRestoreArrayRead(U,&u);
68: VecRestoreArrayRead(actx->V,&v);
69: VecRestoreArrayRead(actx->W,&w);
71: rowcol[0] = 2*step;
72: rowcol[1] = 2*step+1;
73: MatSetValues(A,2,rows,2,rowcol,&Jp[0][0],INSERT_VALUES);
75: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
76: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
77: return(0);
78: }
80: static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
81: {
83: return(0);
84: }
86: static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
87: {
89: return(0);
90: }
92: static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
93: {
95: return(0);
96: }
98: static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
99: {
100: Aircraft actx = (Aircraft)ctx;
101: const PetscScalar *v,*w,*vl,*vr,*u;
102: PetscScalar *vhv;
103: PetscScalar dJpdP[2][2][2]={{{0}}};
104: PetscInt step,i,j,k;
105: PetscErrorCode ierr;
108: TSGetStepNumber(ts,&step);
109: VecGetArrayRead(U,&u);
110: VecGetArrayRead(actx->V,&v);
111: VecGetArrayRead(actx->W,&w);
112: VecGetArrayRead(Vl[0],&vl);
113: VecGetArrayRead(Vr,&vr);
114: VecSet(VHV[0],0.0);
115: VecGetArray(VHV[0],&vhv);
117: dJpdP[0][0][1] = -PetscSinReal(w[step]);
118: dJpdP[0][1][0] = -PetscSinReal(w[step]);
119: dJpdP[0][1][1] = -v[step]*PetscCosReal(w[step]);
120: dJpdP[1][0][1] = PetscCosReal(w[step]);
121: dJpdP[1][1][0] = PetscCosReal(w[step]);
122: dJpdP[1][1][1] = -v[step]*PetscSinReal(w[step]);
124: for (j=0; j<2; j++) {
125: vhv[2*step+j] = 0;
126: for (k=0; k<2; k++)
127: for (i=0; i<2; i++)
128: vhv[2*step+j] += vl[i]*dJpdP[i][j][k]*vr[2*step+k];
129: }
130: VecRestoreArrayRead(U,&u);
131: VecRestoreArrayRead(Vl[0],&vl);
132: VecRestoreArrayRead(Vr,&vr);
133: VecRestoreArray(VHV[0],&vhv);
134: return(0);
135: }
137: /* Vl in NULL,updates to VHV must be added */
138: static PetscErrorCode IntegrandHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
139: {
140: Aircraft actx = (Aircraft)ctx;
141: const PetscScalar *v,*w,*vr,*u;
142: PetscScalar *vhv;
143: PetscScalar dRudU[2][2]={{0}};
144: PetscInt step,j,k;
145: PetscErrorCode ierr;
148: TSGetStepNumber(ts,&step);
149: VecGetArrayRead(U,&u);
150: VecGetArrayRead(actx->V,&v);
151: VecGetArrayRead(actx->W,&w);
152: VecGetArrayRead(Vr,&vr);
153: VecGetArray(VHV[0],&vhv);
155: dRudU[0][0] = 2.0;
156: dRudU[1][1] = 2.0;
158: for (j=0; j<2; j++) {
159: vhv[j] = 0;
160: for (k=0; k<2; k++)
161: vhv[j] += dRudU[j][k]*vr[k];
162: }
163: VecRestoreArrayRead(U,&u);
164: VecRestoreArrayRead(Vr,&vr);
165: VecRestoreArray(VHV[0],&vhv);
166: return(0);
167: }
169: static PetscErrorCode IntegrandHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
170: {
172: return(0);
173: }
175: static PetscErrorCode IntegrandHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
176: {
178: return(0);
179: }
181: static PetscErrorCode IntegrandHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx)
182: {
184: return(0);
185: }
188: static PetscErrorCode CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,void *ctx)
189: {
190: Aircraft actx = (Aircraft)ctx;
191: PetscScalar *r;
192: PetscReal dx,dy;
193: const PetscScalar *u;
194: PetscErrorCode ierr;
197: VecGetArrayRead(U,&u);
198: VecGetArray(R,&r);
199: dx = u[0] - actx->lv*t*PetscCosReal(actx->lw);
200: dy = u[1] - actx->lv*t*PetscSinReal(actx->lw);
201: r[0] = dx*dx+dy*dy;
202: VecRestoreArray(R,&r);
203: VecRestoreArrayRead(U,&u);
204: return(0);
205: }
207: static PetscErrorCode DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,void *ctx)
208: {
209: Aircraft actx = (Aircraft)ctx;
210: PetscScalar drdu[2][1];
211: const PetscScalar *u;
212: PetscReal dx,dy;
213: PetscInt row[] = {0,1},col[] = {0};
214: PetscErrorCode ierr;
217: VecGetArrayRead(U,&u);
218: dx = u[0] - actx->lv*t*PetscCosReal(actx->lw);
219: dy = u[1] - actx->lv*t*PetscSinReal(actx->lw);
220: drdu[0][0] = 2.*dx;
221: drdu[1][0] = 2.*dy;
222: MatSetValues(DRDU,2,row,1,col,&drdu[0][0],INSERT_VALUES);
223: VecRestoreArrayRead(U,&u);
224: MatAssemblyBegin(DRDU,MAT_FINAL_ASSEMBLY);
225: MatAssemblyEnd(DRDU,MAT_FINAL_ASSEMBLY);
226: return(0);
227: }
229: static PetscErrorCode DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,void *ctx)
230: {
234: MatZeroEntries(DRDP);
235: MatAssemblyBegin(DRDP,MAT_FINAL_ASSEMBLY);
236: MatAssemblyEnd(DRDP,MAT_FINAL_ASSEMBLY);
237: return(0);
238: }
240: int main(int argc,char **argv)
241: {
242: Vec P,PL,PU;
243: struct _n_aircraft aircraft;
244: PetscMPIInt size;
245: Tao tao;
246: KSP ksp;
247: PC pc;
248: PetscScalar *u,*p;
249: PetscInt i;
250: PetscErrorCode ierr;
252: /* Initialize program */
253: PetscInitialize(&argc,&argv,NULL,NULL);if(ierr) return ierr;
254: MPI_Comm_size(PETSC_COMM_WORLD,&size);
255: if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!");
257: /* Parameter settings */
258: aircraft.ftime = 1.; /* time interval in hour */
259: aircraft.nsteps = 10; /* number of steps */
260: aircraft.lv = 2.0; /* leader speed in kmph */
261: aircraft.lw = PETSC_PI/4.; /* leader heading angle */
263: PetscOptionsGetReal(NULL,NULL,"-ftime",&aircraft.ftime,NULL);
264: PetscOptionsGetInt(NULL,NULL,"-nsteps",&aircraft.nsteps,NULL);
265: PetscOptionsHasName(NULL,NULL,"-matrixfree",&aircraft.mf);
266: PetscOptionsHasName(NULL,NULL,"-exacthessian",&aircraft.eh);
268: /* Create TAO solver and set desired solution method */
269: TaoCreate(PETSC_COMM_WORLD,&tao);
270: TaoSetType(tao,TAOBQNLS);
272: /* Create necessary matrix and vectors, solve same ODE on every process */
273: MatCreate(PETSC_COMM_WORLD,&aircraft.A);
274: MatSetSizes(aircraft.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
275: MatSetFromOptions(aircraft.A);
276: MatSetUp(aircraft.A);
277: MatAssemblyBegin(aircraft.A,MAT_FINAL_ASSEMBLY);
278: MatAssemblyEnd(aircraft.A,MAT_FINAL_ASSEMBLY);
279: MatShift(aircraft.A,1);
280: MatShift(aircraft.A,-1);
282: MatCreate(PETSC_COMM_WORLD,&aircraft.Jacp);
283: MatSetSizes(aircraft.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,2*aircraft.nsteps);
284: MatSetFromOptions(aircraft.Jacp);
285: MatSetUp(aircraft.Jacp);
286: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,1,NULL,&aircraft.DRDP);
287: MatSetUp(aircraft.DRDP);
288: MatCreateDense(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,&aircraft.DRDU);
289: MatSetUp(aircraft.DRDU);
291: /* Create timestepping solver context */
292: TSCreate(PETSC_COMM_WORLD,&aircraft.ts);
293: TSSetType(aircraft.ts,TSRK);
294: TSSetRHSFunction(aircraft.ts,NULL,RHSFunction,&aircraft);
295: TSSetRHSJacobian(aircraft.ts,aircraft.A,aircraft.A,TSComputeRHSJacobianConstant,&aircraft);
296: TSSetRHSJacobianP(aircraft.ts,aircraft.Jacp,RHSJacobianP,&aircraft);
297: TSSetExactFinalTime(aircraft.ts,TS_EXACTFINALTIME_MATCHSTEP);
298: TSSetEquationType(aircraft.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
300: /* Set initial conditions */
301: MatCreateVecs(aircraft.A,&aircraft.U,NULL);
302: TSSetSolution(aircraft.ts,aircraft.U);
303: VecGetArray(aircraft.U,&u);
304: u[0] = 1.5;
305: u[1] = 0;
306: VecRestoreArray(aircraft.U,&u);
307: VecCreate(PETSC_COMM_WORLD,&aircraft.V);
308: VecSetSizes(aircraft.V,PETSC_DECIDE,aircraft.nsteps);
309: VecSetUp(aircraft.V);
310: VecDuplicate(aircraft.V,&aircraft.W);
311: VecSet(aircraft.V,1.);
312: VecSet(aircraft.W,PETSC_PI/4.);
314: /* Save trajectory of solution so that TSAdjointSolve() may be used */
315: TSSetSaveTrajectory(aircraft.ts);
317: /* Set sensitivity context */
318: TSCreateQuadratureTS(aircraft.ts,PETSC_FALSE,&aircraft.quadts);
319: TSSetRHSFunction(aircraft.quadts,NULL,(TSRHSFunction)CostIntegrand,&aircraft);
320: TSSetRHSJacobian(aircraft.quadts,aircraft.DRDU,aircraft.DRDU,(TSRHSJacobian)DRDUJacobianTranspose,&aircraft);
321: TSSetRHSJacobianP(aircraft.quadts,aircraft.DRDP,(TSRHSJacobianP)DRDPJacobianTranspose,&aircraft);
322: MatCreateVecs(aircraft.A,&aircraft.Lambda[0],NULL);
323: MatCreateVecs(aircraft.Jacp,&aircraft.Mup[0],NULL);
324: if (aircraft.eh) {
325: MatCreateVecs(aircraft.A,&aircraft.rhshp1[0],NULL);
326: MatCreateVecs(aircraft.A,&aircraft.rhshp2[0],NULL);
327: MatCreateVecs(aircraft.Jacp,&aircraft.rhshp3[0],NULL);
328: MatCreateVecs(aircraft.Jacp,&aircraft.rhshp4[0],NULL);
329: MatCreateVecs(aircraft.DRDU,&aircraft.inthp1[0],NULL);
330: MatCreateVecs(aircraft.DRDU,&aircraft.inthp2[0],NULL);
331: MatCreateVecs(aircraft.DRDP,&aircraft.inthp3[0],NULL);
332: MatCreateVecs(aircraft.DRDP,&aircraft.inthp4[0],NULL);
333: MatCreateVecs(aircraft.Jacp,&aircraft.Dir,NULL);
334: TSSetRHSHessianProduct(aircraft.ts,aircraft.rhshp1,RHSHessianProductUU,aircraft.rhshp2,RHSHessianProductUP,aircraft.rhshp3,RHSHessianProductPU,aircraft.rhshp4,RHSHessianProductPP,&aircraft);
335: TSSetRHSHessianProduct(aircraft.quadts,aircraft.inthp1,IntegrandHessianProductUU,aircraft.inthp2,IntegrandHessianProductUP,aircraft.inthp3,IntegrandHessianProductPU,aircraft.inthp4,IntegrandHessianProductPP,&aircraft);
336: MatCreateVecs(aircraft.A,&aircraft.Lambda2[0],NULL);
337: MatCreateVecs(aircraft.Jacp,&aircraft.Mup2[0],NULL);
338: }
339: TSSetFromOptions(aircraft.ts);
340: TSSetMaxTime(aircraft.ts,aircraft.ftime);
341: TSSetTimeStep(aircraft.ts,aircraft.ftime/aircraft.nsteps);
343: /* Set initial solution guess */
344: MatCreateVecs(aircraft.Jacp,&P,NULL);
345: VecGetArray(P,&p);
346: for (i=0; i<aircraft.nsteps; i++) {
347: p[2*i] = 2.0;
348: p[2*i+1] = PETSC_PI/2.0;
349: }
350: VecRestoreArray(P,&p);
351: VecDuplicate(P,&PU);
352: VecDuplicate(P,&PL);
353: VecGetArray(PU,&p);
354: for (i=0; i<aircraft.nsteps; i++) {
355: p[2*i] = 2.0;
356: p[2*i+1] = PETSC_PI;
357: }
358: VecRestoreArray(PU,&p);
359: VecGetArray(PL,&p);
360: for (i=0; i<aircraft.nsteps; i++) {
361: p[2*i] = 0.0;
362: p[2*i+1] = -PETSC_PI;
363: }
364: VecRestoreArray(PL,&p);
366: TaoSetInitialVector(tao,P);
367: TaoSetVariableBounds(tao,PL,PU);
368: /* Set routine for function and gradient evaluation */
369: TaoSetObjectiveAndGradientRoutine(tao,FormObjFunctionGradient,(void *)&aircraft);
371: if (aircraft.eh) {
372: if (aircraft.mf) {
373: MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,2*aircraft.nsteps,2*aircraft.nsteps,(void*)&aircraft,&aircraft.H);
374: MatShellSetOperation(aircraft.H,MATOP_MULT,(void(*)(void))MyMatMult);
375: MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
376: TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,MatrixFreeObjHessian,(void*)&aircraft);
377: } else {
378: MatCreateDense(MPI_COMM_WORLD,PETSC_DETERMINE,PETSC_DETERMINE,2*aircraft.nsteps,2*aircraft.nsteps,NULL,&(aircraft.H));
379: MatSetOption(aircraft.H,MAT_SYMMETRIC,PETSC_TRUE);
380: TaoSetHessianRoutine(tao,aircraft.H,aircraft.H,FormObjHessian,(void *)&aircraft);
381: }
382: }
384: /* Check for any TAO command line options */
385: TaoGetKSP(tao,&ksp);
386: if (ksp) {
387: KSPGetPC(ksp,&pc);
388: PCSetType(pc,PCNONE);
389: }
390: TaoSetFromOptions(tao);
392: TaoSolve(tao);
393: VecView(P,PETSC_VIEWER_STDOUT_WORLD);
395: /* Free TAO data structures */
396: TaoDestroy(&tao);
398: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
399: Free work space. All PETSc objects should be destroyed when they
400: are no longer needed.
401: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
402: TSDestroy(&aircraft.ts);
403: MatDestroy(&aircraft.A);
404: VecDestroy(&aircraft.U);
405: VecDestroy(&aircraft.V);
406: VecDestroy(&aircraft.W);
407: VecDestroy(&P);
408: VecDestroy(&PU);
409: VecDestroy(&PL);
410: MatDestroy(&aircraft.Jacp);
411: MatDestroy(&aircraft.DRDU);
412: MatDestroy(&aircraft.DRDP);
413: VecDestroy(&aircraft.Lambda[0]);
414: VecDestroy(&aircraft.Mup[0]);
415: VecDestroy(&P);
416: if (aircraft.eh) {
417: VecDestroy(&aircraft.Lambda2[0]);
418: VecDestroy(&aircraft.Mup2[0]);
419: VecDestroy(&aircraft.Dir);
420: VecDestroy(&aircraft.rhshp1[0]);
421: VecDestroy(&aircraft.rhshp2[0]);
422: VecDestroy(&aircraft.rhshp3[0]);
423: VecDestroy(&aircraft.rhshp4[0]);
424: VecDestroy(&aircraft.inthp1[0]);
425: VecDestroy(&aircraft.inthp2[0]);
426: VecDestroy(&aircraft.inthp3[0]);
427: VecDestroy(&aircraft.inthp4[0]);
428: MatDestroy(&aircraft.H);
429: }
430: PetscFinalize();
431: return ierr;
432: }
434: /*
435: FormObjFunctionGradient - Evaluates the function and corresponding gradient.
437: Input Parameters:
438: tao - the Tao context
439: P - the input vector
440: ctx - optional aircraft-defined context, as set by TaoSetObjectiveAndGradientRoutine()
442: Output Parameters:
443: f - the newly evaluated function
444: G - the newly evaluated gradient
445: */
446: PetscErrorCode FormObjFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
447: {
448: Aircraft actx = (Aircraft)ctx;
449: TS ts = actx->ts;
450: Vec Q;
451: const PetscScalar *p,*q;
452: PetscScalar *u,*v,*w;
453: PetscInt i;
454: PetscErrorCode ierr;
457: VecGetArrayRead(P,&p);
458: VecGetArray(actx->V,&v);
459: VecGetArray(actx->W,&w);
460: for (i=0; i<actx->nsteps; i++) {
461: v[i] = p[2*i];
462: w[i] = p[2*i+1];
463: }
464: VecRestoreArrayRead(P,&p);
465: VecRestoreArray(actx->V,&v);
466: VecRestoreArray(actx->W,&w);
468: TSSetTime(ts,0.0);
469: TSSetStepNumber(ts,0);
470: TSSetFromOptions(ts);
471: TSSetTimeStep(ts,actx->ftime/actx->nsteps);
473: /* reinitialize system state */
474: VecGetArray(actx->U,&u);
475: u[0] = 2.0;
476: u[1] = 0;
477: VecRestoreArray(actx->U,&u);
479: /* reinitialize the integral value */
480: TSGetCostIntegral(ts,&Q);
481: VecSet(Q,0.0);
483: TSSolve(ts,actx->U);
485: /* Reset initial conditions for the adjoint integration */
486: VecSet(actx->Lambda[0],0.0);
487: VecSet(actx->Mup[0],0.0);
488: TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);
490: TSAdjointSolve(ts);
491: VecCopy(actx->Mup[0],G);
492: TSGetCostIntegral(ts,&Q);
493: VecGetArrayRead(Q,&q);
494: *f = q[0];
495: VecRestoreArrayRead(Q,&q);
496: return(0);
497: }
499: PetscErrorCode FormObjHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx)
500: {
501: Aircraft actx = (Aircraft)ctx;
502: const PetscScalar *p;
503: PetscScalar *harr,*v,*w,one = 1.0;
504: PetscInt ind[1];
505: PetscInt *cols,i;
506: Vec Dir;
507: PetscErrorCode ierr;
510: /* set up control parameters */
511: VecGetArrayRead(P,&p);
512: VecGetArray(actx->V,&v);
513: VecGetArray(actx->W,&w);
514: for (i=0; i<actx->nsteps; i++) {
515: v[i] = p[2*i];
516: w[i] = p[2*i+1];
517: }
518: VecRestoreArrayRead(P,&p);
519: VecRestoreArray(actx->V,&v);
520: VecRestoreArray(actx->W,&w);
522: PetscMalloc1(2*actx->nsteps,&harr);
523: PetscMalloc1(2*actx->nsteps,&cols);
524: for (i=0; i<2*actx->nsteps; i++) cols[i] = i;
525: VecDuplicate(P,&Dir);
526: for (i=0; i<2*actx->nsteps; i++) {
527: ind[0] = i;
528: VecSet(Dir,0.0);
529: VecSetValues(Dir,1,ind,&one,INSERT_VALUES);
530: VecAssemblyBegin(Dir);
531: VecAssemblyEnd(Dir);
532: ComputeObjHessianWithSOA(Dir,harr,actx);
533: MatSetValues(H,1,ind,2*actx->nsteps,cols,harr,INSERT_VALUES);
534: MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY);
535: MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY);
536: if (H != Hpre) {
537: MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY);
538: MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY);
539: }
540: }
541: PetscFree(cols);
542: PetscFree(harr);
543: VecDestroy(&Dir);
544: return(0);
545: }
547: PetscErrorCode MatrixFreeObjHessian(Tao tao, Vec P, Mat H, Mat Hpre, void *ctx)
548: {
549: Aircraft actx = (Aircraft)ctx;
550: PetscScalar *v,*w;
551: const PetscScalar *p;
552: PetscInt i;
553: PetscErrorCode ierr;
556: VecGetArrayRead(P,&p);
557: VecGetArray(actx->V,&v);
558: VecGetArray(actx->W,&w);
559: for (i=0; i<actx->nsteps; i++) {
560: v[i] = p[2*i];
561: w[i] = p[2*i+1];
562: }
563: VecRestoreArrayRead(P,&p);
564: VecRestoreArray(actx->V,&v);
565: VecRestoreArray(actx->W,&w);
566: return(0);
567: }
569: PetscErrorCode MyMatMult(Mat H_shell, Vec X, Vec Y)
570: {
571: PetscScalar *y;
572: void *ptr;
576: MatShellGetContext(H_shell,&ptr);
577: VecGetArray(Y,&y);
578: ComputeObjHessianWithSOA(X,y,(Aircraft)ptr);
579: VecRestoreArray(Y,&y);
580: return(0);
581: }
583: PetscErrorCode ComputeObjHessianWithSOA(Vec Dir,PetscScalar arr[],Aircraft actx)
584: {
585: TS ts = actx->ts;
586: const PetscScalar *z_ptr;
587: PetscScalar *u;
588: Vec Q;
589: PetscInt i;
590: PetscErrorCode ierr;
593: /* Reset TSAdjoint so that AdjointSetUp will be called again */
594: TSAdjointReset(ts);
596: TSSetTime(ts,0.0);
597: TSSetStepNumber(ts,0);
598: TSSetFromOptions(ts);
599: TSSetTimeStep(ts,actx->ftime/actx->nsteps);
600: TSSetCostHessianProducts(actx->ts,1,actx->Lambda2,actx->Mup2,Dir);
602: /* reinitialize system state */
603: VecGetArray(actx->U,&u);
604: u[0] = 2.0;
605: u[1] = 0;
606: VecRestoreArray(actx->U,&u);
608: /* reinitialize the integral value */
609: TSGetCostIntegral(ts,&Q);
610: VecSet(Q,0.0);
612: /* initialize tlm variable */
613: MatZeroEntries(actx->Jacp);
614: MatAssemblyBegin(actx->Jacp,MAT_FINAL_ASSEMBLY);
615: MatAssemblyEnd(actx->Jacp,MAT_FINAL_ASSEMBLY);
616: TSAdjointSetForward(ts,actx->Jacp);
618: TSSolve(ts,actx->U);
620: /* Set terminal conditions for first- and second-order adjonts */
621: VecSet(actx->Lambda[0],0.0);
622: VecSet(actx->Mup[0],0.0);
623: VecSet(actx->Lambda2[0],0.0);
624: VecSet(actx->Mup2[0],0.0);
625: TSSetCostGradients(ts,1,actx->Lambda,actx->Mup);
627: TSGetCostIntegral(ts,&Q);
629: /* Reset initial conditions for the adjoint integration */
630: TSAdjointSolve(ts);
632: /* initial condition does not depend on p, so that lambda is not needed to assemble G */
633: VecGetArrayRead(actx->Mup2[0],&z_ptr);
634: for (i=0; i<2*actx->nsteps; i++) arr[i] = z_ptr[i];
635: VecRestoreArrayRead(actx->Mup2[0],&z_ptr);
637: /* Disable second-order adjoint mode */
638: TSAdjointReset(ts);
639: TSAdjointResetForward(ts);
640: return(0);
641: }
643: /*TEST
644: build:
645: requires: !complex !single
647: test:
648: args: -ts_adapt_type none -ts_type rk -ts_rk_type 3 -viewer_binary_skip_info -tao_monitor -tao_gatol 1e-7
650: test:
651: suffix: 2
652: 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
654: test:
655: suffix: 3
656: 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
657: TEST*/