1 static char help[] = "Solves a ODE-constrained optimization problem -- finding the optimal initial conditions for the van der Pol equation.\n"; 2 3 /* 4 Concepts: TS^time-dependent nonlinear problems 5 Concepts: TS^van der Pol equation DAE equivalent 6 Concepts: TS^Optimization using adjoint sensitivity analysis 7 Processors: 1 8 */ 9 /* 10 Notes: 11 This code demonstrates how to solve an ODE-constrained optimization problem with TAO, TSAdjoint and TS. 12 The nonlinear problem is written in an ODE equivalent form. 13 The objective is to minimize the difference between observation and model prediction by finding optimal values for initial conditions. 14 The gradient is computed with the discrete adjoint of an implicit method or an explicit method, see ex20adj.c for details. 15 */ 16 17 #include <petsctao.h> 18 #include <petscts.h> 19 20 typedef struct _n_User *User; 21 struct _n_User { 22 TS ts; 23 PetscReal mu; 24 PetscReal next_output; 25 26 /* Sensitivity analysis support */ 27 PetscInt steps; 28 PetscReal ftime; 29 Mat A; /* Jacobian matrix for ODE */ 30 Mat Jacp; /* JacobianP matrix for ODE*/ 31 Mat H; /* Hessian matrix for optimization */ 32 Vec U,Lambda[1],Mup[1]; /* first-order adjoint variables */ 33 Vec Lambda2[2]; /* second-order adjoint variables */ 34 Vec Ihp1[1]; /* working space for Hessian evaluations */ 35 Vec Dir; /* direction vector */ 36 PetscReal ob[2]; /* observation used by the cost function */ 37 PetscBool implicitform; /* implicit ODE? */ 38 }; 39 PetscErrorCode Adjoint2(Vec,PetscScalar[],User); 40 41 /* ----------------------- Explicit form of the ODE -------------------- */ 42 43 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 44 { 45 User user = (User)ctx; 46 PetscScalar *f; 47 const PetscScalar *u; 48 49 PetscFunctionBeginUser; 50 CHKERRQ(VecGetArrayRead(U,&u)); 51 CHKERRQ(VecGetArray(F,&f)); 52 f[0] = u[1]; 53 f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]); 54 CHKERRQ(VecRestoreArrayRead(U,&u)); 55 CHKERRQ(VecRestoreArray(F,&f)); 56 PetscFunctionReturn(0); 57 } 58 59 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 60 { 61 User user = (User)ctx; 62 PetscReal mu = user->mu; 63 PetscInt rowcol[] = {0,1}; 64 PetscScalar J[2][2]; 65 const PetscScalar *u; 66 67 PetscFunctionBeginUser; 68 CHKERRQ(VecGetArrayRead(U,&u)); 69 J[0][0] = 0; 70 J[1][0] = -mu*(2.0*u[1]*u[0]+1.); 71 J[0][1] = 1.0; 72 J[1][1] = mu*(1.0-u[0]*u[0]); 73 CHKERRQ(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 74 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 75 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 76 if (A != B) { 77 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 78 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 79 } 80 CHKERRQ(VecRestoreArrayRead(U,&u)); 81 PetscFunctionReturn(0); 82 } 83 84 static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 85 { 86 const PetscScalar *vl,*vr,*u; 87 PetscScalar *vhv; 88 PetscScalar dJdU[2][2][2]={{{0}}}; 89 PetscInt i,j,k; 90 User user = (User)ctx; 91 92 PetscFunctionBeginUser; 93 CHKERRQ(VecGetArrayRead(U,&u)); 94 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 95 CHKERRQ(VecGetArrayRead(Vr,&vr)); 96 CHKERRQ(VecGetArray(VHV[0],&vhv)); 97 98 dJdU[1][0][0] = -2.*user->mu*u[1]; 99 dJdU[1][1][0] = -2.*user->mu*u[0]; 100 dJdU[1][0][1] = -2.*user->mu*u[0]; 101 for (j=0;j<2;j++) { 102 vhv[j] = 0; 103 for (k=0;k<2;k++) 104 for (i=0;i<2;i++) 105 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 106 } 107 CHKERRQ(VecRestoreArrayRead(U,&u)); 108 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 109 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 110 CHKERRQ(VecRestoreArray(VHV[0],&vhv)); 111 PetscFunctionReturn(0); 112 } 113 114 /* ----------------------- Implicit form of the ODE -------------------- */ 115 116 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 117 { 118 User user = (User)ctx; 119 const PetscScalar *u,*udot; 120 PetscScalar *f; 121 122 PetscFunctionBeginUser; 123 CHKERRQ(VecGetArrayRead(U,&u)); 124 CHKERRQ(VecGetArrayRead(Udot,&udot)); 125 CHKERRQ(VecGetArray(F,&f)); 126 f[0] = udot[0] - u[1]; 127 f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ; 128 CHKERRQ(VecRestoreArrayRead(U,&u)); 129 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 130 CHKERRQ(VecRestoreArray(F,&f)); 131 PetscFunctionReturn(0); 132 } 133 134 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx) 135 { 136 User user = (User)ctx; 137 PetscInt rowcol[] = {0,1}; 138 PetscScalar J[2][2]; 139 const PetscScalar *u; 140 141 PetscFunctionBeginUser; 142 CHKERRQ(VecGetArrayRead(U,&u)); 143 J[0][0] = a; J[0][1] = -1.0; 144 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]); 145 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 146 CHKERRQ(VecRestoreArrayRead(U,&u)); 147 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 148 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 149 if (A != B) { 150 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 151 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 152 } 153 PetscFunctionReturn(0); 154 } 155 156 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 157 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec U,void *ctx) 158 { 159 PetscErrorCode ierr; 160 const PetscScalar *u; 161 PetscReal tfinal, dt; 162 User user = (User)ctx; 163 Vec interpolatedU; 164 165 PetscFunctionBeginUser; 166 CHKERRQ(TSGetTimeStep(ts,&dt)); 167 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 168 169 while (user->next_output <= t && user->next_output <= tfinal) { 170 CHKERRQ(VecDuplicate(U,&interpolatedU)); 171 CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedU)); 172 CHKERRQ(VecGetArrayRead(interpolatedU,&u)); 173 ierr = PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n", 174 (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]), 175 (double)PetscRealPart(u[1]));CHKERRQ(ierr); 176 CHKERRQ(VecRestoreArrayRead(interpolatedU,&u)); 177 CHKERRQ(VecDestroy(&interpolatedU)); 178 user->next_output += 0.1; 179 } 180 PetscFunctionReturn(0); 181 } 182 183 static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 184 { 185 const PetscScalar *vl,*vr,*u; 186 PetscScalar *vhv; 187 PetscScalar dJdU[2][2][2]={{{0}}}; 188 PetscInt i,j,k; 189 User user = (User)ctx; 190 191 PetscFunctionBeginUser; 192 CHKERRQ(VecGetArrayRead(U,&u)); 193 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 194 CHKERRQ(VecGetArrayRead(Vr,&vr)); 195 CHKERRQ(VecGetArray(VHV[0],&vhv)); 196 dJdU[1][0][0] = 2.*user->mu*u[1]; 197 dJdU[1][1][0] = 2.*user->mu*u[0]; 198 dJdU[1][0][1] = 2.*user->mu*u[0]; 199 for (j=0;j<2;j++) { 200 vhv[j] = 0; 201 for (k=0;k<2;k++) 202 for (i=0;i<2;i++) 203 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 204 } 205 CHKERRQ(VecRestoreArrayRead(U,&u)); 206 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 207 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 208 CHKERRQ(VecRestoreArray(VHV[0],&vhv)); 209 PetscFunctionReturn(0); 210 } 211 212 /* ------------------ User-defined routines for TAO -------------------------- */ 213 214 static PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx) 215 { 216 User user_ptr = (User)ctx; 217 TS ts = user_ptr->ts; 218 const PetscScalar *x_ptr; 219 PetscScalar *y_ptr; 220 221 PetscFunctionBeginUser; 222 CHKERRQ(VecCopy(IC,user_ptr->U)); /* set up the initial condition */ 223 224 CHKERRQ(TSSetTime(ts,0.0)); 225 CHKERRQ(TSSetStepNumber(ts,0)); 226 CHKERRQ(TSSetTimeStep(ts,0.001)); /* can be overwritten by command line options */ 227 CHKERRQ(TSSetFromOptions(ts)); 228 CHKERRQ(TSSolve(ts,user_ptr->U)); 229 230 CHKERRQ(VecGetArrayRead(user_ptr->U,&x_ptr)); 231 CHKERRQ(VecGetArray(user_ptr->Lambda[0],&y_ptr)); 232 *f = (x_ptr[0]-user_ptr->ob[0])*(x_ptr[0]-user_ptr->ob[0])+(x_ptr[1]-user_ptr->ob[1])*(x_ptr[1]-user_ptr->ob[1]); 233 y_ptr[0] = 2.*(x_ptr[0]-user_ptr->ob[0]); 234 y_ptr[1] = 2.*(x_ptr[1]-user_ptr->ob[1]); 235 CHKERRQ(VecRestoreArray(user_ptr->Lambda[0],&y_ptr)); 236 CHKERRQ(VecRestoreArrayRead(user_ptr->U,&x_ptr)); 237 238 CHKERRQ(TSSetCostGradients(ts,1,user_ptr->Lambda,NULL)); 239 CHKERRQ(TSAdjointSolve(ts)); 240 CHKERRQ(VecCopy(user_ptr->Lambda[0],G)); 241 PetscFunctionReturn(0); 242 } 243 244 static PetscErrorCode FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx) 245 { 246 User user_ptr = (User)ctx; 247 PetscScalar harr[2]; 248 PetscScalar *x_ptr; 249 const PetscInt rows[2] = {0,1}; 250 PetscInt col; 251 252 PetscFunctionBeginUser; 253 CHKERRQ(VecCopy(U,user_ptr->U)); 254 CHKERRQ(VecGetArray(user_ptr->Dir,&x_ptr)); 255 x_ptr[0] = 1.; 256 x_ptr[1] = 0.; 257 CHKERRQ(VecRestoreArray(user_ptr->Dir,&x_ptr)); 258 CHKERRQ(Adjoint2(user_ptr->U,harr,user_ptr)); 259 col = 0; 260 CHKERRQ(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 261 262 CHKERRQ(VecCopy(U,user_ptr->U)); 263 CHKERRQ(VecGetArray(user_ptr->Dir,&x_ptr)); 264 x_ptr[0] = 0.; 265 x_ptr[1] = 1.; 266 CHKERRQ(VecRestoreArray(user_ptr->Dir,&x_ptr)); 267 CHKERRQ(Adjoint2(user_ptr->U,harr,user_ptr)); 268 col = 1; 269 CHKERRQ(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 270 271 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 272 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 273 if (H != Hpre) { 274 CHKERRQ(MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY)); 275 CHKERRQ(MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY)); 276 } 277 PetscFunctionReturn(0); 278 } 279 280 static PetscErrorCode MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx) 281 { 282 User user_ptr = (User)ctx; 283 284 PetscFunctionBeginUser; 285 CHKERRQ(VecCopy(U,user_ptr->U)); 286 PetscFunctionReturn(0); 287 } 288 289 /* ------------ Routines calculating second-order derivatives -------------- */ 290 291 /* 292 Compute the Hessian-vector product for the cost function using Second-order adjoint 293 */ 294 PetscErrorCode Adjoint2(Vec U,PetscScalar arr[],User ctx) 295 { 296 TS ts = ctx->ts; 297 PetscScalar *x_ptr,*y_ptr; 298 Mat tlmsen; 299 300 PetscFunctionBeginUser; 301 CHKERRQ(TSAdjointReset(ts)); 302 303 CHKERRQ(TSSetTime(ts,0.0)); 304 CHKERRQ(TSSetStepNumber(ts,0)); 305 CHKERRQ(TSSetTimeStep(ts,0.001)); 306 CHKERRQ(TSSetFromOptions(ts)); 307 CHKERRQ(TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir)); 308 CHKERRQ(TSAdjointSetForward(ts,NULL)); 309 CHKERRQ(TSSolve(ts,U)); 310 311 /* Set terminal conditions for first- and second-order adjonts */ 312 CHKERRQ(VecGetArray(U,&x_ptr)); 313 CHKERRQ(VecGetArray(ctx->Lambda[0],&y_ptr)); 314 y_ptr[0] = 2.*(x_ptr[0]-ctx->ob[0]); 315 y_ptr[1] = 2.*(x_ptr[1]-ctx->ob[1]); 316 CHKERRQ(VecRestoreArray(ctx->Lambda[0],&y_ptr)); 317 CHKERRQ(VecRestoreArray(U,&x_ptr)); 318 CHKERRQ(TSForwardGetSensitivities(ts,NULL,&tlmsen)); 319 CHKERRQ(MatDenseGetColumn(tlmsen,0,&x_ptr)); 320 CHKERRQ(VecGetArray(ctx->Lambda2[0],&y_ptr)); 321 y_ptr[0] = 2.*x_ptr[0]; 322 y_ptr[1] = 2.*x_ptr[1]; 323 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 324 CHKERRQ(MatDenseRestoreColumn(tlmsen,&x_ptr)); 325 326 CHKERRQ(TSSetCostGradients(ts,1,ctx->Lambda,NULL)); 327 if (ctx->implicitform) { 328 CHKERRQ(TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 329 } else { 330 CHKERRQ(TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 331 } 332 CHKERRQ(TSAdjointSolve(ts)); 333 334 CHKERRQ(VecGetArray(ctx->Lambda2[0],&x_ptr)); 335 arr[0] = x_ptr[0]; 336 arr[1] = x_ptr[1]; 337 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&x_ptr)); 338 339 CHKERRQ(TSAdjointReset(ts)); 340 CHKERRQ(TSAdjointResetForward(ts)); 341 PetscFunctionReturn(0); 342 } 343 344 PetscErrorCode FiniteDiff(Vec U,PetscScalar arr[],User ctx) 345 { 346 Vec Up,G,Gp; 347 const PetscScalar eps = PetscRealConstant(1e-7); 348 PetscScalar *u; 349 Tao tao = NULL; 350 PetscReal f; 351 352 PetscFunctionBeginUser; 353 CHKERRQ(VecDuplicate(U,&Up)); 354 CHKERRQ(VecDuplicate(U,&G)); 355 CHKERRQ(VecDuplicate(U,&Gp)); 356 357 CHKERRQ(FormFunctionGradient(tao,U,&f,G,ctx)); 358 359 CHKERRQ(VecCopy(U,Up)); 360 CHKERRQ(VecGetArray(Up,&u)); 361 u[0] += eps; 362 CHKERRQ(VecRestoreArray(Up,&u)); 363 CHKERRQ(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 364 CHKERRQ(VecAXPY(Gp,-1,G)); 365 CHKERRQ(VecScale(Gp,1./eps)); 366 CHKERRQ(VecGetArray(Gp,&u)); 367 arr[0] = u[0]; 368 arr[1] = u[1]; 369 CHKERRQ(VecRestoreArray(Gp,&u)); 370 371 CHKERRQ(VecCopy(U,Up)); 372 CHKERRQ(VecGetArray(Up,&u)); 373 u[1] += eps; 374 CHKERRQ(VecRestoreArray(Up,&u)); 375 CHKERRQ(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 376 CHKERRQ(VecAXPY(Gp,-1,G)); 377 CHKERRQ(VecScale(Gp,1./eps)); 378 CHKERRQ(VecGetArray(Gp,&u)); 379 arr[2] = u[0]; 380 arr[3] = u[1]; 381 CHKERRQ(VecRestoreArray(Gp,&u)); 382 383 CHKERRQ(VecDestroy(&G)); 384 CHKERRQ(VecDestroy(&Gp)); 385 CHKERRQ(VecDestroy(&Up)); 386 PetscFunctionReturn(0); 387 } 388 389 static PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y) 390 { 391 User user_ptr; 392 PetscScalar *y_ptr; 393 394 PetscFunctionBeginUser; 395 CHKERRQ(MatShellGetContext(mat,&user_ptr)); 396 CHKERRQ(VecCopy(svec,user_ptr->Dir)); 397 CHKERRQ(VecGetArray(y,&y_ptr)); 398 CHKERRQ(Adjoint2(user_ptr->U,y_ptr,user_ptr)); 399 CHKERRQ(VecRestoreArray(y,&y_ptr)); 400 PetscFunctionReturn(0); 401 } 402 403 int main(int argc,char **argv) 404 { 405 PetscBool monitor = PETSC_FALSE,mf = PETSC_TRUE; 406 PetscInt mode = 0; 407 PetscMPIInt size; 408 struct _n_User user; 409 Vec x; /* working vector for TAO */ 410 PetscScalar *x_ptr,arr[4]; 411 PetscScalar ic1 = 2.2,ic2 = -0.7; /* initial guess for TAO */ 412 Tao tao; 413 KSP ksp; 414 PC pc; 415 PetscErrorCode ierr; 416 417 /* Initialize program */ 418 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 419 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 420 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 421 422 /* Set runtime options */ 423 user.next_output = 0.0; 424 user.mu = 1.0e3; 425 user.steps = 0; 426 user.ftime = 0.5; 427 user.implicitform = PETSC_TRUE; 428 429 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 430 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 431 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL)); 432 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL)); 433 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL)); 434 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL)); /* matrix-free hessian for optimization */ 435 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL)); 436 437 /* Create necessary matrix and vectors, solve same ODE on every process */ 438 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.A)); 439 CHKERRQ(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 440 CHKERRQ(MatSetFromOptions(user.A)); 441 CHKERRQ(MatSetUp(user.A)); 442 CHKERRQ(MatCreateVecs(user.A,&user.U,NULL)); 443 CHKERRQ(MatCreateVecs(user.A,&user.Dir,NULL)); 444 CHKERRQ(MatCreateVecs(user.A,&user.Lambda[0],NULL)); 445 CHKERRQ(MatCreateVecs(user.A,&user.Lambda2[0],NULL)); 446 CHKERRQ(MatCreateVecs(user.A,&user.Ihp1[0],NULL)); 447 448 /* Create timestepping solver context */ 449 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&user.ts)); 450 CHKERRQ(TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 451 if (user.implicitform) { 452 CHKERRQ(TSSetIFunction(user.ts,NULL,IFunction,&user)); 453 CHKERRQ(TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user)); 454 CHKERRQ(TSSetType(user.ts,TSCN)); 455 } else { 456 CHKERRQ(TSSetRHSFunction(user.ts,NULL,RHSFunction,&user)); 457 CHKERRQ(TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user)); 458 CHKERRQ(TSSetType(user.ts,TSRK)); 459 } 460 CHKERRQ(TSSetMaxTime(user.ts,user.ftime)); 461 CHKERRQ(TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP)); 462 463 if (monitor) { 464 CHKERRQ(TSMonitorSet(user.ts,Monitor,&user,NULL)); 465 } 466 467 /* Set ODE initial conditions */ 468 CHKERRQ(VecGetArray(user.U,&x_ptr)); 469 x_ptr[0] = 2.0; 470 x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu); 471 CHKERRQ(VecRestoreArray(user.U,&x_ptr)); 472 473 /* Set runtime options */ 474 CHKERRQ(TSSetFromOptions(user.ts)); 475 476 /* Obtain the observation */ 477 CHKERRQ(TSSolve(user.ts,user.U)); 478 CHKERRQ(VecGetArray(user.U,&x_ptr)); 479 user.ob[0] = x_ptr[0]; 480 user.ob[1] = x_ptr[1]; 481 CHKERRQ(VecRestoreArray(user.U,&x_ptr)); 482 483 CHKERRQ(VecDuplicate(user.U,&x)); 484 CHKERRQ(VecGetArray(x,&x_ptr)); 485 x_ptr[0] = ic1; 486 x_ptr[1] = ic2; 487 CHKERRQ(VecRestoreArray(x,&x_ptr)); 488 489 /* Save trajectory of solution so that TSAdjointSolve() may be used */ 490 CHKERRQ(TSSetSaveTrajectory(user.ts)); 491 492 /* Compare finite difference and second-order adjoint. */ 493 switch (mode) { 494 case 2 : 495 CHKERRQ(FiniteDiff(x,arr,&user)); 496 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n")); 497 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g %g\n%g %g\n",(double)arr[0],(double)arr[1],(double)arr[2],(double)arr[3])); 498 break; 499 case 1 : /* Compute the Hessian column by column */ 500 CHKERRQ(VecCopy(x,user.U)); 501 CHKERRQ(VecGetArray(user.Dir,&x_ptr)); 502 x_ptr[0] = 1.; 503 x_ptr[1] = 0.; 504 CHKERRQ(VecRestoreArray(user.Dir,&x_ptr)); 505 CHKERRQ(Adjoint2(user.U,arr,&user)); 506 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n")); 507 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1])); 508 CHKERRQ(VecCopy(x,user.U)); 509 CHKERRQ(VecGetArray(user.Dir,&x_ptr)); 510 x_ptr[0] = 0.; 511 x_ptr[1] = 1.; 512 CHKERRQ(VecRestoreArray(user.Dir,&x_ptr)); 513 CHKERRQ(Adjoint2(user.U,arr,&user)); 514 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n")); 515 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1])); 516 break; 517 case 0 : 518 /* Create optimization context and set up */ 519 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 520 CHKERRQ(TaoSetType(tao,TAOBLMVM)); 521 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 522 523 if (mf) { 524 CHKERRQ(MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H)); 525 CHKERRQ(MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat)); 526 CHKERRQ(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE)); 527 CHKERRQ(TaoSetHessian(tao,user.H,user.H,MatrixFreeHessian,(void *)&user)); 528 } else { /* Create Hessian matrix */ 529 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.H)); 530 CHKERRQ(MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2)); 531 CHKERRQ(MatSetUp(user.H)); 532 CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 533 } 534 535 /* Not use any preconditioner */ 536 CHKERRQ(TaoGetKSP(tao,&ksp)); 537 if (ksp) { 538 CHKERRQ(KSPGetPC(ksp,&pc)); 539 CHKERRQ(PCSetType(pc,PCNONE)); 540 } 541 542 /* Set initial solution guess */ 543 CHKERRQ(TaoSetSolution(tao,x)); 544 CHKERRQ(TaoSetFromOptions(tao)); 545 CHKERRQ(TaoSolve(tao)); 546 CHKERRQ(TaoDestroy(&tao)); 547 CHKERRQ(MatDestroy(&user.H)); 548 break; 549 default: 550 break; 551 } 552 553 /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */ 554 CHKERRQ(MatDestroy(&user.A)); 555 CHKERRQ(VecDestroy(&user.U)); 556 CHKERRQ(VecDestroy(&user.Lambda[0])); 557 CHKERRQ(VecDestroy(&user.Lambda2[0])); 558 CHKERRQ(VecDestroy(&user.Ihp1[0])); 559 CHKERRQ(VecDestroy(&user.Dir)); 560 CHKERRQ(TSDestroy(&user.ts)); 561 CHKERRQ(VecDestroy(&x)); 562 ierr = PetscFinalize(); 563 return(ierr); 564 } 565 566 /*TEST 567 build: 568 requires: !complex !single 569 570 test: 571 args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125 572 output_file: output/ex20opt_ic_1.out 573 574 test: 575 suffix: 2 576 args: -ts_type beuler -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none 577 output_file: output/ex20opt_ic_2.out 578 579 test: 580 suffix: 3 581 args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 100 -ts_dt 0.01 -tao_type bntr -tao_bnk_pc_type none 582 output_file: output/ex20opt_ic_3.out 583 584 test: 585 suffix: 4 586 args: -implicitform 0 -ts_dt 0.01 -ts_max_steps 2 -ts_rhs_jacobian_test_mult_transpose -mat_shell_test_mult_transpose_view -ts_rhs_jacobian_test_mult -mat_shell_test_mult_view -mode 1 -my_tao_mf 587 TEST*/ 588