Actual source code: ex20opt_p.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1: #define c11 1.0
  2: #define c12 0
  3: #define c21 2.0
  4: #define c22 1.0
  5: #define rescale 10

  7: static char help[] = "Solves the van der Pol equation.\n\
  8: Input parameters include:\n";

 10: /*
 11:    Concepts: TS^time-dependent nonlinear problems
 12:    Concepts: TS^van der Pol equation DAE equivalent
 13:    Concepts: Optimization using adjoint sensitivity analysis
 14:    Processors: 1
 15: */
 16: /* ------------------------------------------------------------------------

 18:   Notes:
 19:   This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS.
 20:   The nonlinear problem is written in a DAE equivalent form.
 21:   The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu.
 22:   The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details.
 23:   ------------------------------------------------------------------------- */
 24:  #include <petsctao.h>
 25:  #include <petscts.h>

 27: typedef struct _n_User *User;
 28: struct _n_User {
 29:   PetscReal mu;
 30:   PetscReal next_output;

 32:   /* Sensitivity analysis support */
 33:   PetscReal ftime,x_ob[2];
 34:   Mat       A;                       /* Jacobian matrix */
 35:   Mat       Jacp;                    /* JacobianP matrix */
 36:   Vec       x,lambda[2],mup[2];  /* adjoint variables */
 37: };

 39: PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*);

 41: /*
 42: *  User-defined routines
 43: */
 44: static PetscErrorCode IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
 45: {
 46:   PetscErrorCode    ierr;
 47:   User              user = (User)ctx;
 48:   PetscScalar       *f;
 49:   const PetscScalar *x,*xdot;

 52:   VecGetArrayRead(X,&x);
 53:   VecGetArrayRead(Xdot,&xdot);
 54:   VecGetArray(F,&f);
 55:   f[0] = xdot[0] - x[1];
 56:   f[1] = c21*(xdot[0]-x[1]) + xdot[1] - user->mu*((1.0-x[0]*x[0])*x[1] - x[0]) ;
 57:   VecRestoreArrayRead(X,&x);
 58:   VecRestoreArrayRead(Xdot,&xdot);
 59:   VecRestoreArray(F,&f);
 60:   return(0);
 61: }

 63: static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,void *ctx)
 64: {
 65:   PetscErrorCode    ierr;
 66:   User              user     = (User)ctx;
 67:   PetscInt          rowcol[] = {0,1};
 68:   PetscScalar       J[2][2];
 69:   const PetscScalar *x;
 71:   VecGetArrayRead(X,&x);

 73:   J[0][0] = a;     J[0][1] =  -1.0;
 74:   J[1][0] = c21*a + user->mu*(1.0 + 2.0*x[0]*x[1]);   J[1][1] = -c21 + a - user->mu*(1.0-x[0]*x[0]);

 76:   MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES);
 77:   VecRestoreArrayRead(X,&x);

 79:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 81:   if (A != B) {
 82:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 83:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 84:   }
 85:   return(0);
 86: }

 88: static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,void *ctx)
 89: {
 90:   PetscErrorCode    ierr;
 91:   PetscInt          row[] = {0,1},col[]={0};
 92:   PetscScalar       J[2][1];
 93:   const PetscScalar *x;
 95:   VecGetArrayRead(X,&x);

 97:   J[0][0] = 0;
 98:   J[1][0] = (1.-x[0]*x[0])*x[1]-x[0];
 99:   MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES);
100:   VecRestoreArrayRead(X,&x);

102:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
103:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
104:   return(0);
105: }

107: /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */
108: static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx)
109: {
110:   PetscErrorCode    ierr;
111:   const PetscScalar *x;
112:   PetscReal         tfinal, dt;
113:   User              user = (User)ctx;
114:   Vec               interpolatedX;

117:   TSGetTimeStep(ts,&dt);
118:   TSGetMaxTime(ts,&tfinal);

120:   while (user->next_output <= t && user->next_output <= tfinal) {
121:     VecDuplicate(X,&interpolatedX);
122:     TSInterpolate(ts,user->next_output,interpolatedX);
123:     VecGetArrayRead(interpolatedX,&x);
124:     PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",
125:                        user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),
126:                        (double)PetscRealPart(x[1]));
127:     VecRestoreArrayRead(interpolatedX,&x);
128:     VecDestroy(&interpolatedX);
129:     user->next_output += 0.1;
130:   }
131:   return(0);
132: }

