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