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 /* Set values for independent variables and parameters */ 133 PetscCall(VecGetArray(X, &x)); 134 x_a[0].setValue(x[0]); 135 x_a[1].setValue(x[1]); 136 mu_a.setValue(user->mu); 137 PetscCall(VecRestoreArray(X, &x)); 138 139 /* Set seed matrix as 3x3 identity matrix */ 140 x_a[0].setADValue(0, 1.); 141 x_a[0].setADValue(1, 0.); 142 x_a[0].setADValue(2, 0.); 143 x_a[1].setADValue(0, 0.); 144 x_a[1].setADValue(1, 1.); 145 x_a[1].setADValue(2, 0.); 146 mu_a.setADValue(0, 0.); 147 mu_a.setADValue(1, 0.); 148 mu_a.setADValue(2, 1.); 149 150 /* Evaluate residual (on active variables) */ 151 PetscCall(EvaluateResidual(x_a, mu_a, f_a)); 152 153 /* Extract derivatives */ 154 PetscCall(PetscMalloc1(2, &J)); 155 J[0] = (PetscScalar *)f_a[0].getADValue(); 156 J[1] = (PetscScalar *)f_a[1].getADValue(); 157 158 /* Set matrix values */ 159 for (i = 0; i < user->adctx->m; i++) PetscCall(MatSetValues(A, 1, &i, 1, &j, &J[i][user->adctx->n], INSERT_VALUES)); 160 PetscCall(PetscFree(J)); 161 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 162 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 163 PetscFunctionReturn(PETSC_SUCCESS); 164 } 165 166 /* 167 Monitor timesteps and use interpolation to output at integer multiples of 0.1 168 */ 169 static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal t, Vec X, void *ctx) 170 { 171 const PetscScalar *x; 172 PetscReal tfinal, dt, tprev; 173 User user = (User)ctx; 174 175 PetscFunctionBeginUser; 176 PetscCall(TSGetTimeStep(ts, &dt)); 177 PetscCall(TSGetMaxTime(ts, &tfinal)); 178 PetscCall(TSGetPrevTime(ts, &tprev)); 179 PetscCall(VecGetArrayRead(X, &x)); 180 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]))); 181 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "t %.6f (tprev = %.6f) \n", (double)t, (double)tprev)); 182 PetscCall(VecRestoreArrayRead(X, &x)); 183 PetscFunctionReturn(PETSC_SUCCESS); 184 } 185 186 int main(int argc, char **argv) 187 { 188 TS ts; /* nonlinear solver */ 189 Vec x; /* solution, residual vectors */ 190 Mat A; /* Jacobian matrix */ 191 Mat Jacp; /* JacobianP matrix */ 192 PetscInt steps; 193 PetscReal ftime = 0.5; 194 PetscBool monitor = PETSC_FALSE; 195 PetscScalar *x_ptr; 196 PetscMPIInt size; 197 struct _n_User user; 198 AdolcCtx *adctx; 199 Vec lambda[2], mu[2]; 200 201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202 Initialize program 203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204 PetscFunctionBeginUser; 205 PetscCall(PetscInitialize(&argc, &argv, NULL, help)); 206 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 207 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only!"); 208 209 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 210 Set runtime options and create AdolcCtx 211 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 212 PetscCall(PetscNew(&adctx)); 213 user.mu = 1; 214 user.next_output = 0.0; 215 adctx->m = 2; 216 adctx->n = 2; 217 adctx->p = 2; 218 user.adctx = adctx; 219 adtl::setNumDir(adctx->n + 1); /* #indep. variables, plus parameters */ 220 221 PetscCall(PetscOptionsGetReal(NULL, NULL, "-mu", &user.mu, NULL)); 222 PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitor", &monitor, NULL)); 223 224 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 225 Create necessary matrix and vectors, solve same ODE on every process 226 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 227 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 228 PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, 2, 2)); 229 PetscCall(MatSetFromOptions(A)); 230 PetscCall(MatSetUp(A)); 231 PetscCall(MatCreateVecs(A, &x, NULL)); 232 233 PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp)); 234 PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1)); 235 PetscCall(MatSetFromOptions(Jacp)); 236 PetscCall(MatSetUp(Jacp)); 237 238 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239 Create timestepping solver context 240 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 241 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 242 PetscCall(TSSetType(ts, TSRK)); 243 PetscCall(TSSetRHSFunction(ts, NULL, RHSFunctionPassive, &user)); 244 245 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 246 Set initial conditions 247 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 248 PetscCall(VecGetArray(x, &x_ptr)); 249 x_ptr[0] = 2; 250 x_ptr[1] = 0.66666654321; 251 PetscCall(VecRestoreArray(x, &x_ptr)); 252 253 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 254 Set RHS Jacobian for the adjoint integration 255 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 256 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, &user)); 257 PetscCall(TSSetMaxTime(ts, ftime)); 258 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP)); 259 if (monitor) PetscCall(TSMonitorSet(ts, Monitor, &user, NULL)); 260 PetscCall(TSSetTimeStep(ts, .001)); 261 262 /* 263 Have the TS save its trajectory so that TSAdjointSolve() may be used 264 */ 265 PetscCall(TSSetSaveTrajectory(ts)); 266 267 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 268 Set runtime options 269 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 270 PetscCall(TSSetFromOptions(ts)); 271 272 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 273 Solve nonlinear system 274 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 275 PetscCall(TSSolve(ts, x)); 276 PetscCall(TSGetSolveTime(ts, &ftime)); 277 PetscCall(TSGetStepNumber(ts, &steps)); 278 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "mu %g, steps %" PetscInt_FMT ", ftime %g\n", (double)user.mu, steps, (double)ftime)); 279 PetscCall(VecView(x, PETSC_VIEWER_STDOUT_WORLD)); 280 281 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 282 Start the Adjoint model 283 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 284 PetscCall(MatCreateVecs(A, &lambda[0], NULL)); 285 PetscCall(MatCreateVecs(A, &lambda[1], NULL)); 286 /* Reset initial conditions for the adjoint integration */ 287 PetscCall(VecGetArray(lambda[0], &x_ptr)); 288 x_ptr[0] = 1.0; 289 x_ptr[1] = 0.0; 290 PetscCall(VecRestoreArray(lambda[0], &x_ptr)); 291 PetscCall(VecGetArray(lambda[1], &x_ptr)); 292 x_ptr[0] = 0.0; 293 x_ptr[1] = 1.0; 294 PetscCall(VecRestoreArray(lambda[1], &x_ptr)); 295 296 PetscCall(MatCreateVecs(Jacp, &mu[0], NULL)); 297 PetscCall(MatCreateVecs(Jacp, &mu[1], NULL)); 298 PetscCall(VecGetArray(mu[0], &x_ptr)); 299 x_ptr[0] = 0.0; 300 PetscCall(VecRestoreArray(mu[0], &x_ptr)); 301 PetscCall(VecGetArray(mu[1], &x_ptr)); 302 x_ptr[0] = 0.0; 303 PetscCall(VecRestoreArray(mu[1], &x_ptr)); 304 PetscCall(TSSetCostGradients(ts, 2, lambda, mu)); 305 306 /* Set RHS JacobianP */ 307 PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &user)); 308 309 PetscCall(TSAdjointSolve(ts)); 310 311 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD)); 312 PetscCall(VecView(lambda[1], PETSC_VIEWER_STDOUT_WORLD)); 313 PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD)); 314 PetscCall(VecView(mu[1], PETSC_VIEWER_STDOUT_WORLD)); 315 316 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 317 Free work space. All PETSc objects should be destroyed when they 318 are no longer needed. 319 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 320 PetscCall(MatDestroy(&A)); 321 PetscCall(MatDestroy(&Jacp)); 322 PetscCall(VecDestroy(&x)); 323 PetscCall(VecDestroy(&lambda[0])); 324 PetscCall(VecDestroy(&lambda[1])); 325 PetscCall(VecDestroy(&mu[0])); 326 PetscCall(VecDestroy(&mu[1])); 327 PetscCall(TSDestroy(&ts)); 328 PetscCall(PetscFree(adctx)); 329 PetscCall(PetscFinalize()); 330 return 0; 331 } 332 333 /*TEST 334 335 build: 336 requires: double !complex adolc 337 338 test: 339 suffix: 1 340 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor 341 output_file: output/ex16adj_tl_1.out 342 343 test: 344 suffix: 2 345 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -mu 5 346 output_file: output/ex16adj_tl_2.out 347 348 test: 349 suffix: 3 350 args: -ts_max_steps 10 -monitor 351 output_file: output/ex16adj_tl_3.out 352 353 test: 354 suffix: 4 355 args: -ts_max_steps 10 -monitor -mu 5 356 output_file: output/ex16adj_tl_4.out 357 358 TEST*/ 359