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