134: int main(int argc,char **argv)
135: {
136:   TS                 ts;            /* nonlinear solver */
137:   Vec                p;
138:   PetscBool          monitor = PETSC_FALSE;
139:   PetscScalar        *x_ptr;
140:   const PetscScalar  *y_ptr;
141:   PetscMPIInt        size;
142:   struct _n_User     user;
143:   PetscErrorCode     ierr;
144:   Tao                tao;
145:   KSP                ksp;
146:   PC                 pc;

148:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149:      Initialize program
150:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
151:   PetscInitialize(&argc,&argv,NULL,help);
152:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
153:   if (size != 1) SETERRQ(PETSC_COMM_SELF,1,"This is a uniprocessor example only!");

155:   /* Create TAO solver and set desired solution method */
156:   TaoCreate(PETSC_COMM_WORLD,&tao);
157:   TaoSetType(tao,TAOBLMVM);

159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160:     Set runtime options
161:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162:   user.next_output = 0.0;
163:   user.mu          = 1.0;
164:   user.ftime       = 1.0;
165:   PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL);
166:   PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL);

168:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169:     Create necessary matrix and vectors, solve same ODE on every process
170:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
171:   MatCreate(PETSC_COMM_WORLD,&user.A);
172:   MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2);
173:   MatSetFromOptions(user.A);
174:   MatSetUp(user.A);
175:   MatCreateVecs(user.A,&user.x,NULL);

177:   MatCreate(PETSC_COMM_WORLD,&user.Jacp);
178:   MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1);
179:   MatSetFromOptions(user.Jacp);
180:   MatSetUp(user.Jacp);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Create timestepping solver context
184:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   TSCreate(PETSC_COMM_WORLD,&ts);
186:   TSSetType(ts,TSCN);
187:   TSSetIFunction(ts,NULL,IFunction,&user);
188:   TSSetIJacobian(ts,user.A,user.A,IJacobian,&user);
189:   TSSetMaxTime(ts,user.ftime);
190:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
191:   if (monitor) {
192:     TSMonitorSet(ts,Monitor,&user,NULL);
193:   }

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196:      Set initial conditions
197:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198:   VecGetArray(user.x,&x_ptr);
199:   x_ptr[0] = 2.0;   x_ptr[1] = -0.66666654321;
200:   VecRestoreArray(user.x,&x_ptr);
201:   TSSetTimeStep(ts,0.03125);

203:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204:      Set runtime options
205:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206:   TSSetFromOptions(ts);

208:   TSSolve(ts,user.x);

210:   VecGetArrayRead(user.x,&y_ptr);
211:   user.x_ob[0] = y_ptr[0];
212:   user.x_ob[1] = y_ptr[1];
213:   VecRestoreArrayRead(user.x,&y_ptr);

215:   /* Create sensitivity variable */
216:   MatCreateVecs(user.A,&user.lambda[0],NULL);
217:   MatCreateVecs(user.Jacp,&user.mup[0],NULL);

219:   /*
220:      Optimization starts
221:   */
222:   /* Set initial solution guess */
223:   MatCreateVecs(user.Jacp,&p,NULL);
224:   VecGetArray(p,&x_ptr);
225:   x_ptr[0] = 1.2;
226:   VecRestoreArray(p,&x_ptr);
227:   TaoSetInitialVector(tao,p);

229:   /* Set routine for function and gradient evaluation */
230:   TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void *)&user);

232:   /* Check for any TAO command line options */
233:   TaoSetFromOptions(tao);
234:   TaoGetKSP(tao,&ksp);
235:   if (ksp) {
236:     KSPGetPC(ksp,&pc);
237:     PCSetType(pc,PCNONE);
238:   }

