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