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 CHKERRQ(VecGetArrayRead(U,&u)); 58 CHKERRQ(VecGetArray(F,&f)); 59 f[0] = u[1]; 60 f[1] = user->mu*((1.-u[0]*u[0])*u[1]-u[0]); 61 CHKERRQ(VecRestoreArrayRead(U,&u)); 62 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(A,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 82 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 83 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 84 if (B && A != B) { 85 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 86 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 87 } 88 89 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 103 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 104 CHKERRQ(VecGetArrayRead(Vr,&vr)); 105 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 118 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 119 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 120 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 133 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 134 CHKERRQ(VecGetArrayRead(Vr,&vr)); 135 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 147 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 148 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 149 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 162 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 163 CHKERRQ(VecGetArrayRead(Vr,&vr)); 164 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 176 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 177 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 178 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 198 CHKERRQ(VecGetArrayRead(Udot,&udot)); 199 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 205 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 206 CHKERRQ(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 CHKERRQ(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 CHKERRQ(MatSetValues(B,2,rowcol,2,rowcol,&J[0][0],INSERT_VALUES)); 223 CHKERRQ(VecRestoreArrayRead(U,&u)); 224 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 225 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 226 if (A != B) { 227 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 228 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 229 } 230 231 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 243 244 J[0][0] = 0; 245 J[1][0] = (1.-u[0]*u[0])*u[1]-u[0]; 246 CHKERRQ(MatSetValues(A,2,row,1,col,&J[0][0],INSERT_VALUES)); 247 CHKERRQ(VecRestoreArrayRead(U,&u)); 248 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 249 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 250 251 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 265 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 266 CHKERRQ(VecGetArrayRead(Vr,&vr)); 267 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 280 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 281 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 282 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 295 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 296 CHKERRQ(VecGetArrayRead(Vr,&vr)); 297 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 309 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 310 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 311 CHKERRQ(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 CHKERRQ(VecGetArrayRead(U,&u)); 324 CHKERRQ(VecGetArrayRead(Vl[0],&vl)); 325 CHKERRQ(VecGetArrayRead(Vr,&vr)); 326 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(U,&u)); 338 CHKERRQ(VecRestoreArrayRead(Vl[0],&vl)); 339 CHKERRQ(VecRestoreArrayRead(Vr,&vr)); 340 CHKERRQ(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 CHKERRQ(TSGetTimeStep(ts,&dt)); 361 CHKERRQ(TSGetMaxTime(ts,&tfinal)); 362 363 while (user->next_output <= t && user->next_output <= tfinal) { 364 CHKERRQ(VecDuplicate(X,&interpolatedX)); 365 CHKERRQ(TSInterpolate(ts,user->next_output,interpolatedX)); 366 CHKERRQ(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]));CHKERRQ(ierr); 370 CHKERRQ(VecRestoreArrayRead(interpolatedX,&x)); 371 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&argv,NULL,help)); 393 CHKERRMPI(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 CHKERRQ(TaoCreate(PETSC_COMM_WORLD,&tao)); 398 CHKERRQ(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 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-monitor",&monitor,NULL)); 409 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-mu",&user.mu,NULL)); 410 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-implicitform",&user.implicitform,NULL)); 411 412 /* Create necessary matrix and vectors, solve same ODE on every process */ 413 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.A)); 414 CHKERRQ(MatSetSizes(user.A,PETSC_DECIDE,PETSC_DECIDE,2,2)); 415 CHKERRQ(MatSetFromOptions(user.A)); 416 CHKERRQ(MatSetUp(user.A)); 417 CHKERRQ(MatCreateVecs(user.A,&user.U,NULL)); 418 CHKERRQ(MatCreateVecs(user.A,&user.Lambda[0],NULL)); 419 CHKERRQ(MatCreateVecs(user.A,&user.Lambda2[0],NULL)); 420 CHKERRQ(MatCreateVecs(user.A,&user.Ihp1[0],NULL)); 421 CHKERRQ(MatCreateVecs(user.A,&user.Ihp2[0],NULL)); 422 423 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.Jacp)); 424 CHKERRQ(MatSetSizes(user.Jacp,PETSC_DECIDE,PETSC_DECIDE,2,1)); 425 CHKERRQ(MatSetFromOptions(user.Jacp)); 426 CHKERRQ(MatSetUp(user.Jacp)); 427 CHKERRQ(MatCreateVecs(user.Jacp,&user.Dir,NULL)); 428 CHKERRQ(MatCreateVecs(user.Jacp,&user.Mup[0],NULL)); 429 CHKERRQ(MatCreateVecs(user.Jacp,&user.Mup2[0],NULL)); 430 CHKERRQ(MatCreateVecs(user.Jacp,&user.Ihp3[0],NULL)); 431 CHKERRQ(MatCreateVecs(user.Jacp,&user.Ihp4[0],NULL)); 432 433 /* Create timestepping solver context */ 434 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&user.ts)); 435 CHKERRQ(TSSetEquationType(user.ts,TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 436 if (user.implicitform) { 437 CHKERRQ(TSSetIFunction(user.ts,NULL,IFunction,&user)); 438 CHKERRQ(TSSetIJacobian(user.ts,user.A,user.A,IJacobian,&user)); 439 CHKERRQ(TSSetType(user.ts,TSCN)); 440 } else { 441 CHKERRQ(TSSetRHSFunction(user.ts,NULL,RHSFunction,&user)); 442 CHKERRQ(TSSetRHSJacobian(user.ts,user.A,user.A,RHSJacobian,&user)); 443 CHKERRQ(TSSetType(user.ts,TSRK)); 444 } 445 CHKERRQ(TSSetRHSJacobianP(user.ts,user.Jacp,RHSJacobianP,&user)); 446 CHKERRQ(TSSetMaxTime(user.ts,user.ftime)); 447 CHKERRQ(TSSetExactFinalTime(user.ts,TS_EXACTFINALTIME_MATCHSTEP)); 448 449 if (monitor) { 450 CHKERRQ(TSMonitorSet(user.ts,Monitor,&user,NULL)); 451 } 452 453 /* Set ODE initial conditions */ 454 CHKERRQ(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 CHKERRQ(VecRestoreArray(user.U,&x_ptr)); 458 CHKERRQ(TSSetTimeStep(user.ts,PetscRealConstant(0.001))); 459 460 /* Set runtime options */ 461 CHKERRQ(TSSetFromOptions(user.ts)); 462 463 CHKERRQ(TSSolve(user.ts,user.U)); 464 CHKERRQ(VecGetArrayRead(user.U,&y_ptr)); 465 user.ob[0] = y_ptr[0]; 466 user.ob[1] = y_ptr[1]; 467 CHKERRQ(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 CHKERRQ(TSSetSaveTrajectory(user.ts)); 473 474 /* Optimization starts */ 475 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&user.H)); 476 CHKERRQ(MatSetSizes(user.H,PETSC_DECIDE,PETSC_DECIDE,1,1)); 477 CHKERRQ(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 CHKERRQ(MatCreateVecs(user.Jacp,&P,NULL)); 481 CHKERRQ(VecGetArray(P,&x_ptr)); 482 x_ptr[0] = PetscRealConstant(1.2); 483 CHKERRQ(VecRestoreArray(P,&x_ptr)); 484 CHKERRQ(TaoSetSolution(tao,P)); 485 486 /* Set routine for function and gradient evaluation */ 487 CHKERRQ(TaoSetObjectiveAndGradient(tao,NULL,FormFunctionGradient,(void *)&user)); 488 CHKERRQ(TaoSetHessian(tao,user.H,user.H,FormHessian,(void *)&user)); 489 490 /* Check for any TAO command line options */ 491 CHKERRQ(TaoGetKSP(tao,&ksp)); 492 if (ksp) { 493 CHKERRQ(KSPGetPC(ksp,&pc)); 494 CHKERRQ(PCSetType(pc,PCNONE)); 495 } 496 CHKERRQ(TaoSetFromOptions(tao)); 497 498 CHKERRQ(TaoSolve(tao)); 499 500 CHKERRQ(VecView(P,PETSC_VIEWER_STDOUT_WORLD)); 501 /* Free TAO data structures */ 502 CHKERRQ(TaoDestroy(&tao)); 503 504 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 505 Free work space. All PETSc objects should be destroyed when they 506 are no longer needed. 507 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 508 CHKERRQ(MatDestroy(&user.H)); 509 CHKERRQ(MatDestroy(&user.A)); 510 CHKERRQ(VecDestroy(&user.U)); 511 CHKERRQ(MatDestroy(&user.Jacp)); 512 CHKERRQ(VecDestroy(&user.Lambda[0])); 513 CHKERRQ(VecDestroy(&user.Mup[0])); 514 CHKERRQ(VecDestroy(&user.Lambda2[0])); 515 CHKERRQ(VecDestroy(&user.Mup2[0])); 516 CHKERRQ(VecDestroy(&user.Ihp1[0])); 517 CHKERRQ(VecDestroy(&user.Ihp2[0])); 518 CHKERRQ(VecDestroy(&user.Ihp3[0])); 519 CHKERRQ(VecDestroy(&user.Ihp4[0])); 520 CHKERRQ(VecDestroy(&user.Dir)); 521 CHKERRQ(TSDestroy(&user.ts)); 522 CHKERRQ(VecDestroy(&P)); 523 CHKERRQ(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 CHKERRQ(VecGetArrayRead(P,&y_ptr)); 549 user_ptr->mu = y_ptr[0]; 550 CHKERRQ(VecRestoreArrayRead(P,&y_ptr)); 551 552 CHKERRQ(TSSetTime(ts,0.0)); 553 CHKERRQ(TSSetStepNumber(ts,0)); 554 CHKERRQ(TSSetTimeStep(ts,PetscRealConstant(0.001))); /* can be overwritten by command line options */ 555 CHKERRQ(TSSetFromOptions(ts)); 556 CHKERRQ(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 CHKERRQ(VecRestoreArray(user_ptr->U,&x_ptr)); 560 561 CHKERRQ(TSSolve(ts,user_ptr->U)); 562 563 CHKERRQ(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 CHKERRQ(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 CHKERRQ(VecRestoreArrayRead(user_ptr->U,&y_ptr)); 571 CHKERRQ(VecRestoreArray(user_ptr->Lambda[0],&x_ptr)); 572 573 CHKERRQ(VecGetArray(user_ptr->Mup[0],&x_ptr)); 574 x_ptr[0] = 0.0; 575 CHKERRQ(VecRestoreArray(user_ptr->Mup[0],&x_ptr)); 576 CHKERRQ(TSSetCostGradients(ts,1,user_ptr->Lambda,user_ptr->Mup)); 577 578 CHKERRQ(TSAdjointSolve(ts)); 579 580 CHKERRQ(VecGetArray(user_ptr->Mup[0],&x_ptr)); 581 CHKERRQ(VecGetArrayRead(user_ptr->Lambda[0],&y_ptr)); 582 CHKERRQ(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 CHKERRQ(VecRestoreArray(user_ptr->Mup[0],&x_ptr)); 585 CHKERRQ(VecRestoreArrayRead(user_ptr->Lambda[0],&y_ptr)); 586 CHKERRQ(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 CHKERRQ(Adjoint2(P,harr,user_ptr)); 599 CHKERRQ(MatSetValues(H,1,rows,1,&col,harr,INSERT_VALUES)); 600 601 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 602 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 603 if (H != Hpre) { 604 CHKERRQ(MatAssemblyBegin(Hpre,MAT_FINAL_ASSEMBLY)); 605 CHKERRQ(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 CHKERRQ(TSAdjointReset(ts)); 620 621 /* The directional vector should be 1 since it is one-dimensional */ 622 CHKERRQ(VecGetArray(ctx->Dir,&x_ptr)); 623 x_ptr[0] = 1.; 624 CHKERRQ(VecRestoreArray(ctx->Dir,&x_ptr)); 625 626 CHKERRQ(VecGetArrayRead(P,&z_ptr)); 627 ctx->mu = z_ptr[0]; 628 CHKERRQ(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 CHKERRQ(TSSetTime(ts,0.0)); 634 CHKERRQ(TSSetStepNumber(ts,0)); 635 CHKERRQ(TSSetTimeStep(ts,PetscRealConstant(0.001))); /* can be overwritten by command line options */ 636 CHKERRQ(TSSetFromOptions(ts)); 637 CHKERRQ(TSSetCostHessianProducts(ts,1,ctx->Lambda2,ctx->Mup2,ctx->Dir)); 638 639 CHKERRQ(MatZeroEntries(ctx->Jacp)); 640 CHKERRQ(MatSetValue(ctx->Jacp,1,0,dzdp,INSERT_VALUES)); 641 CHKERRQ(MatAssemblyBegin(ctx->Jacp,MAT_FINAL_ASSEMBLY)); 642 CHKERRQ(MatAssemblyEnd(ctx->Jacp,MAT_FINAL_ASSEMBLY)); 643 644 CHKERRQ(TSAdjointSetForward(ts,ctx->Jacp)); 645 CHKERRQ(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 CHKERRQ(VecRestoreArray(ctx->U,&y_ptr)); 649 CHKERRQ(TSSolve(ts,ctx->U)); 650 651 /* Set terminal conditions for first- and second-order adjonts */ 652 CHKERRQ(VecGetArrayRead(ctx->U,&z_ptr)); 653 CHKERRQ(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 CHKERRQ(VecRestoreArray(ctx->Lambda[0],&y_ptr)); 657 CHKERRQ(VecRestoreArrayRead(ctx->U,&z_ptr)); 658 CHKERRQ(VecGetArray(ctx->Mup[0],&y_ptr)); 659 y_ptr[0] = 0.0; 660 CHKERRQ(VecRestoreArray(ctx->Mup[0],&y_ptr)); 661 CHKERRQ(TSForwardGetSensitivities(ts,NULL,&tlmsen)); 662 CHKERRQ(MatDenseGetColumn(tlmsen,0,&x_ptr)); 663 CHKERRQ(VecGetArray(ctx->Lambda2[0],&y_ptr)); 664 y_ptr[0] = 2.*x_ptr[0]; 665 y_ptr[1] = 2.*x_ptr[1]; 666 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 667 CHKERRQ(VecGetArray(ctx->Mup2[0],&y_ptr)); 668 y_ptr[0] = 0.0; 669 CHKERRQ(VecRestoreArray(ctx->Mup2[0],&y_ptr)); 670 CHKERRQ(MatDenseRestoreColumn(tlmsen,&x_ptr)); 671 CHKERRQ(TSSetCostGradients(ts,1,ctx->Lambda,ctx->Mup)); 672 if (ctx->implicitform) { 673 CHKERRQ(TSSetIHessianProduct(ts,ctx->Ihp1,IHessianProductUU,ctx->Ihp2,IHessianProductUP,ctx->Ihp3,IHessianProductPU,ctx->Ihp4,IHessianProductPP,ctx)); 674 } else { 675 CHKERRQ(TSSetRHSHessianProduct(ts,ctx->Ihp1,RHSHessianProductUU,ctx->Ihp2,RHSHessianProductUP,ctx->Ihp3,RHSHessianProductPU,ctx->Ihp4,RHSHessianProductPP,ctx)); 676 } 677 CHKERRQ(TSAdjointSolve(ts)); 678 679 CHKERRQ(VecGetArray(ctx->Lambda[0],&x_ptr)); 680 CHKERRQ(VecGetArray(ctx->Lambda2[0],&y_ptr)); 681 CHKERRQ(VecGetArrayRead(ctx->Mup2[0],&z_ptr)); 682 683 arr[0] = x_ptr[1]*dzdp2 + y_ptr[1]*dzdp2 + z_ptr[0]; 684 685 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&x_ptr)); 686 CHKERRQ(VecRestoreArray(ctx->Lambda2[0],&y_ptr)); 687 CHKERRQ(VecRestoreArrayRead(ctx->Mup2[0],&z_ptr)); 688 689 /* Disable second-order adjoint mode */ 690 CHKERRQ(TSAdjointReset(ts)); 691 CHKERRQ(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