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