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