1 2 static char help[] = "Solves the van der Pol equation.\n\ 3 Input parameters include:\n"; 4 5 /* 6 Concepts: TS^time-dependent nonlinear problems 7 Concepts: TS^van der Pol equation DAE equivalent 8 Concepts: TS^Optimization using adjoint sensitivity analysis 9 Processors: 1 10 */ 11 /* ------------------------------------------------------------------------ 12 13 Notes: 14 This code demonstrates how to solve a DAE-constrained optimization problem with TAO, TSAdjoint and TS. 15 The nonlinear problem is written in a DAE equivalent form. 16 The objective is to minimize the difference between observation and model prediction by finding an optimal value for parameter \mu. 17 The gradient is computed with the discrete adjoint of an implicit theta method, see ex20adj.c for details. 18 ------------------------------------------------------------------------- */ 19 #include <petsctao.h> 20 #include <petscts.h> 21 22 typedef struct _n_User *User; 23 struct _n_User { 24 TS ts; 25 PetscReal mu; 26 PetscReal next_output; 27 28 /* Sensitivity analysis support */ 29 PetscReal ftime; 30 Mat A; /* Jacobian matrix */ 31 Mat Jacp; /* JacobianP matrix */ 32 Mat H; /* Hessian matrix for optimization */ 33 Vec U,Lambda[1],Mup[1]; /* adjoint variables */ 34 Vec Lambda2[1],Mup2[1]; /* second-order adjoint variables */ 35 Vec Ihp1[1]; /* working space for Hessian evaluations */ 36 Vec Ihp2[1]; /* working space for Hessian evaluations */ 37 Vec Ihp3[1]; /* working space for Hessian evaluations */ 38 Vec Ihp4[1]; /* working space for Hessian evaluations */ 39 Vec Dir; /* direction vector */ 40 PetscReal ob[2]; /* observation used by the cost function */ 41 PetscBool implicitform; /* implicit ODE? */ 42 }; 43 44 PetscErrorCode FormFunctionGradient(Tao,Vec,PetscReal*,Vec,void*); 45 PetscErrorCode FormHessian(Tao,Vec,Mat,Mat,void*); 46 PetscErrorCode Adjoint2(Vec,PetscScalar[],User); 47 48 /* ----------------------- Explicit form of the ODE -------------------- */ 49 50 static PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec U,Vec F,void *ctx) 51 { 52 User user = (User)ctx; 53 PetscScalar *f; 54 const PetscScalar *u; 55 56 PetscFunctionBeginUser; 57 PetscCall(VecGetArrayRead(U,&u)); 58 PetscCall(VecGetArray(F,&f)); 59 f[0] = u[1]; 60 f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]); 61 PetscCall(VecRestoreArrayRead(U,&u)); 62 PetscCall(VecRestoreArray(F,&f)); 63 PetscFunctionReturn(0); 64 } 65 66 static PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,void *ctx) 67 { 68 User user = (User)ctx; 69 PetscReal mu = user->mu; 70 PetscInt rowcol[] = {0,1}; 71 PetscScalar J[2][2]; 72 const PetscScalar *u; 73 74 PetscFunctionBeginUser; 75 PetscCall(VecGetArrayRead(U,&u)); 76 77 J[0][0] = 0; 78 J[1][0] = -mu*(2.0*u[1]*u[0]+1.); 79 J[0][1] = 1.0; 80 J[1][1] = mu*(1.0-u[0]*u[0]); 81 PetscCall(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 82 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 83 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 84 if (B && A != B) { 85 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 86 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 87 } 88 89 PetscCall(VecRestoreArrayRead(U,&u)); 90 PetscFunctionReturn(0); 91 } 92 93 static PetscErrorCode RHSHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 94 { 95 const PetscScalar *vl,*vr,*u; 96 PetscScalar *vhv; 97 PetscScalar dJdU[2][2][2]={{{0}}}; 98 PetscInt i,j,k; 99 User user = (User)ctx; 100 101 PetscFunctionBeginUser; 102 PetscCall(VecGetArrayRead(U,&u)); 103 PetscCall(VecGetArrayRead(Vl[0],&vl)); 104 PetscCall(VecGetArrayRead(Vr,&vr)); 105 PetscCall(VecGetArray(VHV[0],&vhv)); 106 107 dJdU[1][0][0] = -2.*user->mu*u[1]; 108 dJdU[1][1][0] = -2.*user->mu*u[0]; 109 dJdU[1][0][1] = -2.*user->mu*u[0]; 110 for (j=0; j<2; j++) { 111 vhv[j] = 0; 112 for (k=0; k<2; k++) 113 for (i=0; i<2; i++) 114 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 115 } 116 117 PetscCall(VecRestoreArrayRead(U,&u)); 118 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 119 PetscCall(VecRestoreArrayRead(Vr,&vr)); 120 PetscCall(VecRestoreArray(VHV[0],&vhv)); 121 PetscFunctionReturn(0); 122 } 123 124 static PetscErrorCode RHSHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 125 { 126 const PetscScalar *vl,*vr,*u; 127 PetscScalar *vhv; 128 PetscScalar dJdP[2][2][1]={{{0}}}; 129 PetscInt i,j,k; 130 131 PetscFunctionBeginUser; 132 PetscCall(VecGetArrayRead(U,&u)); 133 PetscCall(VecGetArrayRead(Vl[0],&vl)); 134 PetscCall(VecGetArrayRead(Vr,&vr)); 135 PetscCall(VecGetArray(VHV[0],&vhv)); 136 137 dJdP[1][0][0] = -(1.+2.*u[0]*u[1]); 138 dJdP[1][1][0] = 1.-u[0]*u[0]; 139 for (j=0; j<2; j++) { 140 vhv[j] = 0; 141 for (k=0; k<1; k++) 142 for (i=0; i<2; i++) 143 vhv[j] += vl[i]*dJdP[i][j][k]*vr[k]; 144 } 145 146 PetscCall(VecRestoreArrayRead(U,&u)); 147 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 148 PetscCall(VecRestoreArrayRead(Vr,&vr)); 149 PetscCall(VecRestoreArray(VHV[0],&vhv)); 150 PetscFunctionReturn(0); 151 } 152 153 static PetscErrorCode RHSHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 154 { 155 const PetscScalar *vl,*vr,*u; 156 PetscScalar *vhv; 157 PetscScalar dJdU[2][1][2]={{{0}}}; 158 PetscInt i,j,k; 159 160 PetscFunctionBeginUser; 161 PetscCall(VecGetArrayRead(U,&u)); 162 PetscCall(VecGetArrayRead(Vl[0],&vl)); 163 PetscCall(VecGetArrayRead(Vr,&vr)); 164 PetscCall(VecGetArray(VHV[0],&vhv)); 165 166 dJdU[1][0][0] = -1.-2.*u[1]*u[0]; 167 dJdU[1][0][1] = 1.-u[0]*u[0]; 168 for (j=0; j<1; j++) { 169 vhv[j] = 0; 170 for (k=0; k<2; k++) 171 for (i=0; i<2; i++) 172 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 173 } 174 175 PetscCall(VecRestoreArrayRead(U,&u)); 176 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 177 PetscCall(VecRestoreArrayRead(Vr,&vr)); 178 PetscCall(VecRestoreArray(VHV[0],&vhv)); 179 PetscFunctionReturn(0); 180 } 181 182 static PetscErrorCode RHSHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 183 { 184 PetscFunctionBeginUser; 185 PetscFunctionReturn(0); 186 } 187 188 /* ----------------------- Implicit form of the ODE -------------------- */ 189 190 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,void *ctx) 191 { 192 User user = (User)ctx; 193 PetscScalar *f; 194 const PetscScalar *u,*udot; 195 196 PetscFunctionBeginUser; 197 PetscCall(VecGetArrayRead(U,&u)); 198 PetscCall(VecGetArrayRead(Udot,&udot)); 199 PetscCall(VecGetArray(F,&f)); 200 201 f[0] = udot[0] - u[1]; 202 f[1] = udot[1] - user->mu*((1.0-u[0]*u[0])*u[1] - u[0]) ; 203 204 PetscCall(VecRestoreArrayRead(U,&u)); 205 PetscCall(VecRestoreArrayRead(Udot,&udot)); 206 PetscCall(VecRestoreArray(F,&f)); 207 PetscFunctionReturn(0); 208 } 209 210 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,void *ctx) 211 { 212 User user = (User)ctx; 213 PetscInt rowcol[] = {0,1}; 214 PetscScalar J[2][2]; 215 const PetscScalar *u; 216 217 PetscFunctionBeginUser; 218 PetscCall(VecGetArrayRead(U,&u)); 219 220 J[0][0] = a; J[0][1] = -1.0; 221 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]); 222 PetscCall(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 223 PetscCall(VecRestoreArrayRead(U,&u)); 224 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 225 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 226 if (A != B) { 227 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 228 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 229 } 230 231 PetscCall(VecRestoreArrayRead(U,&u)); 232 PetscFunctionReturn(0); 233 } 234 235 static PetscErrorCode RHSJacobianP(TS ts,PetscReal t,Vec U,Mat A,void *ctx) 236 { 237 PetscInt row[] = {0,1},col[]={0}; 238 PetscScalar J[2][1]; 239 const PetscScalar *u; 240 241 PetscFunctionBeginUser; 242 PetscCall(VecGetArrayRead(U,&u)); 243 244 J[0][0] = 0; 245 J[1][0] = (1.-u[0]*u[0])*u[1]-u[0]; 246 PetscCall(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 247 PetscCall(VecRestoreArrayRead(U,&u)); 248 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 249 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 250 251 PetscCall(VecRestoreArrayRead(U,&u)); 252 PetscFunctionReturn(0); 253 } 254 255 static PetscErrorCode IHessianProductUU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 256 { 257 const PetscScalar *vl,*vr,*u; 258 PetscScalar *vhv; 259 PetscScalar dJdU[2][2][2]={{{0}}}; 260 PetscInt i,j,k; 261 User user = (User)ctx; 262 263 PetscFunctionBeginUser; 264 PetscCall(VecGetArrayRead(U,&u)); 265 PetscCall(VecGetArrayRead(Vl[0],&vl)); 266 PetscCall(VecGetArrayRead(Vr,&vr)); 267 PetscCall(VecGetArray(VHV[0],&vhv)); 268 269 dJdU[1][0][0] = 2.*user->mu*u[1]; 270 dJdU[1][1][0] = 2.*user->mu*u[0]; 271 dJdU[1][0][1] = 2.*user->mu*u[0]; 272 for (j=0; j<2; j++) { 273 vhv[j] = 0; 274 for (k=0; k<2; k++) 275 for (i=0; i<2; i++) 276 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 277 } 278 279 PetscCall(VecRestoreArrayRead(U,&u)); 280 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 281 PetscCall(VecRestoreArrayRead(Vr,&vr)); 282 PetscCall(VecRestoreArray(VHV[0],&vhv)); 283 PetscFunctionReturn(0); 284 } 285 286 static PetscErrorCode IHessianProductUP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 287 { 288 const PetscScalar *vl,*vr,*u; 289 PetscScalar *vhv; 290 PetscScalar dJdP[2][2][1]={{{0}}}; 291 PetscInt i,j,k; 292 293 PetscFunctionBeginUser; 294 PetscCall(VecGetArrayRead(U,&u)); 295 PetscCall(VecGetArrayRead(Vl[0],&vl)); 296 PetscCall(VecGetArrayRead(Vr,&vr)); 297 PetscCall(VecGetArray(VHV[0],&vhv)); 298 299 dJdP[1][0][0] = 1.+2.*u[0]*u[1]; 300 dJdP[1][1][0] = u[0]*u[0]-1.; 301 for (j=0; j<2; j++) { 302 vhv[j] = 0; 303 for (k=0; k<1; k++) 304 for (i=0; i<2; i++) 305 vhv[j] += vl[i]*dJdP[i][j][k]*vr[k]; 306 } 307 308 PetscCall(VecRestoreArrayRead(U,&u)); 309 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 310 PetscCall(VecRestoreArrayRead(Vr,&vr)); 311 PetscCall(VecRestoreArray(VHV[0],&vhv)); 312 PetscFunctionReturn(0); 313 } 314 315 static PetscErrorCode IHessianProductPU(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 316 { 317 const PetscScalar *vl,*vr,*u; 318 PetscScalar *vhv; 319 PetscScalar dJdU[2][1][2]={{{0}}}; 320 PetscInt i,j,k; 321 322 PetscFunctionBeginUser; 323 PetscCall(VecGetArrayRead(U,&u)); 324 PetscCall(VecGetArrayRead(Vl[0],&vl)); 325 PetscCall(VecGetArrayRead(Vr,&vr)); 326 PetscCall(VecGetArray(VHV[0],&vhv)); 327 328 dJdU[1][0][0] = 1.+2.*u[1]*u[0]; 329 dJdU[1][0][1] = u[0]*u[0]-1.; 330 for (j=0; j<1; j++) { 331 vhv[j] = 0; 332 for (k=0; k<2; k++) 333 for (i=0; i<2; i++) 334 vhv[j] += vl[i]*dJdU[i][j][k]*vr[k]; 335 } 336 337 PetscCall(VecRestoreArrayRead(U,&u)); 338 PetscCall(VecRestoreArrayRead(Vl[0],&vl)); 339 PetscCall(VecRestoreArrayRead(Vr,&vr)); 340 PetscCall(VecRestoreArray(VHV[0],&vhv)); 341 PetscFunctionReturn(0); 342 } 343 344 static PetscErrorCode IHessianProductPP(TS ts,PetscReal t,Vec U,Vec *Vl,Vec Vr,Vec *VHV,void *ctx) 345 { 346 PetscFunctionBeginUser; 347 PetscFunctionReturn(0); 348 } 349 350 /* Monitor timesteps and use interpolation to output at integer multiples of 0.1 */ 351 static PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal t,Vec X,void *ctx) 352 { 353 PetscErrorCode ierr; 354 const PetscScalar *x; 355 PetscReal tfinal, dt; 356 User user = (User)ctx; 357 Vec interpolatedX; 358 359 PetscFunctionBeginUser; 360 PetscCall(TSGetTimeStep(ts,&dt)); 361 PetscCall(TSGetMaxTime(ts,&tfinal)); 362 363 while (user->next_output <= t && user->next_output <= tfinal) { 364 PetscCall(VecDuplicate(X,&interpolatedX)); 365 PetscCall(TSInterpolate(ts,user->next_output,interpolatedX)); 366 PetscCall(VecGetArrayRead(interpolatedX,&x)); 367 ierr = PetscPrintf(PETSC_COMM_WORLD,"[%g] %D TS %g (dt = %g) X %g %g\n", 368 (double)user->next_output,step,(double)t,(double)dt,(double)PetscRealPart(x[0]), 369 (double)PetscRealPart(x[1]));PetscCall(ierr); 370 PetscCall(VecRestoreArrayRead(interpolatedX,&x)); 371 PetscCall(VecDestroy(&interpolatedX)); 372 user->next_output += PetscRealConstant(0.1); 373 } 374 PetscFunctionReturn(0); 375 } 376 377 int main(int argc,char **argv) 378 { 379 Vec P; 380 PetscBool monitor = PETSC_FALSE; 381 PetscScalar *x_ptr; 382 const PetscScalar *y_ptr; 383 PetscMPIInt size; 384 struct _n_User user; 385 Tao tao; 386 KSP ksp; 387 PC pc; 388 389 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 390 Initialize program 391 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 392 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 393 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 394 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 395 396 /* Create TAO solver and set desired solution method */ 397 PetscCall(TaoCreate(PETSC_COMM_WORLD,&tao)); 398 PetscCall(TaoSetType(tao,TAOBQNLS)); 399 400 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 401 Set runtime options 402 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 403 user.next_output = 0.0; 404 user.mu = PetscRealConstant(1.0e3); 405 user.ftime = PetscRealConstant(0.5); 406 user.implicitform = PETSC_TRUE; 407 408 PetscCall(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 409 PetscCall(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 410 PetscCall(PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL)); 411 412 /* Create necessary matrix and vectors, solve same ODE on every process */ 413 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.A)); 414 PetscCall(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 415 PetscCall(MatSetFromOptions(user.A)); 416 PetscCall(MatSetUp(user.A)); 417 PetscCall(MatCreateVecs(user.A,&user.U,NULL)); 418 PetscCall(MatCreateVecs(user.A,&user.Lambda[0],NULL)); 419 PetscCall(MatCreateVecs(user.A,&user.Lambda2[0],NULL)); 420 PetscCall(MatCreateVecs(user.A,&user.Ihp1[0],NULL)); 421 PetscCall(MatCreateVecs(user.A,&user.Ihp2[0],NULL)); 422 423 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 424 PetscCall(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1)); 425 PetscCall(MatSetFromOptions(user.Jacp)); 426 PetscCall(MatSetUp(user.Jacp)); 427 PetscCall(MatCreateVecs(user.Jacp,&user.Dir,NULL)); 428 PetscCall(MatCreateVecs(user.Jacp,&user.Mup[0],NULL)); 429 PetscCall(MatCreateVecs(user.Jacp,&user.Mup2[0],NULL)); 430 PetscCall(MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL)); 431 PetscCall(MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL)); 432 433 /* Create timestepping solver context */ 434 PetscCall(TSCreate(PETSC_COMM_WORLD,&user.ts)); 435 PetscCall(TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 436 if (user.implicitform) { 437 PetscCall(TSSetIFunction(user.ts,NULL,IFunction,&user)); 438 PetscCall(TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user)); 439 PetscCall(TSSetType(user.ts,TSCN)); 440 } else { 441 PetscCall(TSSetRHSFunction(user.ts,NULL,RHSFunction,&user)); 442 PetscCall(TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user)); 443 PetscCall(TSSetType(user.ts,TSRK)); 444 } 445 PetscCall(TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user)); 446 PetscCall(TSSetMaxTime(user.ts,user.ftime)); 447 PetscCall(TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP)); 448 449 if (monitor) { 450 PetscCall(TSMonitorSet(user.ts,Monitor,&user,NULL)); 451 } 452 453 /* Set ODE initial conditions */ 454 PetscCall(VecGetArray(user.U,&x_ptr)); 455 x_ptr[0] = 2.0; 456 x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user.mu) - 292.0/(2187.0*user.mu*user.mu); 457 PetscCall(VecRestoreArray(user.U,&x_ptr)); 458 PetscCall(TSSetTimeStep(user.ts,PetscRealConstant(0.001))); 459 460 /* Set runtime options */ 461 PetscCall(TSSetFromOptions(user.ts)); 462 463 PetscCall(TSSolve(user.ts,user.U)); 464 PetscCall(VecGetArrayRead(user.U,&y_ptr)); 465 user.ob[0] = y_ptr[0]; 466 user.ob[1] = y_ptr[1]; 467 PetscCall(VecRestoreArrayRead(user.U,&y_ptr)); 468 469 /* Save trajectory of solution so that TSAdjointSolve() may be used. 470 Skip checkpointing for the first TSSolve since no adjoint run follows it. 471 */ 472 PetscCall(TSSetSaveTrajectory(user.ts)); 473 474 /* Optimization starts */ 475 PetscCall(MatCreate(PETSC_COMM_WORLD,&user.H)); 476 PetscCall(MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1)); 477 PetscCall(MatSetUp(user.H)); /* Hessian should be symmetric. Do we need to do MatSetOption(user.H,MAT_SYMMETRIC,PETSC_TRUE) ? */ 478 479 /* Set initial solution guess */ 480 PetscCall(MatCreateVecs(user.Jacp,&P,NULL)); 481 PetscCall(VecGetArray(P,&x_ptr)); 482 x_ptr[0] = PetscRealConstant(1.2); 483 PetscCall(VecRestoreArray(P,&x_ptr)); 484 PetscCall(TaoSetSolution(tao,P)); 485 486 /* Set routine for function and gradient evaluation */ 487 PetscCall(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 488 PetscCall(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 489 490 /* Check for any TAO command line options */ 491 PetscCall(TaoGetKSP(tao,&ksp)); 492 if (ksp) { 493 PetscCall(KSPGetPC(ksp,&pc)); 494 PetscCall(PCSetType(pc,PCNONE)); 495 } 496 PetscCall(TaoSetFromOptions(tao)); 497 498 PetscCall(TaoSolve(tao)); 499 500 PetscCall(VecView(P,PETSC_VIEWER_STDOUT_WORLD)); 501 /* Free TAO data structures */ 502 PetscCall(TaoDestroy(&tao)); 503 504 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 505 Free work space. All PETSc objects should be destroyed when they 506 are no longer needed. 507 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 508 PetscCall(MatDestroy(&user.H)); 509 PetscCall(MatDestroy(&user.A)); 510 PetscCall(VecDestroy(&user.U)); 511 PetscCall(MatDestroy(&user.Jacp)); 512 PetscCall(VecDestroy(&user.Lambda[0])); 513 PetscCall(VecDestroy(&user.Mup[0])); 514 PetscCall(VecDestroy(&user.Lambda2[0])); 515 PetscCall(VecDestroy(&user.Mup2[0])); 516 PetscCall(VecDestroy(&user.Ihp1[0])); 517 PetscCall(VecDestroy(&user.Ihp2[0])); 518 PetscCall(VecDestroy(&user.Ihp3[0])); 519 PetscCall(VecDestroy(&user.Ihp4[0])); 520 PetscCall(VecDestroy(&user.Dir)); 521 PetscCall(TSDestroy(&user.ts)); 522 PetscCall(VecDestroy(&P)); 523 PetscCall(PetscFinalize()); 524 return 0; 525 } 526 527 /* ------------------------------------------------------------------ */ 528 /* 529 FormFunctionGradient - Evaluates the function and corresponding gradient. 530 531 Input Parameters: 532 tao - the Tao context 533 X - the input vector 534 ptr - optional user-defined context, as set by TaoSetObjectiveAndGradient() 535 536 Output Parameters: 537 f - the newly evaluated function 538 G - the newly evaluated gradient 539 */ 540 PetscErrorCode FormFunctionGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx) 541 { 542 User user_ptr = (User)ctx; 543 TS ts = user_ptr->ts; 544 PetscScalar *x_ptr,*g; 545 const PetscScalar *y_ptr; 546 547 PetscFunctionBeginUser; 548 PetscCall(VecGetArrayRead(P,&y_ptr)); 549 user_ptr->mu = y_ptr[0]; 550 PetscCall(VecRestoreArrayRead(P,&y_ptr)); 551 552 PetscCall(TSSetTime(ts,0.0)); 553 PetscCall(TSSetStepNumber(ts,0)); 554 PetscCall(TSSetTimeStep(ts,PetscRealConstant(0.001))); /* can be overwritten by command line options */ 555 PetscCall(TSSetFromOptions(ts)); 556 PetscCall(VecGetArray(user_ptr->U,&x_ptr)); 557 x_ptr[0] = 2.0; 558 x_ptr[1] = -2.0/3.0 + 10.0/(81.0*user_ptr->mu) - 292.0/(2187.0*user_ptr->mu*user_ptr->mu); 559 PetscCall(VecRestoreArray(user_ptr->U,&x_ptr)); 560 561 PetscCall(TSSolve(ts,user_ptr->U)); 562 563 PetscCall(VecGetArrayRead(user_ptr->U,&y_ptr)); 564 *f = (y_ptr[0]-user_ptr->ob[0])*(y_ptr[0]-user_ptr->ob[0])+(y_ptr[1]-user_ptr->ob[1])*(y_ptr[1]-user_ptr->ob[1]); 565 566 /* Reset initial conditions for the adjoint integration */ 567 PetscCall(VecGetArray(user_ptr->Lambda[0],&x_ptr)); 568 x_ptr[0] = 2.*(y_ptr[0]-user_ptr->ob[0]); 569 x_ptr[1] = 2.*(y_ptr[1]-user_ptr->ob[1]); 570 PetscCall(VecRestoreArrayRead(user_ptr->U,&y_ptr)); 571 PetscCall(VecRestoreArray(user_ptr->Lambda[0],&x_ptr)); 572 573 PetscCall(VecGetArray(user_ptr->Mup[0],&x_ptr)); 574 x_ptr[0] = 0.0; 575 PetscCall(VecRestoreArray(user_ptr->Mup[0],&x_ptr)); 576 PetscCall(TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup)); 577 578 PetscCall(TSAdjointSolve(ts)); 579 580 PetscCall(VecGetArray(user_ptr->Mup[0],&x_ptr)); 581 PetscCall(VecGetArrayRead(user_ptr->Lambda[0],&y_ptr)); 582 PetscCall(VecGetArray(G,&g)); 583 g[0] = y_ptr[1]*(-10.0/(81.0*user_ptr->mu*user_ptr->mu)+2.0*292.0/(2187.0*user_ptr->mu*user_ptr->mu*user_ptr->mu))+x_ptr[0]; 584 PetscCall(VecRestoreArray(user_ptr->Mup[0],&x_ptr)); 585 PetscCall(VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr)); 586 PetscCall(VecRestoreArray(G,&g)); 587 PetscFunctionReturn(0); 588 } 589 590 PetscErrorCode FormHessian(Tao tao,Vec P,Mat H,Mat Hpre,void *ctx) 591 { 592 User user_ptr = (User)ctx; 593 PetscScalar harr[1]; 594 const PetscInt rows[1] = {0}; 595 PetscInt col = 0; 596 597 PetscFunctionBeginUser; 598 PetscCall(Adjoint2(P,harr,user_ptr)); 599 PetscCall(MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES)); 600 601 PetscCall(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 602 PetscCall(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 603 if (H != Hpre) { 604 PetscCall(MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY)); 605 PetscCall(MatAssemblyEnd(Hpre,MAT_FINAL_ASSEMBLY)); 606 } 607 PetscFunctionReturn(0); 608 } 609 610 PetscErrorCode Adjoint2(Vec P,PetscScalar arr[],User ctx) 611 { 612 TS ts = ctx->ts; 613 const PetscScalar *z_ptr; 614 PetscScalar *x_ptr,*y_ptr,dzdp,dzdp2; 615 Mat tlmsen; 616 617 PetscFunctionBeginUser; 618 /* Reset TSAdjoint so that AdjointSetUp will be called again */ 619 PetscCall(TSAdjointReset(ts)); 620 621 /* The directional vector should be 1 since it is one-dimensional */ 622 PetscCall(VecGetArray(ctx->Dir,&x_ptr)); 623 x_ptr[0] = 1.; 624 PetscCall(VecRestoreArray(ctx->Dir,&x_ptr)); 625 626 PetscCall(VecGetArrayRead(P,&z_ptr)); 627 ctx->mu = z_ptr[0]; 628 PetscCall(VecRestoreArrayRead(P,&z_ptr)); 629 630 dzdp = -10.0/(81.0*ctx->mu*ctx->mu) + 2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu); 631 dzdp2 = 2.*10.0/(81.0*ctx->mu*ctx->mu*ctx->mu) - 3.0*2.0*292.0/(2187.0*ctx->mu*ctx->mu*ctx->mu*ctx->mu); 632 633 PetscCall(TSSetTime(ts,0.0)); 634 PetscCall(TSSetStepNumber(ts,0)); 635 PetscCall(TSSetTimeStep(ts,PetscRealConstant(0.001))); /* can be overwritten by command line options */ 636 PetscCall(TSSetFromOptions(ts)); 637 PetscCall(TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir)); 638 639 PetscCall(MatZeroEntries(ctx->Jacp)); 640 PetscCall(MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES)); 641 PetscCall(MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY)); 642 PetscCall(MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY)); 643 644 PetscCall(TSAdjointSetForward(ts,ctx->Jacp)); 645 PetscCall(VecGetArray(ctx->U,&y_ptr)); 646 y_ptr[0] = 2.0; 647 y_ptr[1] = -2.0/3.0 + 10.0/(81.0*ctx->mu) - 292.0/(2187.0*ctx->mu*ctx->mu); 648 PetscCall(VecRestoreArray(ctx->U,&y_ptr)); 649 PetscCall(TSSolve(ts,ctx->U)); 650 651 /* Set terminal conditions for first- and second-order adjonts */ 652 PetscCall(VecGetArrayRead(ctx->U,&z_ptr)); 653 PetscCall(VecGetArray(ctx->Lambda[0],&y_ptr)); 654 y_ptr[0] = 2.*(z_ptr[0]-ctx->ob[0]); 655 y_ptr[1] = 2.*(z_ptr[1]-ctx->ob[1]); 656 PetscCall(VecRestoreArray(ctx->Lambda[0],&y_ptr)); 657 PetscCall(VecRestoreArrayRead(ctx->U,&z_ptr)); 658 PetscCall(VecGetArray(ctx->Mup[0],&y_ptr)); 659 y_ptr[0] = 0.0; 660 PetscCall(VecRestoreArray(ctx->Mup[0],&y_ptr)); 661 PetscCall(TSForwardGetSensitivities(ts,NULL,&tlmsen)); 662 PetscCall(MatDenseGetColumn(tlmsen,0,&x_ptr)); 663 PetscCall(VecGetArray(ctx->Lambda2[0],&y_ptr)); 664 y_ptr[0] = 2.*x_ptr[0]; 665 y_ptr[1] = 2.*x_ptr[1]; 666 PetscCall(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 667 PetscCall(VecGetArray(ctx->Mup2[0],&y_ptr)); 668 y_ptr[0] = 0.0; 669 PetscCall(VecRestoreArray(ctx->Mup2[0],&y_ptr)); 670 PetscCall(MatDenseRestoreColumn(tlmsen,&x_ptr)); 671 PetscCall(TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup)); 672 if (ctx->implicitform) { 673 PetscCall(TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx)); 674 } else { 675 PetscCall(TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx)); 676 } 677 PetscCall(TSAdjointSolve(ts)); 678 679 PetscCall(VecGetArray(ctx->Lambda[0],&x_ptr)); 680 PetscCall(VecGetArray(ctx->Lambda2[0],&y_ptr)); 681 PetscCall(VecGetArrayRead(ctx->Mup2[0],&z_ptr)); 682 683 arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0]; 684 685 PetscCall(VecRestoreArray(ctx->Lambda2[0],&x_ptr)); 686 PetscCall(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 687 PetscCall(VecRestoreArrayRead(ctx->Mup2[0],&z_ptr)); 688 689 /* Disable second-order adjoint mode */ 690 PetscCall(TSAdjointReset(ts)); 691 PetscCall(TSAdjointResetForward(ts)); 692 PetscFunctionReturn(0); 693 } 694 695 /*TEST 696 build: 697 requires: !complex !single 698 test: 699 args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.1 -viewer_binary_skip_info -tao_monitor -tao_view 700 output_file: output/ex20opt_p_1.out 701 702 test: 703 suffix: 2 704 args: -implicitform 0 -ts_type rk -ts_adapt_type none -mu 10 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none 705 output_file: output/ex20opt_p_2.out 706 707 test: 708 suffix: 3 709 args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_view 710 output_file: output/ex20opt_p_3.out 711 712 test: 713 suffix: 4 714 args: -ts_type cn -ts_adapt_type none -mu 100 -ts_dt 0.01 -viewer_binary_skip_info -tao_monitor -tao_type bntr -tao_bnk_pc_type none 715 output_file: output/ex20opt_p_4.out 716 717 TEST*/ 718