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 PetscCall(VecGetArrayRead(U,&u)); 51 PetscCall(VecGetArray(F,&f)); 52 f[0] = u[1]; 53 f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]); 54 PetscCall(VecRestoreArrayRead(U,&u)); 55 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 74 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 75 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 76 if (A != B) { 77 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 78 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 79 } 80 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 94 PetscCall(VecGetArrayRead(Vl[0],&vl)); 95 PetscCall(VecGetArrayRead(Vr,&vr)); 96 PetscCall(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 PetscCall(VecRestoreArrayRead(U,&u)); 108 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 109 PetscCall(VecRestoreArrayRead(Vr,&vr)); 110 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 124 PetscCall(VecGetArrayRead(Udot,&udot)); 125 PetscCall(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 PetscCall(VecRestoreArrayRead(U,&u)); 129 PetscCall(VecRestoreArrayRead(Udot,&udot)); 130 PetscCall(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 PetscCall(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 PetscCall(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 146 PetscCall(VecRestoreArrayRead(U,&u)); 147 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 149 if (A != B) { 150 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 151 PetscCall(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 PetscCall(TSGetTimeStep(ts,&dt)); 166 PetscCall(TSGetMaxTime(ts,&tfinal)); 167 168 while (user->next_output <= t && user->next_output <= tfinal) { 169 PetscCall(VecDuplicate(U,&interpolatedU)); 170 PetscCall(TSInterpolate(ts,user->next_output,interpolatedU)); 171 PetscCall(VecGetArrayRead(interpolatedU,&u)); 172 PetscCall(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 PetscCall(VecRestoreArrayRead(interpolatedU,&u)); 176 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 192 PetscCall(VecGetArrayRead(Vl[0],&vl)); 193 PetscCall(VecGetArrayRead(Vr,&vr)); 194 PetscCall(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 PetscCall(VecRestoreArrayRead(U,&u)); 205 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 206 PetscCall(VecRestoreArrayRead(Vr,&vr)); 207 PetscCall(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 PetscCall(VecCopy(IC,user_ptr->U)); /* set up the initial condition */ 222 223 PetscCall(TSSetTime(ts,0.0)); 224 PetscCall(TSSetStepNumber(ts,0)); 225 PetscCall(TSSetTimeStep(ts,0.001)); /* can be overwritten by command line options */ 226 PetscCall(TSSetFromOptions(ts)); 227 PetscCall(TSSolve(ts,user_ptr->U)); 228 229 PetscCall(VecGetArrayRead(user_ptr->U,&x_ptr)); 230 PetscCall(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 PetscCall(VecRestoreArray(user_ptr->Lambda[0],&y_ptr)); 235 PetscCall(VecRestoreArrayRead(user_ptr->U,&x_ptr)); 236 237 PetscCall(TSSetCostGradients(ts,1,user_ptr->Lambda,NULL)); 238 PetscCall(TSAdjointSolve(ts)); 239 PetscCall(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 PetscCall(VecCopy(U,user_ptr->U)); 253 PetscCall(VecGetArray(user_ptr->Dir,&x_ptr)); 254 x_ptr[0] = 1.; 255 x_ptr[1] = 0.; 256 PetscCall(VecRestoreArray(user_ptr->Dir,&x_ptr)); 257 PetscCall(Adjoint2(user_ptr->U,harr,user_ptr)); 258 col = 0; 259 PetscCall(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 260 261 PetscCall(VecCopy(U,user_ptr->U)); 262 PetscCall(VecGetArray(user_ptr->Dir,&x_ptr)); 263 x_ptr[0] = 0.; 264 x_ptr[1] = 1.; 265 PetscCall(VecRestoreArray(user_ptr->Dir,&x_ptr)); 266 PetscCall(Adjoint2(user_ptr->U,harr,user_ptr)); 267 col = 1; 268 PetscCall(MatSetValues(H,2,rows,1,&col,harr,INSERT_VALUES)); 269 270 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 271 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 272 if (H != Hpre) { 273 PetscCall(MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY)); 274 PetscCall(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 PetscCall(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 PetscCall(TSAdjointReset(ts)); 301 302 PetscCall(TSSetTime(ts,0.0)); 303 PetscCall(TSSetStepNumber(ts,0)); 304 PetscCall(TSSetTimeStep(ts,0.001)); 305 PetscCall(TSSetFromOptions(ts)); 306 PetscCall(TSSetCostHessianProducts(ts,1,ctx->Lambda2,NULL,ctx->Dir)); 307 PetscCall(TSAdjointSetForward(ts,NULL)); 308 PetscCall(TSSolve(ts,U)); 309 310 /* Set terminal conditions for first- and second-order adjonts */ 311 PetscCall(VecGetArray(U,&x_ptr)); 312 PetscCall(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 PetscCall(VecRestoreArray(ctx->Lambda[0],&y_ptr)); 316 PetscCall(VecRestoreArray(U,&x_ptr)); 317 PetscCall(TSForwardGetSensitivities(ts,NULL,&tlmsen)); 318 PetscCall(MatDenseGetColumn(tlmsen,0,&x_ptr)); 319 PetscCall(VecGetArray(ctx->Lambda2[0],&y_ptr)); 320 y_ptr[0] = 2.*x_ptr[0]; 321 y_ptr[1] = 2.*x_ptr[1]; 322 PetscCall(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 323 PetscCall(MatDenseRestoreColumn(tlmsen,&x_ptr)); 324 325 PetscCall(TSSetCostGradients(ts,1,ctx->Lambda,NULL)); 326 if (ctx->implicitform) { 327 PetscCall(TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 328 } else { 329 PetscCall(TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,NULL,NULL,NULL,NULL,NULL,NULL,ctx)); 330 } 331 PetscCall(TSAdjointSolve(ts)); 332 333 PetscCall(VecGetArray(ctx->Lambda2[0],&x_ptr)); 334 arr[0] = x_ptr[0]; 335 arr[1] = x_ptr[1]; 336 PetscCall(VecRestoreArray(ctx->Lambda2[0],&x_ptr)); 337 338 PetscCall(TSAdjointReset(ts)); 339 PetscCall(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 PetscCall(VecDuplicate(U,&Up)); 353 PetscCall(VecDuplicate(U,&G)); 354 PetscCall(VecDuplicate(U,&Gp)); 355 356 PetscCall(FormFunctionGradient(tao,U,&f,G,ctx)); 357 358 PetscCall(VecCopy(U,Up)); 359 PetscCall(VecGetArray(Up,&u)); 360 u[0] += eps; 361 PetscCall(VecRestoreArray(Up,&u)); 362 PetscCall(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 363 PetscCall(VecAXPY(Gp,-1,G)); 364 PetscCall(VecScale(Gp,1./eps)); 365 PetscCall(VecGetArray(Gp,&u)); 366 arr[0] = u[0]; 367 arr[1] = u[1]; 368 PetscCall(VecRestoreArray(Gp,&u)); 369 370 PetscCall(VecCopy(U,Up)); 371 PetscCall(VecGetArray(Up,&u)); 372 u[1] += eps; 373 PetscCall(VecRestoreArray(Up,&u)); 374 PetscCall(FormFunctionGradient(tao,Up,&f,Gp,ctx)); 375 PetscCall(VecAXPY(Gp,-1,G)); 376 PetscCall(VecScale(Gp,1./eps)); 377 PetscCall(VecGetArray(Gp,&u)); 378 arr[2] = u[0]; 379 arr[3] = u[1]; 380 PetscCall(VecRestoreArray(Gp,&u)); 381 382 PetscCall(VecDestroy(&G)); 383 PetscCall(VecDestroy(&Gp)); 384 PetscCall(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 PetscCall(MatShellGetContext(mat,&user_ptr)); 395 PetscCall(VecCopy(svec,user_ptr->Dir)); 396 PetscCall(VecGetArray(y,&y_ptr)); 397 PetscCall(Adjoint2(user_ptr->U,y_ptr,user_ptr)); 398 PetscCall(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 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 417 PetscCallMPI(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 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 428 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 429 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL)); 430 PetscCall(PetscOptionsGetReal(NULL,NULL,"-ic1",&ic1,NULL)); 431 PetscCall(PetscOptionsGetReal(NULL,NULL,"-ic2",&ic2,NULL)); 432 PetscCall(PetscOptionsGetBool(NULL,NULL,"-my_tao_mf",&mf,NULL)); /* matrix-free hessian for optimization */ 433 PetscCall(PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL)); 434 435 /* Create necessary matrix and vectors, solve same ODE on every process */ 436 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.A)); 437 PetscCall(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 438 PetscCall(MatSetFromOptions(user.A)); 439 PetscCall(MatSetUp(user.A)); 440 PetscCall(MatCreateVecs(user.A,&user.U,NULL)); 441 PetscCall(MatCreateVecs(user.A,&user.Dir,NULL)); 442 PetscCall(MatCreateVecs(user.A,&user.Lambda[0],NULL)); 443 PetscCall(MatCreateVecs(user.A,&user.Lambda2[0],NULL)); 444 PetscCall(MatCreateVecs(user.A,&user.Ihp1[0],NULL)); 445 446 /* Create timestepping solver context */ 447 PetscCall(TSCreate(PETSC_COMM_WORLD,&user.ts)); 448 PetscCall(TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 449 if (user.implicitform) { 450 PetscCall(TSSetIFunction(user.ts,NULL,IFunction,&user)); 451 PetscCall(TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user)); 452 PetscCall(TSSetType(user.ts,TSCN)); 453 } else { 454 PetscCall(TSSetRHSFunction(user.ts,NULL,RHSFunction,&user)); 455 PetscCall(TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user)); 456 PetscCall(TSSetType(user.ts,TSRK)); 457 } 458 PetscCall(TSSetMaxTime(user.ts,user.ftime)); 459 PetscCall(TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP)); 460 461 if (monitor) { 462 PetscCall(TSMonitorSet(user.ts,Monitor,&user,NULL)); 463 } 464 465 /* Set ODE initial conditions */ 466 PetscCall(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 PetscCall(VecRestoreArray(user.U,&x_ptr)); 470 471 /* Set runtime options */ 472 PetscCall(TSSetFromOptions(user.ts)); 473 474 /* Obtain the observation */ 475 PetscCall(TSSolve(user.ts,user.U)); 476 PetscCall(VecGetArray(user.U,&x_ptr)); 477 user.ob[0] = x_ptr[0]; 478 user.ob[1] = x_ptr[1]; 479 PetscCall(VecRestoreArray(user.U,&x_ptr)); 480 481 PetscCall(VecDuplicate(user.U,&x)); 482 PetscCall(VecGetArray(x,&x_ptr)); 483 x_ptr[0] = ic1; 484 x_ptr[1] = ic2; 485 PetscCall(VecRestoreArray(x,&x_ptr)); 486 487 /* Save trajectory of solution so that TSAdjointSolve() may be used */ 488 PetscCall(TSSetSaveTrajectory(user.ts)); 489 490 /* Compare finite difference and second-order adjoint. */ 491 switch (mode) { 492 case 2 : 493 PetscCall(FiniteDiff(x,arr,&user)); 494 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n Finite difference approximation of the Hessian\n")); 495 PetscCall(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 PetscCall(VecCopy(x,user.U)); 499 PetscCall(VecGetArray(user.Dir,&x_ptr)); 500 x_ptr[0] = 1.; 501 x_ptr[1] = 0.; 502 PetscCall(VecRestoreArray(user.Dir,&x_ptr)); 503 PetscCall(Adjoint2(user.U,arr,&user)); 504 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nFirst column of the Hessian\n")); 505 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%g\n%g\n",(double)arr[0],(double)arr[1])); 506 PetscCall(VecCopy(x,user.U)); 507 PetscCall(VecGetArray(user.Dir,&x_ptr)); 508 x_ptr[0] = 0.; 509 x_ptr[1] = 1.; 510 PetscCall(VecRestoreArray(user.Dir,&x_ptr)); 511 PetscCall(Adjoint2(user.U,arr,&user)); 512 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSecond column of the Hessian\n")); 513 PetscCall(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 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 518 PetscCall(TaoSetType(tao,TAOBLMVM)); 519 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 520 521 if (mf) { 522 PetscCall(MatCreateShell(PETSC_COMM_SELF,2,2,2,2,(void*)&user,&user.H)); 523 PetscCall(MatShellSetOperation(user.H,MATOP_MULT,(void(*)(void))HessianProductMat)); 524 PetscCall(MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE)); 525 PetscCall(TaoSetHessian(tao,user.H,user.H,MatrixFreeHessian,(void *)&user)); 526 } else { /* Create Hessian matrix */ 527 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.H)); 528 PetscCall(MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,2,2)); 529 PetscCall(MatSetUp(user.H)); 530 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 531 } 532 533 /* Not use any preconditioner */ 534 PetscCall(TaoGetKSP(tao,&ksp)); 535 if (ksp) { 536 PetscCall(KSPGetPC(ksp,&pc)); 537 PetscCall(PCSetType(pc,PCNONE)); 538 } 539 540 /* Set initial solution guess */ 541 PetscCall(TaoSetSolution(tao,x)); 542 PetscCall(TaoSetFromOptions(tao)); 543 PetscCall(TaoSolve(tao)); 544 PetscCall(TaoDestroy(&tao)); 545 PetscCall(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 PetscCall(MatDestroy(&user.A)); 553 PetscCall(VecDestroy(&user.U)); 554 PetscCall(VecDestroy(&user.Lambda[0])); 555 PetscCall(VecDestroy(&user.Lambda2[0])); 556 PetscCall(VecDestroy(&user.Ihp1[0])); 557 PetscCall(VecDestroy(&user.Dir)); 558 PetscCall(TSDestroy(&user.ts)); 559 PetscCall(VecDestroy(&x)); 560 PetscCall(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