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