240:   TaoSolve(tao);

242:   VecView(p,PETSC_VIEWER_STDOUT_WORLD);
243:   /* Free TAO data structures */
244:   TaoDestroy(&tao);

246:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
247:      Free work space.  All PETSc objects should be destroyed when they
248:      are no longer needed.
249:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
250:   MatDestroy(&user.A);
251:   MatDestroy(&user.Jacp);
252:   VecDestroy(&user.x);
253:   VecDestroy(&user.lambda[0]);
254:   VecDestroy(&user.mup[0]);
255:   TSDestroy(&ts);
256:   VecDestroy(&p);
257:   PetscFinalize();
258:   return ierr;
259: }


262: /* ------------------------------------------------------------------ */
263: /*
264:    FormFunctionGradient - Evaluates the function and corresponding gradient.

266:    Input Parameters:
267:    tao - the Tao context
268:    X   - the input vector
269:    ptr - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

271:    Output Parameters:
272:    f   - the newly evaluated function
273:    G   - the newly evaluated gradient
274: */
275: PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
276: {
277:   User              user_ptr = (User)ctx;
278:   TS                ts;
279:   PetscScalar       *x_ptr;
280:   const PetscScalar *y_ptr;
281:   PetscErrorCode    ierr;

284:   VecGetArrayRead(P,&y_ptr);
285:   user_ptr->mu = y_ptr[0];
286:   VecRestoreArrayRead(P,&y_ptr);

288:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289:      Create timestepping solver context
290:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291:   TSCreate(PETSC_COMM_WORLD,&ts);
292:   TSSetType(ts,TSCN);
293:   TSSetIFunction(ts,NULL,IFunction,user_ptr);
294:   TSSetIJacobian(ts,user_ptr->A,user_ptr->A,IJacobian,user_ptr);
295:   TSSetRHSJacobianP(ts,user_ptr->Jacp,RHSJacobianP,user_ptr);

297:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
298:      Set time
299:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
300:   TSSetTime(ts,0.0);
301:   TSSetMaxTime(ts,user_ptr->ftime);
302:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);

304:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
305:     Save trajectory of solution so that TSAdjointSolve() may be used
306:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
307:   TSSetSaveTrajectory(ts);

309:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310:      Set initial conditions
311:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312:   VecGetArray(user_ptr->x,&x_ptr);
313:   x_ptr[0] = 2.0;   x_ptr[1] = -0.66666654321;
314:   VecRestoreArray(user_ptr->x,&x_ptr);

316:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
317:      Set runtime options
318:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
319:   TSSetFromOptions(ts);

321:   TSSolve(ts,user_ptr->x);
322:   VecGetArrayRead(user_ptr->x,&y_ptr);
323:   *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]);
324:   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));
325:   /*   Redet initial conditions for the adjoint integration */
326:   VecGetArray(user_ptr->lambda[0],&x_ptr);
327:   x_ptr[0] = rescale*2.*(y_ptr[0]-user_ptr->x_ob[0]);
328:   x_ptr[1] = rescale*2.*(y_ptr[1]-user_ptr->x_ob[1]);
329:   VecRestoreArrayRead(user_ptr->x,&y_ptr);
330:   VecRestoreArray(user_ptr->lambda[0],&x_ptr);

332:   VecGetArray(user_ptr->mup[0],&x_ptr);
333:   x_ptr[0] = 0.0;
334:   VecRestoreArray(user_ptr->mup[0],&x_ptr);
335:   TSSetCostGradients(ts,1,user_ptr->lambda,user_ptr->mup);

337:   TSAdjointSolve(ts);
338:   VecCopy(user_ptr->mup[0],G);
339:   TSDestroy(&ts);
340:   return(0);
341: }

343: /*TEST
344:     build:
345:       requires: !complex !single
346:     test:
347:       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
348:       output_file: output/ex20opt_p_1.out

350: TEST*/