1 static char help[] = "Demonstrates automatic Jacobian generation using ADOL-C for an ODE-constrained optimization problem.\n\ 2 Input parameters include:\n\ 3 -mu : stiffness parameter\n\n"; 4 5 /* 6 Concepts: TS^time-dependent nonlinear problems 7 Concepts: TS^van der Pol equation 8 Concepts: Optimization using adjoint sensitivities 9 Concepts: Automatic differentation using ADOL-C 10 Processors: 1 11 */ 12 /* 13 REQUIRES configuration of PETSc with option --download-adolc. 14 15 For documentation on ADOL-C, see 16 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 17 */ 18 /* ------------------------------------------------------------------------ 19 See ex16opt_ic for a description of the problem being solved. 20 ------------------------------------------------------------------------- */ 21 #include <petsctao.h> 22 #include <petscts.h> 23 #include <petscmat.h> 24 #include "adolc-utils/drivers.cxx" 25 #include <adolc/adolc.h> 26 27 typedef struct _n_User *User; 28 struct _n_User { 29 PetscReal mu; 30 PetscReal next_output; 31 PetscInt steps; 32 33 /* Sensitivity analysis support */ 34 PetscReal ftime,x_ob[2]; 35 Mat A; /* Jacobian matrix */ 36 Vec x,lambda[2]; /* adjoint variables */ 37 38 /* Automatic differentiation support */ 39 AdolcCtx *adctx; 40 }; 41 42 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 43 44 /* 45 'Passive' RHS function, used in residual evaluations during the time integration. 46 */ 47 static PetscErrorCode RHSFunctionPassive(TS ts,PetscReal t,Vec X,Vec F,void *ctx) 48 { 49 User user = (User)ctx; 50 PetscScalar *f; 51 const PetscScalar *x; 52 53 PetscFunctionBeginUser; 54 CHKERRQ(VecGetArrayRead(X,&x)); 55 CHKERRQ(VecGetArray(F,&f)); 56 f[0] = x[1]; 57 f[1] = user->mu*(1.-x[0]*x[0])*x[1]-x[0]; 58 CHKERRQ(VecRestoreArrayRead(X,&x)); 59 CHKERRQ(VecRestoreArray(F,&f)); 60 PetscFunctionReturn(0); 61 } 62 63 /* 64 Trace RHS to mark on tape 1 the dependence of f upon x. This tape is used in generating the 65 Jacobian transform. 66 */ 67 static PetscErrorCode RHSFunctionActive(TS ts,PetscReal t,Vec X,Vec F,void *ctx) 68 { 69 User user = (User)ctx; 70 PetscReal mu = user->mu; 71 PetscScalar *f; 72 const PetscScalar *x; 73 74 adouble f_a[2]; /* adouble for dependent variables */ 75 adouble x_a[2]; /* adouble for independent variables */ 76 77 PetscFunctionBeginUser; 78 CHKERRQ(VecGetArrayRead(X,&x)); 79 CHKERRQ(VecGetArray(F,&f)); 80 81 trace_on(1); /* Start of active section */ 82 x_a[0] <<= x[0]; x_a[1] <<= x[1]; /* Mark as independent */ 83 f_a[0] = x_a[1]; 84 f_a[1] = mu*(1.-x_a[0]*x_a[0])*x_a[1]-x_a[0]; 85 f_a[0] >>= f[0]; f_a[1] >>= f[1]; /* Mark as dependent */ 86 trace_off(1); /* End of active section */ 87 88 CHKERRQ(VecRestoreArrayRead(X,&x)); 89 CHKERRQ(VecRestoreArray(F,&f)); 90 PetscFunctionReturn(0); 91 } 92 93 /* 94 Compute the Jacobian w.r.t. x using PETSc-ADOL-C driver. 95 */ 96 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec X,Mat A,Mat B,void *ctx) 97 { 98 User user=(User)ctx; 99 const PetscScalar *x; 100 101 PetscFunctionBeginUser; 102 CHKERRQ(VecGetArrayRead(X,&x)); 103 CHKERRQ(PetscAdolcComputeRHSJacobian(1,A,x,user->adctx)); 104 CHKERRQ(VecRestoreArrayRead(X,&x)); 105 PetscFunctionReturn(0); 106 } 107 108 /* 109 Monitor timesteps and use interpolation to output at integer multiples of 0.1 110 */ 111 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx) 112 { 113 const PetscScalar *x; 114 PetscReal tfinal, dt, tprev; 115 User user = (User)ctx; 116 117 PetscFunctionBeginUser; 118 CHKERRQ(TSGetTimeStep(ts,&dt)); 119 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 120 CHKERRQ(TSGetPrevTime(ts,&tprev)); 121 CHKERRQ(VecGetArrayRead(X,&x)); 122 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[%.1f] %D TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n",(double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]),(double)PetscRealPart(x[1]))); 123 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"t %.6f (tprev = %.6f) \n",(double)t,(double)tprev)); 124 CHKERRQ(VecGetArrayRead(X,&x)); 125 PetscFunctionReturn(0); 126 } 127 128 int main(int argc,char **argv) 129 { 130 TS ts = NULL; /* nonlinear solver */ 131 Vec ic,r; 132 PetscBool monitor = PETSC_FALSE; 133 PetscScalar *x_ptr; 134 PetscMPIInt size; 135 struct _n_User user; 136 AdolcCtx *adctx; 137 PetscErrorCode ierr; 138 Tao tao; 139 KSP ksp; 140 PC pc; 141 142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 143 Initialize program 144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 145 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 146 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 147 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 148 149 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 150 Set runtime options and create AdolcCtx 151 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 152 CHKERRQ(PetscNew(&adctx)); 153 user.mu = 1.0; 154 user.next_output = 0.0; 155 user.steps = 0; 156 user.ftime = 0.5; 157 adctx->m = 2;adctx->n = 2;adctx->p = 2; 158 user.adctx = adctx; 159 160 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 161 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 162 163 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164 Create necessary matrix and vectors, solve same ODE on every process 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 166 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.A)); 167 CHKERRQ(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 168 CHKERRQ(MatSetFromOptions(user.A)); 169 CHKERRQ(MatSetUp(user.A)); 170 CHKERRQ(MatCreateVecs(user.A,&user.x,NULL)); 171 172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 Set initial conditions 174 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 175 CHKERRQ(VecGetArray(user.x,&x_ptr)); 176 x_ptr[0] = 2.0; x_ptr[1] = 0.66666654321; 177 CHKERRQ(VecRestoreArray(user.x,&x_ptr)); 178 179 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 180 Trace just once on each tape and put zeros on Jacobian diagonal 181 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 182 CHKERRQ(VecDuplicate(user.x,&r)); 183 CHKERRQ(RHSFunctionActive(ts,0.,user.x,r,&user)); 184 CHKERRQ(VecSet(r,0)); 185 CHKERRQ(MatDiagonalSet(user.A,r,INSERT_VALUES)); 186 CHKERRQ(VecDestroy(&r)); 187 188 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 189 Create timestepping solver context 190 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 191 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 192 CHKERRQ(TSSetType(ts,TSRK)); 193 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunctionPassive,&user)); 194 CHKERRQ(TSSetRHSJacobian(ts,user.A,user.A,RHSJacobian,&user)); 195 CHKERRQ(TSSetMaxTime(ts,user.ftime)); 196 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 197 if (monitor) { 198 CHKERRQ(TSMonitorSet(ts,Monitor,&user,NULL)); 199 } 200 201 CHKERRQ(TSSetTime(ts,0.0)); 202 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,user.steps,(double)(user.ftime))); 203 204 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 205 Save trajectory of solution so that TSAdjointSolve() may be used 206 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 207 CHKERRQ(TSSetSaveTrajectory(ts)); 208 209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 210 Set runtime options 211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 212 CHKERRQ(TSSetFromOptions(ts)); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Solve nonlinear system 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 CHKERRQ(TSSolve(ts,user.x)); 218 CHKERRQ(TSGetSolveTime(ts,&(user.ftime))); 219 CHKERRQ(TSGetStepNumber(ts,&user.steps)); 220 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"mu %g, steps %D, ftime %g\n",(double)user.mu,user.steps,(double)user.ftime)); 221 222 CHKERRQ(VecGetArray(user.x,&x_ptr)); 223 user.x_ob[0] = x_ptr[0]; 224 user.x_ob[1] = x_ptr[1]; 225 CHKERRQ(VecRestoreArray(user.x,&x_ptr)); 226 227 CHKERRQ(MatCreateVecs(user.A,&user.lambda[0],NULL)); 228 229 /* Create TAO solver and set desired solution method */ 230 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 231 CHKERRQ(TaoSetType(tao,TAOCG)); 232 233 /* Set initial solution guess */ 234 CHKERRQ(MatCreateVecs(user.A,&ic,NULL)); 235 CHKERRQ(VecGetArray(ic,&x_ptr)); 236 x_ptr[0] = 2.1; 237 x_ptr[1] = 0.7; 238 CHKERRQ(VecRestoreArray(ic,&x_ptr)); 239 240 CHKERRQ(TaoSetSolution(tao,ic)); 241 242 /* Set routine for function and gradient evaluation */ 243 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 244 245 /* Check for any TAO command line options */ 246 CHKERRQ(TaoSetFromOptions(tao)); 247 CHKERRQ(TaoGetKSP(tao,&ksp)); 248 if (ksp) { 249 CHKERRQ(KSPGetPC(ksp,&pc)); 250 CHKERRQ(PCSetType(pc,PCNONE)); 251 } 252 253 CHKERRQ(TaoSetTolerances(tao,1e-10,PETSC_DEFAULT,PETSC_DEFAULT)); 254 255 /* SOLVE THE APPLICATION */ 256 CHKERRQ(TaoSolve(tao)); 257 258 /* Free TAO data structures */ 259 CHKERRQ(TaoDestroy(&tao)); 260 261 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 262 Free work space. All PETSc objects should be destroyed when they 263 are no longer needed. 264 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 265 CHKERRQ(MatDestroy(&user.A)); 266 CHKERRQ(VecDestroy(&user.x)); 267 CHKERRQ(VecDestroy(&user.lambda[0])); 268 CHKERRQ(TSDestroy(&ts)); 269 CHKERRQ(VecDestroy(&ic)); 270 CHKERRQ(PetscFree(adctx)); 271 ierr = PetscFinalize(); 272 return ierr; 273 } 274 275 /* ------------------------------------------------------------------ */ 276 /* 277 FormFunctionGradient - Evaluates the function and corresponding gradient. 278 279 Input Parameters: 280 tao - the Tao context 281 X - the input vector 282 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 283 284 Output Parameters: 285 f - the newly evaluated function 286 G - the newly evaluated gradient 287 */ 288 PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx) 289 { 290 User user = (User)ctx; 291 TS ts; 292 PetscScalar *x_ptr,*y_ptr; 293 294 PetscFunctionBeginUser; 295 CHKERRQ(VecCopy(IC,user->x)); 296 297 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 298 Create timestepping solver context 299 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 300 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 301 CHKERRQ(TSSetType(ts,TSRK)); 302 CHKERRQ(TSSetRHSFunction(ts,NULL,RHSFunctionPassive,user)); 303 /* Set RHS Jacobian for the adjoint integration */ 304 CHKERRQ(TSSetRHSJacobian(ts,user->A,user->A,RHSJacobian,user)); 305 306 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 307 Set time 308 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 309 CHKERRQ(TSSetTime(ts,0.0)); 310 CHKERRQ(TSSetTimeStep(ts,.001)); 311 CHKERRQ(TSSetMaxTime(ts,0.5)); 312 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP)); 313 314 CHKERRQ(TSSetTolerances(ts,1e-7,NULL,1e-7,NULL)); 315 316 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 317 Save trajectory of solution so that TSAdjointSolve() may be used 318 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 319 CHKERRQ(TSSetSaveTrajectory(ts)); 320 321 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 322 Set runtime options 323 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 324 CHKERRQ(TSSetFromOptions(ts)); 325 326 CHKERRQ(TSSolve(ts,user->x)); 327 CHKERRQ(TSGetSolveTime(ts,&user->ftime)); 328 CHKERRQ(TSGetStepNumber(ts,&user->steps)); 329 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"mu %.6f, steps %D, ftime %g\n",(double)user->mu,user->steps,(double)user->ftime)); 330 331 CHKERRQ(VecGetArray(user->x,&x_ptr)); 332 *f = (x_ptr[0]-user->x_ob[0])*(x_ptr[0]-user->x_ob[0])+(x_ptr[1]-user->x_ob[1])*(x_ptr[1]-user->x_ob[1]); 333 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"Observed value y_ob=[%f; %f], ODE solution y=[%f;%f], Cost function f=%f\n",(double)user->x_ob[0],(double)user->x_ob[1],(double)x_ptr[0],(double)x_ptr[1],(double)(*f))); 334 335 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 336 Adjoint model starts here 337 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 338 /* Redet initial conditions for the adjoint integration */ 339 CHKERRQ(VecGetArray(user->lambda[0],&y_ptr)); 340 y_ptr[0] = 2.*(x_ptr[0]-user->x_ob[0]); 341 y_ptr[1] = 2.*(x_ptr[1]-user->x_ob[1]); 342 CHKERRQ(VecRestoreArray(user->lambda[0],&y_ptr)); 343 CHKERRQ(VecRestoreArray(user->x,&x_ptr)); 344 CHKERRQ(TSSetCostGradients(ts,1,user->lambda,NULL)); 345 346 CHKERRQ(TSAdjointSolve(ts)); 347 348 CHKERRQ(VecCopy(user->lambda[0],G)); 349 350 CHKERRQ(TSDestroy(&ts)); 351 PetscFunctionReturn(0); 352 } 353 354 /*TEST 355 356 build: 357 requires: double !complex adolc 358 359 test: 360 suffix: 1 361 args: -ts_rhs_jacobian_test_mult_transpose FALSE -tao_max_it 2 -ts_rhs_jacobian_test_mult FALSE 362 output_file: output/ex16opt_ic_1.out 363 364 TEST*/ 365