1 static char help[] = "Demonstrates automatic Jacobian generation using ADOL-C for a nonlinear reaction problem from chemistry.\n"; 2 3 /* 4 Concepts: TS^time-dependent nonlinear problems 5 Concepts: Automatic differentiation using ADOL-C 6 Processors: 1 7 */ 8 /* 9 REQUIRES configuration of PETSc with option --download-adolc. 10 11 For documentation on ADOL-C, see 12 $PETSC_ARCH/externalpackages/ADOL-C-2.6.0/ADOL-C/doc/adolc-manual.pdf 13 */ 14 /* ------------------------------------------------------------------------ 15 See ../advection-diffusion-reaction/ex1 for a description of the problem 16 ------------------------------------------------------------------------- */ 17 #include <petscts.h> 18 #include "adolc-utils/drivers.cxx" 19 #include <adolc/adolc.h> 20 21 typedef struct { 22 PetscScalar k; 23 Vec initialsolution; 24 AdolcCtx *adctx; /* Automatic differentiation support */ 25 } AppCtx; 26 27 PetscErrorCode IFunctionView(AppCtx *ctx,PetscViewer v) 28 { 29 PetscErrorCode ierr; 30 31 PetscFunctionBegin; 32 ierr = PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 37 { 38 PetscErrorCode ierr; 39 40 PetscFunctionBegin; 41 ierr = PetscNew(ctx);CHKERRQ(ierr); 42 ierr = PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR);CHKERRQ(ierr); 43 PetscFunctionReturn(0); 44 } 45 46 /* 47 Defines the ODE passed to the ODE solver 48 */ 49 PetscErrorCode IFunctionPassive(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 50 { 51 PetscErrorCode ierr; 52 PetscScalar *f; 53 const PetscScalar *u,*udot; 54 55 PetscFunctionBegin; 56 /* The next three lines allow us to access the entries of the vectors directly */ 57 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 58 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 59 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 60 f[0] = udot[0] + ctx->k*u[0]*u[1]; 61 f[1] = udot[1] + ctx->k*u[0]*u[1]; 62 f[2] = udot[2] - ctx->k*u[0]*u[1]; 63 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 64 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 65 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 66 PetscFunctionReturn(0); 67 } 68 69 /* 70 'Active' ADOL-C annotated version, marking dependence upon u. 71 */ 72 PetscErrorCode IFunctionActive1(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 73 { 74 PetscErrorCode ierr; 75 PetscScalar *f; 76 const PetscScalar *u,*udot; 77 78 adouble f_a[3]; /* 'active' double for dependent variables */ 79 adouble u_a[3]; /* 'active' double for independent variables */ 80 81 PetscFunctionBegin; 82 /* The next three lines allow us to access the entries of the vectors directly */ 83 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 84 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 85 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 86 87 /* Start of active section */ 88 trace_on(1); 89 u_a[0] <<= u[0]; u_a[1] <<= u[1]; u_a[2] <<= u[2]; /* Mark independence */ 90 f_a[0] = udot[0] + ctx->k*u_a[0]*u_a[1]; 91 f_a[1] = udot[1] + ctx->k*u_a[0]*u_a[1]; 92 f_a[2] = udot[2] - ctx->k*u_a[0]*u_a[1]; 93 f_a[0] >>= f[0]; f_a[1] >>= f[1]; f_a[2] >>= f[2]; /* Mark dependence */ 94 trace_off(); 95 /* End of active section */ 96 97 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 98 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 99 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 100 PetscFunctionReturn(0); 101 } 102 103 /* 104 'Active' ADOL-C annotated version, marking dependence upon udot. 105 */ 106 PetscErrorCode IFunctionActive2(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 107 { 108 PetscErrorCode ierr; 109 PetscScalar *f; 110 const PetscScalar *u,*udot; 111 112 adouble f_a[3]; /* 'active' double for dependent variables */ 113 adouble udot_a[3]; /* 'active' double for independent variables */ 114 115 PetscFunctionBegin; 116 /* The next three lines allow us to access the entries of the vectors directly */ 117 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 118 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 119 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 120 121 /* Start of active section */ 122 trace_on(2); 123 udot_a[0] <<= udot[0]; udot_a[1] <<= udot[1]; udot_a[2] <<= udot[2]; /* Mark independence */ 124 f_a[0] = udot_a[0] + ctx->k*u[0]*u[1]; 125 f_a[1] = udot_a[1] + ctx->k*u[0]*u[1]; 126 f_a[2] = udot_a[2] - ctx->k*u[0]*u[1]; 127 f_a[0] >>= f[0]; f_a[1] >>= f[1]; f_a[2] >>= f[2]; /* Mark dependence */ 128 trace_off(); 129 /* End of active section */ 130 131 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 132 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 133 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 134 PetscFunctionReturn(0); 135 } 136 137 /* 138 Defines the Jacobian of the ODE passed to the ODE solver, using the PETSc-ADOL-C driver for 139 implicit TS. 140 */ 141 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 142 { 143 PetscErrorCode ierr; 144 AppCtx *appctx = (AppCtx*)ctx; 145 const PetscScalar *u; 146 147 PetscFunctionBegin; 148 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 149 ierr = PetscAdolcComputeIJacobian(1,2,A,u,a,appctx->adctx);CHKERRQ(ierr); 150 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 151 PetscFunctionReturn(0); 152 } 153 154 /* 155 Defines the exact (analytic) solution to the ODE 156 */ 157 static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 158 { 159 PetscErrorCode ierr; 160 const PetscScalar *uinit; 161 PetscScalar *u,d0,q; 162 163 PetscFunctionBegin; 164 ierr = VecGetArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 165 ierr = VecGetArray(U,&u);CHKERRQ(ierr); 166 d0 = uinit[0] - uinit[1]; 167 if (d0 == 0.0) q = ctx->k*t; 168 else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0; 169 u[0] = uinit[0]/(1.0 + uinit[1]*q); 170 u[1] = u[0] - d0; 171 u[2] = uinit[1] + uinit[2] - u[1]; 172 ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 173 ierr = VecRestoreArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 174 PetscFunctionReturn(0); 175 } 176 177 int main(int argc,char **argv) 178 { 179 TS ts; /* ODE integrator */ 180 Vec U,Udot,R; /* solution, derivative, residual */ 181 Mat A; /* Jacobian matrix */ 182 PetscErrorCode ierr; 183 PetscMPIInt size; 184 PetscInt n = 3; 185 AppCtx ctx; 186 AdolcCtx *adctx; 187 PetscScalar *u; 188 const char * const names[] = {"U1","U2","U3",NULL}; 189 190 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191 Initialize program 192 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 193 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr; 194 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 195 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 196 ierr = PetscNew(&adctx);CHKERRQ(ierr); 197 adctx->m = n;adctx->n = n;adctx->p = n; 198 ctx.adctx = adctx; 199 200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 201 Create necessary matrix and vectors 202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 203 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 204 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 205 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 206 ierr = MatSetUp(A);CHKERRQ(ierr); 207 208 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 209 210 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 211 Set runtime options 212 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 213 ctx.k = .9; 214 ierr = PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL);CHKERRQ(ierr); 215 ierr = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr); 216 ierr = VecGetArray(ctx.initialsolution,&u);CHKERRQ(ierr); 217 u[0] = 1; 218 u[1] = .7; 219 u[2] = 0; 220 ierr = VecRestoreArray(ctx.initialsolution,&u);CHKERRQ(ierr); 221 ierr = PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL);CHKERRQ(ierr); 222 223 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 224 Create timestepping solver context 225 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 226 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 227 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 228 ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 229 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunctionPassive,&ctx);CHKERRQ(ierr); 230 231 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 232 Set initial conditions 233 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 234 ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr); 235 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 236 237 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 238 Trace just once for each tape 239 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 240 ierr = VecDuplicate(U,&Udot);CHKERRQ(ierr); 241 ierr = VecDuplicate(U,&R);CHKERRQ(ierr); 242 ierr = IFunctionActive1(ts,0.,U,Udot,R,&ctx);CHKERRQ(ierr); 243 ierr = IFunctionActive2(ts,0.,U,Udot,R,&ctx);CHKERRQ(ierr); 244 ierr = VecDestroy(&R);CHKERRQ(ierr); 245 ierr = VecDestroy(&Udot);CHKERRQ(ierr); 246 247 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 248 Set Jacobian 249 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 250 ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr); 251 ierr = TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx);CHKERRQ(ierr); 252 253 { 254 DM dm; 255 void *ptr; 256 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 257 ierr = PetscDLSym(NULL,"IFunctionView",&ptr);CHKERRQ(ierr); 258 ierr = PetscDLSym(NULL,"IFunctionLoad",&ptr);CHKERRQ(ierr); 259 ierr = DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 260 ierr = DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 261 } 262 263 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 264 Set solver options 265 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 266 ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 267 ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 268 ierr = TSSetMaxTime(ts,20.0);CHKERRQ(ierr); 269 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 270 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 271 ierr = TSMonitorLGSetVariableNames(ts,names);CHKERRQ(ierr); 272 273 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 274 Solve nonlinear system 275 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 276 ierr = TSSolve(ts,U);CHKERRQ(ierr); 277 278 ierr = TSView(ts,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 279 280 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 281 Free work space. All PETSc objects should be destroyed when they are no longer needed. 282 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 283 ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr); 284 ierr = MatDestroy(&A);CHKERRQ(ierr); 285 ierr = VecDestroy(&U);CHKERRQ(ierr); 286 ierr = TSDestroy(&ts);CHKERRQ(ierr); 287 ierr = PetscFree(adctx);CHKERRQ(ierr); 288 ierr = PetscFinalize(); 289 return ierr; 290 } 291 292 /*TEST 293 294 build: 295 requires: double !complex adolc 296 297 test: 298 suffix: 1 299 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor 300 output_file: output/adr_ex1_1.out 301 302 test: 303 suffix: 2 304 args: -ts_max_steps 1 -snes_test_jacobian 305 output_file: output/adr_ex1_2.out 306 307 TEST*/ 308