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