1 static char help[] = "Demonstrates tapeless automatic Jacobian generation using ADOL-C for an adjoint sensitivity analysis of the van der Pol equation.\n\ 2 Input parameters include:\n\ 3 -mu : stiffness parameter\n\n"; 4 5 /* 6 REQUIRES configuration of PETSc with option --download-adolc. 7 8 For documentation on ADOL-C, see 9 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 10 */ 11 /* ------------------------------------------------------------------------ 12 See ex16adj for a description of the problem being solved. 13 ------------------------------------------------------------------------- */ 14 15 #include <petscts.h> 16 #include <petscmat.h> 17 18 #define ADOLC_TAPELESS 19 #define NUMBER_DIRECTIONS 3 20 #include "adolc-utils/drivers.cxx" 21 #include <adolc/adtl.h> 22 using namespace adtl; 23 24 typedef struct _n_User *User; 25 struct _n_User { 26 PetscReal mu; 27 PetscReal next_output; 28 PetscReal tprev; 29 30 /* Automatic differentiation support */ 31 AdolcCtx *adctx; 32 Vec F; 33 }; 34 35 /* 36 Residual evaluation templated, so as to allow for PetscScalar or adouble 37 arguments. 38 */ 39 template <class T> 40 PetscErrorCode EvaluateResidual(const T *x, T mu, T *f) 41 { 42 PetscFunctionBegin; 43 f[0] = x[1]; 44 f[1] = mu * (1. - x[0] * x[0]) * x[1] - x[0]; 45 PetscFunctionReturn(PETSC_SUCCESS); 46 } 47 48 /* 49 'Passive' RHS function, used in residual evaluations during the time integration. 50 */ 51 static PetscErrorCode RHSFunctionPassive(TS ts, PetscReal t, Vec X, Vec F, void *ctx) 52 { 53 User user = (User)ctx; 54 PetscScalar *f; 55 const PetscScalar *x; 56 57 PetscFunctionBeginUser; 58 PetscCall(VecGetArrayRead(X, &x)); 59 PetscCall(VecGetArray(F, &f)); 60 PetscCall(EvaluateResidual(x, user->mu, f)); 61 PetscCall(VecRestoreArrayRead(X, &x)); 62 PetscCall(VecRestoreArray(F, &f)); 63 PetscFunctionReturn(PETSC_SUCCESS); 64 } 65 66 /* 67 Compute the Jacobian w.r.t. x using tapeless mode of ADOL-C. 68 */ 69 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat A, Mat B, void *ctx) 70 { 71 User user = (User)ctx; 72 PetscScalar **J; 73 const PetscScalar *x; 74 adouble f_a[2]; /* 'active' double for dependent variables */ 75 adouble x_a[2], mu_a; /* 'active' doubles for independent variables */ 76 PetscInt i, j; 77 78 PetscFunctionBeginUser; 79 /* Set values for independent variables and parameters */ 80 PetscCall(VecGetArrayRead(X, &x)); 81 x_a[0].setValue(x[0]); 82 x_a[1].setValue(x[1]); 83 mu_a.setValue(user->mu); 84 PetscCall(VecRestoreArrayRead(X, &x)); 85 86 /* Set seed matrix as 3x3 identity matrix */ 87 x_a[0].setADValue(0, 1.); 88 x_a[0].setADValue(1, 0.); 89 x_a[0].setADValue(2, 0.); 90 x_a[1].setADValue(0, 0.); 91 x_a[1].setADValue(1, 1.); 92 x_a[1].setADValue(2, 0.); 93 mu_a.setADValue(0, 0.); 94 mu_a.setADValue(1, 0.); 95 mu_a.setADValue(2, 1.); 96 97 /* Evaluate residual (on active variables) */ 98 PetscCall(EvaluateResidual(x_a, mu_a, f_a)); 99 100 /* Extract derivatives */ 101 PetscCall(PetscMalloc1(user->adctx->n, &J)); 102 J[0] = (PetscScalar *)f_a[0].getADValue(); 103 J[1] = (PetscScalar *)f_a[1].getADValue(); 104 105 /* Set matrix values */ 106 for (i = 0; i < user->adctx->m; i++) { 107 for (j = 0; j < user->adctx->n; j++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][j], INSERT_VALUES)); 108 } 109 PetscCall(PetscFree(J)); 110 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 111 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 112 if (A != B) { 113 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 114 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 115 } 116 PetscFunctionReturn(PETSC_SUCCESS); 117 } 118 119 /* 120 Compute the Jacobian w.r.t. mu using tapeless mode of ADOL-C. 121 */ 122 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, void *ctx) 123 { 124 User user = (User)ctx; 125 PetscScalar **J; 126 PetscScalar *x; 127 adouble f_a[2]; /* 'active' double for dependent variables */ 128 adouble x_a[2], mu_a; /* 'active' doubles for independent variables */ 129 PetscInt i, j = 0; 130 131 PetscFunctionBeginUser; 132 133 /* Set values for independent variables and parameters */ 134 PetscCall(VecGetArray(X, &x)); 135 x_a[0].setValue(x[0]); 136 x_a[1].setValue(x[1]); 137 mu_a.setValue(user->mu); 138 PetscCall(VecRestoreArray(X, &x)); 139 140 /* Set seed matrix as 3x3 identity matrix */ 141 x_a[0].setADValue(0, 1.); 142 x_a[0].setADValue(1, 0.); 143 x_a[0].setADValue(2, 0.); 144 x_a[1].setADValue(0, 0.); 145 x_a[1].setADValue(1, 1.); 146 x_a[1].setADValue(2, 0.); 147 mu_a.setADValue(0, 0.); 148 mu_a.setADValue(1, 0.); 149 mu_a.setADValue(2, 1.); 150 151 /* Evaluate residual (on active variables) */ 152 PetscCall(EvaluateResidual(x_a, mu_a, f_a)); 153 154 /* Extract derivatives */ 155 PetscCall(PetscMalloc1(2, &J)); 156 J[0] = (PetscScalar *)f_a[0].getADValue(); 157 J[1] = (PetscScalar *)f_a[1].getADValue(); 158 159 /* Set matrix values */ 160 for (i = 0; i < user->adctx->m; i++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][user->adctx->n], INSERT_VALUES)); 161 PetscCall(PetscFree(J)); 162 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 163 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 164 PetscFunctionReturn(PETSC_SUCCESS); 165 } 166 167 /* 168 Monitor timesteps and use interpolation to output at integer multiples of 0.1 169 */ 170 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx) 171 { 172 const PetscScalar *x; 173 PetscReal tfinal, dt, tprev; 174 User user = (User)ctx; 175 176 PetscFunctionBeginUser; 177 PetscCall(TSGetTimeStep(ts, &dt)); 178 PetscCall(TSGetMaxTime(ts, &tfinal)); 179 PetscCall(TSGetPrevTime(ts, &tprev)); 180 PetscCall(VecGetArrayRead(X, &x)); 181 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[%.1f] %" PetscInt_FMT " TS %.6f (dt = %.6f) X % 12.6e % 12.6e\n", (double)user->next_output, step, (double)t, (double)dt, (double)PetscRealPart(x[0]), (double)PetscRealPart(x[1]))); 182 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "t %.6f (tprev = %.6f) \n", (double)t, (double)tprev)); 183 PetscCall(VecRestoreArrayRead(X, &x)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 int main(int argc, char **argv) 188 { 189 TS ts; /* nonlinear solver */ 190 Vec x; /* solution, residual vectors */ 191 Mat A; /* Jacobian matrix */ 192 Mat Jacp; /* JacobianP matrix */ 193 PetscInt steps; 194 PetscReal ftime = 0.5; 195 PetscBool monitor = PETSC_FALSE; 196 PetscScalar *x_ptr; 197 PetscMPIInt size; 198 struct _n_User user; 199 AdolcCtx *adctx; 200 Vec lambda[2], mu[2]; 201 202 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 203 Initialize program 204 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 205 PetscFunctionBeginUser; 206 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 207 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 208 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 209 210 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 211 Set runtime options and create AdolcCtx 212 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 213 PetscCall(PetscNew(&adctx)); 214 user.mu = 1; 215 user.next_output = 0.0; 216 adctx->m = 2; 217 adctx->n = 2; 218 adctx->p = 2; 219 user.adctx = adctx; 220 adtl::setNumDir(adctx->n + 1); /* #indep. variables, plus parameters */ 221 222 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL)); 223 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL)); 224 225 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 226 Create necessary matrix and vectors, solve same ODE on every process 227 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 228 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 229 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 230 PetscCall(MatSetFromOptions(A)); 231 PetscCall(MatSetUp(A)); 232 PetscCall(MatCreateVecs(A, &x, NULL)); 233 234 PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp)); 235 PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1)); 236 PetscCall(MatSetFromOptions(Jacp)); 237 PetscCall(MatSetUp(Jacp)); 238 239 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 240 Create timestepping solver context 241 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 242 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 243 PetscCall(TSSetType(ts, TSRK)); 244 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionPassive, &user)); 245 246 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 247 Set initial conditions 248 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 249 PetscCall(VecGetArray(x, &x_ptr)); 250 x_ptr[0] = 2; 251 x_ptr[1] = 0.66666654321; 252 PetscCall(VecRestoreArray(x, &x_ptr)); 253 254 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 255 Set RHS Jacobian for the adjoint integration 256 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 257 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, &user)); 258 PetscCall(TSSetMaxTime(ts, ftime)); 259 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 260 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL)); 261 PetscCall(TSSetTimeStep(ts, .001)); 262 263 /* 264 Have the TS save its trajectory so that TSAdjointSolve() may be used 265 */ 266 PetscCall(TSSetSaveTrajectory(ts)); 267 268 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 269 Set runtime options 270 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 271 PetscCall(TSSetFromOptions(ts)); 272 273 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 274 Solve nonlinear system 275 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 276 PetscCall(TSSolve(ts, x)); 277 PetscCall(TSGetSolveTime(ts, &ftime)); 278 PetscCall(TSGetStepNumber(ts, &steps)); 279 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, steps, (double)ftime)); 280 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 281 282 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 283 Start the Adjoint model 284 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 285 PetscCall(MatCreateVecs(A, &lambda[0], NULL)); 286 PetscCall(MatCreateVecs(A, &lambda[1], NULL)); 287 /* Reset initial conditions for the adjoint integration */ 288 PetscCall(VecGetArray(lambda[0], &x_ptr)); 289 x_ptr[0] = 1.0; 290 x_ptr[1] = 0.0; 291 PetscCall(VecRestoreArray(lambda[0], &x_ptr)); 292 PetscCall(VecGetArray(lambda[1], &x_ptr)); 293 x_ptr[0] = 0.0; 294 x_ptr[1] = 1.0; 295 PetscCall(VecRestoreArray(lambda[1], &x_ptr)); 296 297 PetscCall(MatCreateVecs(Jacp, &mu[0], NULL)); 298 PetscCall(MatCreateVecs(Jacp, &mu[1], NULL)); 299 PetscCall(VecGetArray(mu[0], &x_ptr)); 300 x_ptr[0] = 0.0; 301 PetscCall(VecRestoreArray(mu[0], &x_ptr)); 302 PetscCall(VecGetArray(mu[1], &x_ptr)); 303 x_ptr[0] = 0.0; 304 PetscCall(VecRestoreArray(mu[1], &x_ptr)); 305 PetscCall(TSSetCostGradients(ts, 2, lambda, mu)); 306 307 /* Set RHS JacobianP */ 308 PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &user)); 309 310 PetscCall(TSAdjointSolve(ts)); 311 312 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD)); 313 PetscCall(VecView(lambda[1], PETSC_VIEWER_STDOUT_WORLD)); 314 PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD)); 315 PetscCall(VecView(mu[1], PETSC_VIEWER_STDOUT_WORLD)); 316 317 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 318 Free work space. All PETSc objects should be destroyed when they 319 are no longer needed. 320 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 321 PetscCall(MatDestroy(&A)); 322 PetscCall(MatDestroy(&Jacp)); 323 PetscCall(VecDestroy(&x)); 324 PetscCall(VecDestroy(&lambda[0])); 325 PetscCall(VecDestroy(&lambda[1])); 326 PetscCall(VecDestroy(&mu[0])); 327 PetscCall(VecDestroy(&mu[1])); 328 PetscCall(TSDestroy(&ts)); 329 PetscCall(PetscFree(adctx)); 330 PetscCall(PetscFinalize()); 331 return 0; 332 } 333 334 /*TEST 335 336 build: 337 requires: double !complex adolc 338 339 test: 340 suffix: 1 341 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor 342 output_file: output/ex16adj_tl_1.out 343 344 test: 345 suffix: 2 346 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -mu 5 347 output_file: output/ex16adj_tl_2.out 348 349 test: 350 suffix: 3 351 args: -ts_max_steps 10 -monitor 352 output_file: output/ex16adj_tl_3.out 353 354 test: 355 suffix: 4 356 args: -ts_max_steps 10 -monitor -mu 5 357 output_file: output/ex16adj_tl_4.out 358 359 TEST*/ 360