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