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