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 const PetscScalar *u; 160 PetscReal tfinal, dt; 161 User user = (User)ctx; 162 Vec interpolatedU; 163 164 PetscFunctionBeginUser; 165 CHKERRQ(TSGetTimeStep(ts,&dt)); 166 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 167 168 while (user->next_output <= t && user->next_output <= tfinal) { 169 CHKERRQ(VecDuplicate(U,&interpolatedU)); 170 CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedU)); 171 CHKERRQ(VecGetArrayRead(interpolatedU,&u)); 172 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n", 173 (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(u[0]), 174 (double)PetscRealPart(u[1]))); 175 CHKERRQ(VecRestoreArrayRead(interpolatedU,&u)); 176 CHKERRQ(VecDestroy(&interpolatedU)); 177 user->next_output += 0.1; 178 } 179 PetscFunctionReturn(0); 180 } 181 182 static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 183 { 184 const PetscScalar *vl,*vr,*u; 185 PetscScalar *vhv; 186 PetscScalar dJdU[2][2][2]={{{0}}}; 187 PetscInt i,j,k; 188 User user = (User)ctx; 189 190 PetscFunctionBeginUser; 191 CHKERRQ(VecGetArrayRead(U,&u)); 192 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 193 CHKERRQ(VecGetArrayRead(Vr,&vr)); 194 CHKERRQ(VecGetArray(VHV[0],&vhv)); 195 dJdU[1][0][0] = 2.*user->mu*u[1]; 196 dJdU[1][1][0] = 2.*user->mu*u[0]; 197 dJdU[1][0][1] = 2.*user->mu*u[0]; 198 for (j=0;j<2;j++) { 199 vhv[j] = 0; 200 for (k=0;k<2;k++) 201 for (i=0;i<2;i++) 202 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 203 } 204 CHKERRQ(VecRestoreArrayRead(U,&u)); 205 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 206 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 207 CHKERRQ(VecRestoreArray(VHV[0],&vhv)); 208 PetscFunctionReturn(0); 209 } 210 211 /* ------------------ User-defined routines for TAO -------------------------- */ 212 213 static PetscErrorCode FormFunctionGradient(Tao tao,Vec IC,PetscReal *f,Vec G,void *ctx) 214 { 215 User user_ptr = (User)ctx; 216 TS ts = user_ptr->ts; 217 const PetscScalar *x_ptr; 218 PetscScalar *y_ptr; 219 220 PetscFunctionBeginUser; 221 CHKERRQ(VecCopy(IC,user_ptr->U)); /* set up the initial condition */ 222 223 CHKERRQ(TSSetTime(ts,0.0)); 224 CHKERRQ(TSSetStepNumber(ts,0)); 225 CHKERRQ(TSSetTimeStep(ts,0.001)); /* can be overwritten by command line options */ 226 CHKERRQ(TSSetFromOptions(ts)); 227 CHKERRQ(TSSolve(ts,user_ptr->U)); 228 229 CHKERRQ(VecGetArrayRead(user_ptr->U,&x_ptr)); 230 CHKERRQ(VecGetArray(user_ptr->Lambda[0],&y_ptr)); 231 *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]); 232 y_ptr[0] = 2.*(x_ptr[0]-user_ptr->ob[0]); 233 y_ptr[1] = 2.*(x_ptr[1]-user_ptr->ob[1]); 234 CHKERRQ(VecRestoreArray(user_ptr->Lambda[0],&y_ptr)); 235 CHKERRQ(VecRestoreArrayRead(user_ptr->U,&x_ptr)); 236 237 CHKERRQ(TSSetCostGradients(ts,1,user_ptr->Lambda,NULL)); 238 CHKERRQ(TSAdjointSolve(ts)); 239 CHKERRQ(VecCopy(user_ptr->Lambda[0],G)); 240 PetscFunctionReturn(0); 241 } 242 243 static PetscErrorCode FormHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx) 244 { 245 User user_ptr = (User)ctx; 246 PetscScalar harr[2]; 247 PetscScalar *x_ptr; 248 const PetscInt rows[2] = {0,1}; 249 PetscInt col; 250 251 PetscFunctionBeginUser; 252 CHKERRQ(VecCopy(U,user_ptr->U)); 253 CHKERRQ(VecGetArray(user_ptr->Dir,&x_ptr)); 254 x_ptr[0] = 1.; 255 x_ptr[1] = 0.; 256 CHKERRQ(VecRestoreArray(user_ptr->Dir,&x_ptr)); 257 CHKERRQ(Adjoint2(user_ptr->U,harr,user_ptr)); 258 col = 0; 259 CHKERRQ(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 260 261 CHKERRQ(VecCopy(U,user_ptr->U)); 262 CHKERRQ(VecGetArray(user_ptr->Dir,&x_ptr)); 263 x_ptr[0] = 0.; 264 x_ptr[1] = 1.; 265 CHKERRQ(VecRestoreArray(user_ptr->Dir,&x_ptr)); 266 CHKERRQ(Adjoint2(user_ptr->U,harr,user_ptr)); 267 col = 1; 268 CHKERRQ(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 269 270 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 271 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 272 if (H != Hpre) { 273 CHKERRQ(MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY)); 274 CHKERRQ(MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY)); 275 } 276 PetscFunctionReturn(0); 277 } 278 279 static PetscErrorCode MatrixFreeHessian(Tao tao,Vec U,Mat H,Mat Hpre,void *ctx) 280 { 281 User user_ptr = (User)ctx; 282 283 PetscFunctionBeginUser; 284 CHKERRQ(VecCopy(U,user_ptr->U)); 285 PetscFunctionReturn(0); 286 } 287 288 /* ------------ Routines calculating second-order derivatives -------------- */ 289 290 /* 291 Compute the Hessian-vector product for the cost function using Second-order adjoint 292 */ 293 PetscErrorCode Adjoint2(Vec U,PetscScalar arr[],User ctx) 294 { 295 TS ts = ctx->ts; 296 PetscScalar *x_ptr,*y_ptr; 297 Mat tlmsen; 298 299 PetscFunctionBeginUser; 300 CHKERRQ(TSAdjointReset(ts)); 301 302 CHKERRQ(TSSetTime(ts,0.0)); 303 CHKERRQ(TSSetStepNumber(ts,0)); 304 CHKERRQ(TSSetTimeStep(ts,0.001)); 305 CHKERRQ(TSSetFromOptions(ts)); 306 CHKERRQ(TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir)); 307 CHKERRQ(TSAdjointSetForward(ts,NULL)); 308 CHKERRQ(TSSolve(ts,U)); 309 310 /* Set terminal conditions for first- and second-order adjonts */ 311 CHKERRQ(VecGetArray(U,&x_ptr)); 312 CHKERRQ(VecGetArray(ctx->Lambda[0],&y_ptr)); 313 y_ptr[0] = 2.*(x_ptr[0]-ctx->ob[0]); 314 y_ptr[1] = 2.*(x_ptr[1]-ctx->ob[1]); 315 CHKERRQ(VecRestoreArray(ctx->Lambda[0],&y_ptr)); 316 CHKERRQ(VecRestoreArray(U,&x_ptr)); 317 CHKERRQ(TSForwardGetSensitivities(ts,NULL,&tlmsen)); 318 CHKERRQ(MatDenseGetColumn(tlmsen,0,&x_ptr)); 319 CHKERRQ(VecGetArray(ctx->Lambda2[0],&y_ptr)); 320 y_ptr[0] = 2.*x_ptr[0]; 321 y_ptr[1] = 2.*x_ptr[1]; 322 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 323 CHKERRQ(MatDenseRestoreColumn(tlmsen,&x_ptr)); 324 325 CHKERRQ(TSSetCostGradients(ts,1,ctx->Lambda,NULL)); 326 if (ctx->implicitform) { 327 CHKERRQ(TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 328 } else { 329 CHKERRQ(TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 330 } 331 CHKERRQ(TSAdjointSolve(ts)); 332 333 CHKERRQ(VecGetArray(ctx->Lambda2[0],&x_ptr)); 334 arr[0] = x_ptr[0]; 335 arr[1] = x_ptr[1]; 336 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&x_ptr)); 337 338 CHKERRQ(TSAdjointReset(ts)); 339 CHKERRQ(TSAdjointResetForward(ts)); 340 PetscFunctionReturn(0); 341 } 342 343 PetscErrorCode FiniteDiff(Vec U,PetscScalar arr[],User ctx) 344 { 345 Vec Up,G,Gp; 346 const PetscScalar eps = PetscRealConstant(1e-7); 347 PetscScalar *u; 348 Tao tao = NULL; 349 PetscReal f; 350 351 PetscFunctionBeginUser; 352 CHKERRQ(VecDuplicate(U,&Up)); 353 CHKERRQ(VecDuplicate(U,&G)); 354 CHKERRQ(VecDuplicate(U,&Gp)); 355 356 CHKERRQ(FormFunctionGradient(tao,U,&f,G,ctx)); 357 358 CHKERRQ(VecCopy(U,Up)); 359 CHKERRQ(VecGetArray(Up,&u)); 360 u[0] += eps; 361 CHKERRQ(VecRestoreArray(Up,&u)); 362 CHKERRQ(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 363 CHKERRQ(VecAXPY(Gp,-1,G)); 364 CHKERRQ(VecScale(Gp,1./eps)); 365 CHKERRQ(VecGetArray(Gp,&u)); 366 arr[0] = u[0]; 367 arr[1] = u[1]; 368 CHKERRQ(VecRestoreArray(Gp,&u)); 369 370 CHKERRQ(VecCopy(U,Up)); 371 CHKERRQ(VecGetArray(Up,&u)); 372 u[1] += eps; 373 CHKERRQ(VecRestoreArray(Up,&u)); 374 CHKERRQ(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 375 CHKERRQ(VecAXPY(Gp,-1,G)); 376 CHKERRQ(VecScale(Gp,1./eps)); 377 CHKERRQ(VecGetArray(Gp,&u)); 378 arr[2] = u[0]; 379 arr[3] = u[1]; 380 CHKERRQ(VecRestoreArray(Gp,&u)); 381 382 CHKERRQ(VecDestroy(&G)); 383 CHKERRQ(VecDestroy(&Gp)); 384 CHKERRQ(VecDestroy(&Up)); 385 PetscFunctionReturn(0); 386 } 387 388 static PetscErrorCode HessianProductMat(Mat mat,Vec svec,Vec y) 389 { 390 User user_ptr; 391 PetscScalar *y_ptr; 392 393 PetscFunctionBeginUser; 394 CHKERRQ(MatShellGetContext(mat,&user_ptr)); 395 CHKERRQ(VecCopy(svec,user_ptr->Dir)); 396 CHKERRQ(VecGetArray(y,&y_ptr)); 397 CHKERRQ(Adjoint2(user_ptr->U,y_ptr,user_ptr)); 398 CHKERRQ(VecRestoreArray(y,&y_ptr)); 399 PetscFunctionReturn(0); 400 } 401 402 int main(int argc,char **argv) 403 { 404 PetscBool monitor = PETSC_FALSE,mf = PETSC_TRUE; 405 PetscInt mode = 0; 406 PetscMPIInt size; 407 struct _n_User user; 408 Vec x; /* working vector for TAO */ 409 PetscScalar *x_ptr,arr[4]; 410 PetscScalar ic1 = 2.2,ic2 = -0.7; /* initial guess for TAO */ 411 Tao tao; 412 KSP ksp; 413 PC pc; 414 415 /* Initialize program */ 416 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 417 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 418 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 419 420 /* Set runtime options */ 421 user.next_output = 0.0; 422 user.mu = 1.0e3; 423 user.steps = 0; 424 user.ftime = 0.5; 425 user.implicitform = PETSC_TRUE; 426 427 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 428 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 429 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL)); 430 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL)); 431 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL)); 432 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL)); /* matrix-free hessian for optimization */ 433 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL)); 434 435 /* Create necessary matrix and vectors, solve same ODE on every process */ 436 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.A)); 437 CHKERRQ(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 438 CHKERRQ(MatSetFromOptions(user.A)); 439 CHKERRQ(MatSetUp(user.A)); 440 CHKERRQ(MatCreateVecs(user.A,&user.U,NULL)); 441 CHKERRQ(MatCreateVecs(user.A,&user.Dir,NULL)); 442 CHKERRQ(MatCreateVecs(user.A,&user.Lambda[0],NULL)); 443 CHKERRQ(MatCreateVecs(user.A,&user.Lambda2[0],NULL)); 444 CHKERRQ(MatCreateVecs(user.A,&user.Ihp1[0],NULL)); 445 446 /* Create timestepping solver context */ 447 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&user.ts)); 448 CHKERRQ(TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 449 if (user.implicitform) { 450 CHKERRQ(TSSetIFunction(user.ts,NULL,IFunction,&user)); 451 CHKERRQ(TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user)); 452 CHKERRQ(TSSetType(user.ts,TSCN)); 453 } else { 454 CHKERRQ(TSSetRHSFunction(user.ts,NULL,RHSFunction,&user)); 455 CHKERRQ(TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user)); 456 CHKERRQ(TSSetType(user.ts,TSRK)); 457 } 458 CHKERRQ(TSSetMaxTime(user.ts,user.ftime)); 459 CHKERRQ(TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP)); 460 461 if (monitor) { 462 CHKERRQ(TSMonitorSet(user.ts,Monitor,&user,NULL)); 463 } 464 465 /* Set ODE initial conditions */ 466 CHKERRQ(VecGetArray(user.U,&x_ptr)); 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 CHKERRQ(VecRestoreArray(user.U,&x_ptr)); 470 471 /* Set runtime options */ 472 CHKERRQ(TSSetFromOptions(user.ts)); 473 474 /* Obtain the observation */ 475 CHKERRQ(TSSolve(user.ts,user.U)); 476 CHKERRQ(VecGetArray(user.U,&x_ptr)); 477 user.ob[0] = x_ptr[0]; 478 user.ob[1] = x_ptr[1]; 479 CHKERRQ(VecRestoreArray(user.U,&x_ptr)); 480 481 CHKERRQ(VecDuplicate(user.U,&x)); 482 CHKERRQ(VecGetArray(x,&x_ptr)); 483 x_ptr[0] = ic1; 484 x_ptr[1] = ic2; 485 CHKERRQ(VecRestoreArray(x,&x_ptr)); 486 487 /* Save trajectory of solution so that TSAdjointSolve() may be used */ 488 CHKERRQ(TSSetSaveTrajectory(user.ts)); 489 490 /* Compare finite difference and second-order adjoint. */ 491 switch (mode) { 492 case 2 : 493 CHKERRQ(FiniteDiff(x,arr,&user)); 494 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n")); 495 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g %g\n%g %g\n",(double)arr[0],(double)arr[1],(double)arr[2],(double)arr[3])); 496 break; 497 case 1 : /* Compute the Hessian column by column */ 498 CHKERRQ(VecCopy(x,user.U)); 499 CHKERRQ(VecGetArray(user.Dir,&x_ptr)); 500 x_ptr[0] = 1.; 501 x_ptr[1] = 0.; 502 CHKERRQ(VecRestoreArray(user.Dir,&x_ptr)); 503 CHKERRQ(Adjoint2(user.U,arr,&user)); 504 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n")); 505 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1])); 506 CHKERRQ(VecCopy(x,user.U)); 507 CHKERRQ(VecGetArray(user.Dir,&x_ptr)); 508 x_ptr[0] = 0.; 509 x_ptr[1] = 1.; 510 CHKERRQ(VecRestoreArray(user.Dir,&x_ptr)); 511 CHKERRQ(Adjoint2(user.U,arr,&user)); 512 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n")); 513 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1])); 514 break; 515 case 0 : 516 /* Create optimization context and set up */ 517 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 518 CHKERRQ(TaoSetType(tao,TAOBLMVM)); 519 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 520 521 if (mf) { 522 CHKERRQ(MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H)); 523 CHKERRQ(MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat)); 524 CHKERRQ(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE)); 525 CHKERRQ(TaoSetHessian(tao,user.H,user.H,MatrixFreeHessian,(void *)&user)); 526 } else { /* Create Hessian matrix */ 527 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.H)); 528 CHKERRQ(MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2)); 529 CHKERRQ(MatSetUp(user.H)); 530 CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 531 } 532 533 /* Not use any preconditioner */ 534 CHKERRQ(TaoGetKSP(tao,&ksp)); 535 if (ksp) { 536 CHKERRQ(KSPGetPC(ksp,&pc)); 537 CHKERRQ(PCSetType(pc,PCNONE)); 538 } 539 540 /* Set initial solution guess */ 541 CHKERRQ(TaoSetSolution(tao,x)); 542 CHKERRQ(TaoSetFromOptions(tao)); 543 CHKERRQ(TaoSolve(tao)); 544 CHKERRQ(TaoDestroy(&tao)); 545 CHKERRQ(MatDestroy(&user.H)); 546 break; 547 default: 548 break; 549 } 550 551 /* Free work space. All PETSc objects should be destroyed when they are no longer needed. */ 552 CHKERRQ(MatDestroy(&user.A)); 553 CHKERRQ(VecDestroy(&user.U)); 554 CHKERRQ(VecDestroy(&user.Lambda[0])); 555 CHKERRQ(VecDestroy(&user.Lambda2[0])); 556 CHKERRQ(VecDestroy(&user.Ihp1[0])); 557 CHKERRQ(VecDestroy(&user.Dir)); 558 CHKERRQ(TSDestroy(&user.ts)); 559 CHKERRQ(VecDestroy(&x)); 560 CHKERRQ(PetscFinalize()); 561 return 0; 562 } 563 564 /*TEST 565 build: 566 requires: !complex !single 567 568 test: 569 args: -ts_type cn -viewer_binary_skip_info -tao_monitor -tao_view -mu 1000 -ts_dt 0.03125 570 output_file: output/ex20opt_ic_1.out 571 572 test: 573 suffix: 2 574 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 575 output_file: output/ex20opt_ic_2.out 576 577 test: 578 suffix: 3 579 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 580 output_file: output/ex20opt_ic_3.out 581 582 test: 583 suffix: 4 584 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 585 TEST*/ 586