Actual source code: ex20opt_p.c
petsc-3.11.4 2019-09-28
1: #define c21 2.0
2: #define rescale 10
4: static char help[] = "Solves the van der Pol equation.\n\
5: Input parameters include:\n";
7: /*
8: Concepts: TS^time-dependent nonlinear problems
9: Concepts: TS^van der Pol equation DAE equivalent
10: Concepts: Optimization using adjoint sensitivity analysis
11: Processors: 1
12: */
13: /* ------------------------------------------------------------------------
15: Notes:
16: This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
17: The nonlinear problem is written in a DAE equivalent form.
18: The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
19: The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
20: ------------------------------------------------------------------------- */
21: #include <petsctao.h>
22: #include <petscts.h>
24: typedef struct _n_User *User;
25: struct _n_User {
26: PetscReal mu;
27: PetscReal next_output;
29: /* Sensitivity analysis support */
30: PetscReal ftime,x_ob[2];
31: Mat A; /* Jacobian matrix */
32: Mat Jacp; /* JacobianP matrix */
33: Vec x,lambda[2],mup[2]; /* adjoint variables */
34: };
36: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);
38: /*
39: * User-defined routines
40: */
41: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
42: {
43: PetscErrorCode ierr;
44: User user = (User)ctx;
45: PetscScalar *f;
46: const PetscScalar *x,*xdot;
49: VecGetArrayRead(X,&x);
50: VecGetArrayRead(Xdot,&xdot);
51: VecGetArray(F,&f);
52: f[0] = xdot[0] - x[1];
53: f[1] = c21*(xdot[0]-x[1]) + xdot[1] - user->mu*((1.0-x[0]*x[0])*x[1] - x[0]) ;
54: VecRestoreArrayRead(X,&x);
55: VecRestoreArrayRead(Xdot,&xdot);
56: VecRestoreArray(F,&f);
57: return(0);
58: }
60: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
61: {
62: PetscErrorCode ierr;
63: User user = (User)ctx;
64: PetscInt rowcol[] = {0,1};
65: PetscScalar J[2][2];
66: const PetscScalar *x;
69: VecGetArrayRead(X,&x);
70: J[0][0] = a;
71: J[0][1] = -1.0;
72: J[1][0] = c21*a + user->mu*(1.0 + 2.0*x[0]*x[1]);
73: J[1][1] = -c21 + a - user->mu*(1.0-x[0]*x[0]);
74: MatZeroEntries(A);
75: MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
76: VecRestoreArrayRead(X,&x);
77: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
78: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
79: if (B && A != B) {
80: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
81: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
82: }
83: return(0);
84: }
86: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx)
87: {
88: PetscErrorCode ierr;
89: PetscInt row[] = {0,1},col[]={0};
90: PetscScalar J[2][1];
91: const PetscScalar *x;
94: MatZeroEntries(A);
95: VecGetArrayRead(X,&x);
96: J[0][0] = 0;
97: J[1][0] = (1.-x[0]*x[0])*x[1]-x[0];
98: MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
99: VecRestoreArrayRead(X,&x);
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
102: return(0);
103: }
105: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
106: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
107: {
108: PetscErrorCode ierr;
109: const PetscScalar *x;
110: PetscReal tfinal, dt;
111: User user = (User)ctx;
112: Vec interpolatedX;
115: TSGetTimeStep(ts,&dt);
116: TSGetMaxTime(ts,&tfinal);
118: while (user->next_output <= t && user->next_output <= tfinal) {
119: VecDuplicate(X,&interpolatedX);
120: TSInterpolate(ts,user->next_output,interpolatedX);
121: VecGetArrayRead(interpolatedX,&x);
122: PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",
123: user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
124: (double)PetscRealPart(x[1]));
125: VecRestoreArrayRead(interpolatedX,&x);
126: VecDestroy(&interpolatedX);
127: user->next_output += 0.1;
128: }
129: return(0);
130: }
132: int main(int argc,char **argv)
133: {
134: TS ts; /* nonlinear solver */
135: Vec p;
136: PetscBool monitor = PETSC_FALSE;
137: PetscScalar *x_ptr;
138: const PetscScalar *y_ptr;
139: PetscMPIInt size;
140: struct _n_User user;
141: PetscErrorCode ierr;
142: Tao tao;
143: KSP ksp;
144: PC pc;
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Initialize program
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
150: MPI_Comm_size(PETSC_COMM_WORLD,&size);
151: if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");
153: /* Create TAO solver and set desired solution method */
154: TaoCreate(PETSC_COMM_WORLD,&tao);
155: TaoSetType(tao,TAOBLMVM);
157: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158: Set runtime options
159: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160: user.next_output = 0.0;
161: user.mu = 1.0;
162: user.ftime = 1.0;
163: PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
164: PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);
166: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: Create necessary matrix and vectors, solve same ODE on every process
168: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: MatCreate(PETSC_COMM_WORLD,&user.A);
170: MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
171: MatSetFromOptions(user.A);
172: MatSetUp(user.A);
173: MatCreateVecs(user.A,&user.x,NULL);
175: MatCreate(PETSC_COMM_WORLD,&user.Jacp);
176: MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
177: MatSetFromOptions(user.Jacp);
178: MatSetUp(user.Jacp);
180: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181: Create timestepping solver context
182: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
183: TSCreate(PETSC_COMM_WORLD,&ts);
184: TSSetType(ts,TSCN);
185: TSSetIFunction(ts,NULL,IFunction,&user);
186: TSSetIJacobian(ts,user.A,user.A,IJacobian,&user);
187: TSSetMaxTime(ts,user.ftime);
188: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
189: if (monitor) {
190: TSMonitorSet(ts,Monitor,&user,NULL);
191: }
193: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: Set initial conditions
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: VecGetArray(user.x,&x_ptr);
197: x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321;
198: VecRestoreArray(user.x,&x_ptr);
199: TSSetTimeStep(ts,0.03125);
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Set runtime options
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: TSSetFromOptions(ts);
206: TSSolve(ts,user.x);
208: VecGetArrayRead(user.x,&y_ptr);
209: user.x_ob[0] = y_ptr[0];
210: user.x_ob[1] = y_ptr[1];
211: VecRestoreArrayRead(user.x,&y_ptr);
213: /* Create sensitivity variable */
214: MatCreateVecs(user.A,&user.lambda[0],NULL);
215: MatCreateVecs(user.Jacp,&user.mup[0],NULL);
217: /*
218: Optimization starts
219: */
220: /* Set initial solution guess */
221: MatCreateVecs(user.Jacp,&p,NULL);
222: VecGetArray(p,&x_ptr);
223: x_ptr[0] = 1.2;
224: VecRestoreArray(p,&x_ptr);
225: TaoSetInitialVector(tao,p);
227: /* Set routine for function and gradient evaluation */
228: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);
230: /* Check for any TAO command line options */
231: TaoSetFromOptions(tao);
232: TaoGetKSP(tao,&ksp);
233: if (ksp) {
234: KSPGetPC(ksp,&pc);
235: PCSetType(pc,PCNONE);
236: }
238: TaoSolve(tao);
240: VecView(p,PETSC_VIEWER_STDOUT_WORLD);
241: /* Free TAO data structures */
242: TaoDestroy(&tao);
244: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
245: Free work space. All PETSc objects should be destroyed when they
246: are no longer needed.
247: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
248: MatDestroy(&user.A);
249: MatDestroy(&user.Jacp);
250: VecDestroy(&user.x);
251: VecDestroy(&user.lambda[0]);
252: VecDestroy(&user.mup[0]);
253: TSDestroy(&ts);
254: VecDestroy(&p);
255: PetscFinalize();
256: return ierr;
257: }
260: /* ------------------------------------------------------------------ */
261: /*
262: FormFunctionGradient - Evaluates the function and corresponding gradient.
264: Input Parameters:
265: tao - the Tao context
266: X - the input vector
267: ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
269: Output Parameters:
270: f - the newly evaluated function
271: G - the newly evaluated gradient
272: */
273: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
274: {
275: User user_ptr = (User)ctx;
276: TS ts;
277: PetscScalar *x_ptr;
278: const PetscScalar *y_ptr;
279: PetscErrorCode ierr;
282: VecGetArrayRead(P,&y_ptr);
283: user_ptr->mu = y_ptr[0];
284: VecRestoreArrayRead(P,&y_ptr);
286: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
287: Create timestepping solver context
288: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
289: TSCreate(PETSC_COMM_WORLD,&ts);
290: TSSetType(ts,TSCN);
291: TSSetIFunction(ts,NULL,IFunction,user_ptr);
292: TSSetIJacobian(ts,user_ptr->A,user_ptr->A,IJacobian,user_ptr);
293: TSSetRHSJacobianP(ts,user_ptr->Jacp,RHSJacobianP,user_ptr);
295: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
296: Set time
297: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
298: TSSetTime(ts,0.0);
299: TSSetMaxTime(ts,user_ptr->ftime);
300: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
302: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
303: Save trajectory of solution so that TSAdjointSolve() may be used
304: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
305: TSSetSaveTrajectory(ts);
307: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308: Set initial conditions
309: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
310: VecGetArray(user_ptr->x,&x_ptr);
311: x_ptr[0] = 2.0; x_ptr[1] = -0.66666654321;
312: VecRestoreArray(user_ptr->x,&x_ptr);
314: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
315: Set runtime options
316: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
317: TSSetFromOptions(ts);
319: TSSolve(ts,user_ptr->x);
320: VecGetArrayRead(user_ptr->x,&y_ptr);
321: *f = rescale*(y_ptr[0]-user_ptr->x_ob[0])*(y_ptr[0]-user_ptr->x_ob[0])+rescale*(y_ptr[1]-user_ptr->x_ob[1])*(y_ptr[1]-user_ptr->x_ob[1]);
322: PetscPrintf(PETSC_COMM_WORLD,"Observed value y_ob=[%g; %g], ODE solution y=[%g;%g], Cost function f=%g\n",(double)user_ptr->x_ob[0],(double)user_ptr->x_ob[1],(double)y_ptr[0],(double)y_ptr[1],(double)(*f));
323: /* Redet initial conditions for the adjoint integration */
324: VecGetArray(user_ptr->lambda[0],&x_ptr);
325: x_ptr[0] = rescale*2.*(y_ptr[0]-user_ptr->x_ob[0]);
326: x_ptr[1] = rescale*2.*(y_ptr[1]-user_ptr->x_ob[1]);
327: VecRestoreArrayRead(user_ptr->x,&y_ptr);
328: VecRestoreArray(user_ptr->lambda[0],&x_ptr);
330: VecGetArray(user_ptr->mup[0],&x_ptr);
331: x_ptr[0] = 0.0;
332: VecRestoreArray(user_ptr->mup[0],&x_ptr);
333: TSSetCostGradients(ts,1,user_ptr->lambda,user_ptr->mup);
335: TSAdjointSolve(ts);
336: VecCopy(user_ptr->mup[0],G);
337: TSDestroy(&ts);
338: return(0);
339: }
341: /*TEST
342: build:
343: requires: !complex !single
344: test:
345: args: -monitor 0 -ts_type theta -ts_theta_endpoint -ts_theta_theta 0.5 -viewer_binary_skip_info -tao_view -ts_trajectory_dirname ex20opt_pdir -tao_blmvm_mat_lmvm_scale_type none
346: output_file: output/ex20opt_p_1.out
348: TEST*/