1 static char help[] = "Demonstrates automatic, matrix-free Jacobian generation using ADOL-C for a time-dependent PDE in 2d, solved using implicit timestepping.\n";
2
3 /*
4 Concepts: TS^time-dependent nonlinear problems
5 Concepts: Automatic differentiation using ADOL-C
6 Concepts: Matrix-free methods
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/ex5 for a description of the problem
16 ------------------------------------------------------------------------- */
17
18 #include <petscdmda.h>
19 #include <petscts.h>
20 #include "adolc-utils/init.cxx"
21 #include "adolc-utils/matfree.cxx"
22 #include <adolc/adolc.h>
23
24 /* (Passive) field for the two variables */
25 typedef struct {
26 PetscScalar u,v;
27 } Field;
28
29 /* Active field for the two variables */
30 typedef struct {
31 adouble u,v;
32 } AField;
33
34 /* Application context */
35 typedef struct {
36 PetscReal D1,D2,gamma,kappa;
37 AField **u_a,**f_a;
38 AdolcCtx *adctx; /* Automatic differentation support */
39 } AppCtx;
40
41 extern PetscErrorCode InitialConditions(DM da,Vec U);
42 extern PetscErrorCode InitializeLambda(DM da,Vec lambda,PetscReal x,PetscReal y);
43 extern PetscErrorCode IFunctionLocalPassive(DMDALocalInfo *info,PetscReal t,Field**u,Field**udot,Field**f,void *ptr);
44 extern PetscErrorCode IFunctionActive(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr);
45 extern PetscErrorCode IJacobianMatFree(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A_shell,Mat B,void *ctx);
46
main(int argc,char ** argv)47 int main(int argc,char **argv)
48 {
49 TS ts; /* ODE integrator */
50 Vec x,r; /* solution, residual */
51 PetscErrorCode ierr;
52 DM da;
53 AppCtx appctx; /* Application context */
54 AdolcMatCtx matctx; /* Matrix (free) context */
55 Vec lambda[1];
56 PetscBool forwardonly=PETSC_FALSE;
57 Mat A; /* (Matrix free) Jacobian matrix */
58 PetscInt gxm,gym;
59
60 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61 Initialize program
62 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63 ierr = PetscInitialize(&argc,&argv,NULL,help);if (ierr) return ierr;
64 ierr = PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);CHKERRQ(ierr);
65 PetscFunctionBeginUser;
66 appctx.D1 = 8.0e-5;
67 appctx.D2 = 4.0e-5;
68 appctx.gamma = .024;
69 appctx.kappa = .06;
70 ierr = PetscLogEventRegister("df/dx forward",MAT_CLASSID,&matctx.event1);CHKERRQ(ierr);
71 ierr = PetscLogEventRegister("df/d(xdot) forward",MAT_CLASSID,&matctx.event2);CHKERRQ(ierr);
72 ierr = PetscLogEventRegister("df/dx reverse",MAT_CLASSID,&matctx.event3);CHKERRQ(ierr);
73 ierr = PetscLogEventRegister("df/d(xdot) reverse",MAT_CLASSID,&matctx.event4);CHKERRQ(ierr);
74
75 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76 Create distributed array (DMDA) to manage parallel grid and vectors
77 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,65,65,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);
79 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
80 ierr = DMSetUp(da);CHKERRQ(ierr);
81 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
82 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);
83
84 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85 Extract global vectors from DMDA; then duplicate for remaining
86 vectors that are the same types
87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
89 ierr = VecDuplicate(x,&r);CHKERRQ(ierr);
90
91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92 Create matrix free context and specify usage of PETSc-ADOL-C drivers
93 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94 ierr = DMSetMatType(da,MATSHELL);CHKERRQ(ierr);
95 ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);
96 ierr = MatShellSetContext(A,&matctx);CHKERRQ(ierr);
97 ierr = MatShellSetOperation(A,MATOP_MULT,(void (*)(void))PetscAdolcIJacobianVectorProductIDMass);CHKERRQ(ierr);
98 ierr = MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void (*)(void))PetscAdolcIJacobianTransposeVectorProductIDMass);CHKERRQ(ierr);
99 ierr = VecDuplicate(x,&matctx.X);CHKERRQ(ierr);
100 ierr = VecDuplicate(x,&matctx.Xdot);CHKERRQ(ierr);
101 ierr = DMGetLocalVector(da,&matctx.localX0);CHKERRQ(ierr);
102
103 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104 Create timestepping solver context
105 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
107 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr);
108 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
109 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
110 ierr = DMDATSSetIFunctionLocal(da,INSERT_VALUES,(DMDATSIFunctionLocal)IFunctionLocalPassive,&appctx);CHKERRQ(ierr);
111
112 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
113 Some data required for matrix-free context
114 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
115 ierr = DMDAGetGhostCorners(da,NULL,NULL,NULL,&gxm,&gym,NULL);CHKERRQ(ierr);
116 matctx.m = 2*gxm*gym;matctx.n = 2*gxm*gym; /* Number of dependent and independent variables */
117 matctx.flg = PETSC_FALSE; /* Flag for reverse mode */
118 matctx.tag1 = 1; /* Tape identifier */
119
120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121 Trace function just once
122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123 ierr = PetscNew(&appctx.adctx);CHKERRQ(ierr);
124 ierr = IFunctionActive(ts,1.,x,matctx.Xdot,r,&appctx);CHKERRQ(ierr);
125 ierr = PetscFree(appctx.adctx);CHKERRQ(ierr);
126
127 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128 Set Jacobian. In this case, IJacobian simply acts to pass context
129 information to the matrix-free Jacobian vector product.
130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131 ierr = TSSetIJacobian(ts,A,A,IJacobianMatFree,&appctx);CHKERRQ(ierr);
132
133 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134 Set initial conditions
135 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136 ierr = InitialConditions(da,x);CHKERRQ(ierr);
137 ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
138
139 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140 Have the TS save its trajectory so that TSAdjointSolve() may be used
141 and set solver options
142 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143 if (!forwardonly) {
144 ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr);
145 ierr = TSSetMaxTime(ts,200.0);CHKERRQ(ierr);
146 ierr = TSSetTimeStep(ts,0.5);CHKERRQ(ierr);
147 } else {
148 ierr = TSSetMaxTime(ts,2000.0);CHKERRQ(ierr);
149 ierr = TSSetTimeStep(ts,10);CHKERRQ(ierr);
150 }
151 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
152 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
153
154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155 Solve ODE system
156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157 ierr = TSSolve(ts,x);CHKERRQ(ierr);
158 if (!forwardonly) {
159 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160 Start the Adjoint model
161 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162 ierr = VecDuplicate(x,&lambda[0]);CHKERRQ(ierr);
163 /* Reset initial conditions for the adjoint integration */
164 ierr = InitializeLambda(da,lambda[0],0.5,0.5);CHKERRQ(ierr);
165 ierr = TSSetCostGradients(ts,1,lambda,NULL);CHKERRQ(ierr);
166 ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
167 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr);
168 }
169
170 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
171 Free work space. All PETSc objects should be destroyed when they
172 are no longer needed.
173 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
174 ierr = DMRestoreLocalVector(da,&matctx.localX0);CHKERRQ(ierr);
175 ierr = VecDestroy(&r);CHKERRQ(ierr);
176 ierr = VecDestroy(&matctx.X);CHKERRQ(ierr);
177 ierr = VecDestroy(&matctx.Xdot);CHKERRQ(ierr);
178 ierr = MatDestroy(&A);CHKERRQ(ierr);
179 ierr = VecDestroy(&x);CHKERRQ(ierr);
180 ierr = TSDestroy(&ts);CHKERRQ(ierr);
181 ierr = DMDestroy(&da);CHKERRQ(ierr);
182
183 ierr = PetscFinalize();
184 return ierr;
185 }
186
InitialConditions(DM da,Vec U)187 PetscErrorCode InitialConditions(DM da,Vec U)
188 {
189 PetscErrorCode ierr;
190 PetscInt i,j,xs,ys,xm,ym,Mx,My;
191 Field **u;
192 PetscReal hx,hy,x,y;
193
194 PetscFunctionBegin;
195 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
196
197 hx = 2.5/(PetscReal)Mx;
198 hy = 2.5/(PetscReal)My;
199
200 /*
201 Get pointers to vector data
202 */
203 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
204
205 /*
206 Get local grid boundaries
207 */
208 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
209
210 /*
211 Compute function over the locally owned part of the grid
212 */
213 for (j=ys; j<ys+ym; j++) {
214 y = j*hy;
215 for (i=xs; i<xs+xm; i++) {
216 x = i*hx;
217 if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0);
218 else u[j][i].v = 0.0;
219
220 u[j][i].u = 1.0 - 2.0*u[j][i].v;
221 }
222 }
223
224 /*
225 Restore vectors
226 */
227 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
228 PetscFunctionReturn(0);
229 }
230
InitializeLambda(DM da,Vec lambda,PetscReal x,PetscReal y)231 PetscErrorCode InitializeLambda(DM da,Vec lambda,PetscReal x,PetscReal y)
232 {
233 PetscInt i,j,Mx,My,xs,ys,xm,ym;
234 PetscErrorCode ierr;
235 Field **l;
236
237 PetscFunctionBegin;
238 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
239 /* locate the global i index for x and j index for y */
240 i = (PetscInt)(x*(Mx-1));
241 j = (PetscInt)(y*(My-1));
242 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
243
244 if (xs <= i && i < xs+xm && ys <= j && j < ys+ym) {
245 /* the i,j vertex is on this process */
246 ierr = DMDAVecGetArray(da,lambda,&l);CHKERRQ(ierr);
247 l[j][i].u = 1.0;
248 l[j][i].v = 1.0;
249 ierr = DMDAVecRestoreArray(da,lambda,&l);CHKERRQ(ierr);
250 }
251 PetscFunctionReturn(0);
252 }
253
IFunctionLocalPassive(DMDALocalInfo * info,PetscReal t,Field ** u,Field ** udot,Field ** f,void * ptr)254 PetscErrorCode IFunctionLocalPassive(DMDALocalInfo *info,PetscReal t,Field**u,Field**udot,Field**f,void *ptr)
255 {
256 AppCtx *appctx = (AppCtx*)ptr;
257 PetscInt i,j,xs,ys,xm,ym;
258 PetscReal hx,hy,sx,sy;
259 PetscScalar uc,uxx,uyy,vc,vxx,vyy;
260 PetscErrorCode ierr;
261
262 PetscFunctionBegin;
263 hx = 2.50/(PetscReal)(info->mx); sx = 1.0/(hx*hx);
264 hy = 2.50/(PetscReal)(info->my); sy = 1.0/(hy*hy);
265
266 /* Get local grid boundaries */
267 xs = info->xs; xm = info->xm; ys = info->ys; ym = info->ym;
268
269 /* Compute function over the locally owned part of the grid */
270 for (j=ys; j<ys+ym; j++) {
271 for (i=xs; i<xs+xm; i++) {
272 uc = u[j][i].u;
273 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
274 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
275 vc = u[j][i].v;
276 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
277 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
278 f[j][i].u = udot[j][i].u - appctx->D1*(uxx + uyy) + uc*vc*vc - appctx->gamma*(1.0 - uc);
279 f[j][i].v = udot[j][i].v - appctx->D2*(vxx + vyy) - uc*vc*vc + (appctx->gamma + appctx->kappa)*vc;
280 }
281 }
282 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
283 PetscFunctionReturn(0);
284 }
285
IFunctionActive(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)286 PetscErrorCode IFunctionActive(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
287 {
288 PetscErrorCode ierr;
289 AppCtx *appctx = (AppCtx*)ptr;
290 DM da;
291 DMDALocalInfo info;
292 Field **u,**f,**udot;
293 Vec localU;
294 PetscInt i,j,xs,ys,xm,ym,gxs,gys,gxm,gym;
295 PetscReal hx,hy,sx,sy;
296 adouble uc,uxx,uyy,vc,vxx,vyy;
297 AField **f_a,*f_c,**u_a,*u_c;
298 PetscScalar dummy;
299
300 PetscFunctionBegin;
301 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
302 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
303 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
304 hx = 2.50/(PetscReal)(info.mx); sx = 1.0/(hx*hx);
305 hy = 2.50/(PetscReal)(info.my); sy = 1.0/(hy*hy);
306 xs = info.xs; xm = info.xm; gxs = info.gxs; gxm = info.gxm;
307 ys = info.ys; ym = info.ym; gys = info.gys; gym = info.gym;
308
309 /*
310 Scatter ghost points to local vector,using the 2-step process
311 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
312 By placing code between these two statements, computations can be
313 done while messages are in transition.
314 */
315 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
316 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
317
318 /*
319 Get pointers to vector data
320 */
321 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
322 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
323 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr);
324
325 /*
326 Create contiguous 1-arrays of AFields
327
328 NOTE: Memory for ADOL-C active variables (such as adouble and AField)
329 cannot be allocated using PetscMalloc, as this does not call the
330 relevant class constructor. Instead, we use the C++ keyword `new`.
331 */
332 u_c = new AField[info.gxm*info.gym];
333 f_c = new AField[info.gxm*info.gym];
334
335 /* Create corresponding 2-arrays of AFields */
336 u_a = new AField*[info.gym];
337 f_a = new AField*[info.gym];
338
339 /* Align indices between array types to endow 2d array with ghost points */
340 ierr = GiveGhostPoints(da,u_c,&u_a);CHKERRQ(ierr);
341 ierr = GiveGhostPoints(da,f_c,&f_a);CHKERRQ(ierr);
342
343 trace_on(1); /* Start of active section on tape 1 */
344
345 /*
346 Mark independence
347
348 NOTE: Ghost points are marked as independent, in place of the points they represent on
349 other processors / on other boundaries.
350 */
351 for (j=gys; j<gys+gym; j++) {
352 for (i=gxs; i<gxs+gxm; i++) {
353 u_a[j][i].u <<= u[j][i].u;
354 u_a[j][i].v <<= u[j][i].v;
355 }
356 }
357
358 /* Compute function over the locally owned part of the grid */
359 for (j=ys; j<ys+ym; j++) {
360 for (i=xs; i<xs+xm; i++) {
361 uc = u_a[j][i].u;
362 uxx = (-2.0*uc + u_a[j][i-1].u + u_a[j][i+1].u)*sx;
363 uyy = (-2.0*uc + u_a[j-1][i].u + u_a[j+1][i].u)*sy;
364 vc = u_a[j][i].v;
365 vxx = (-2.0*vc + u_a[j][i-1].v + u_a[j][i+1].v)*sx;
366 vyy = (-2.0*vc + u_a[j-1][i].v + u_a[j+1][i].v)*sy;
367 f_a[j][i].u = udot[j][i].u - appctx->D1*(uxx + uyy) + uc*vc*vc - appctx->gamma*(1.0 - uc);
368 f_a[j][i].v = udot[j][i].v - appctx->D2*(vxx + vyy) - uc*vc*vc + (appctx->gamma + appctx->kappa)*vc;
369 }
370 }
371
372 /*
373 Mark dependence
374
375 NOTE: Marking dependence of dummy variables makes the index notation much simpler when forming
376 the Jacobian later.
377 */
378 for (j=gys; j<gys+gym; j++) {
379 for (i=gxs; i<gxs+gxm; i++) {
380 if ((i < xs) || (i >= xs+xm) || (j < ys) || (j >= ys+ym)) {
381 f_a[j][i].u >>= dummy;
382 f_a[j][i].v >>= dummy;
383 } else {
384 f_a[j][i].u >>= f[j][i].u;
385 f_a[j][i].v >>= f[j][i].v;
386 }
387 }
388 }
389 trace_off(); /* End of active section */
390 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
391
392 /* Restore vectors */
393 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
394 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
395 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr);
396
397 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
398
399 /* Destroy AFields appropriately */
400 f_a += info.gys;
401 u_a += info.gys;
402 delete[] f_a;
403 delete[] u_a;
404 delete[] f_c;
405 delete[] u_c;
406 PetscFunctionReturn(0);
407 }
408
409 /*
410 Simply acts to pass TS information to the AdolcMatCtx
411 */
IJacobianMatFree(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A_shell,Mat B,void * ctx)412 PetscErrorCode IJacobianMatFree(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A_shell,Mat B,void *ctx)
413 {
414 AdolcMatCtx *mctx;
415 PetscErrorCode ierr;
416 DM da;
417
418 PetscFunctionBeginUser;
419 ierr = MatShellGetContext(A_shell,(void **)&mctx);CHKERRQ(ierr);
420
421 mctx->time = t;
422 mctx->shift = a;
423 if (mctx->ts != ts) mctx->ts = ts;
424 ierr = VecCopy(X,mctx->X);CHKERRQ(ierr);
425 ierr = VecCopy(Xdot,mctx->Xdot);CHKERRQ(ierr);
426 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
427 ierr = DMGlobalToLocalBegin(da,mctx->X,INSERT_VALUES,mctx->localX0);CHKERRQ(ierr);
428 ierr = DMGlobalToLocalEnd(da,mctx->X,INSERT_VALUES,mctx->localX0);CHKERRQ(ierr);
429 PetscFunctionReturn(0);
430 }
431
432 /*TEST
433
434 build:
435 requires: double !complex adolc
436
437 test:
438 suffix: 1
439 args: -ts_max_steps 1 -da_grid_x 12 -da_grid_y 12 -snes_test_jacobian
440 output_file: output/adr_ex5adj_mf_1.out
441
442 test:
443 suffix: 2
444 nsize: 4
445 args: -ts_max_steps 10 -da_grid_x 12 -da_grid_y 12 -ts_monitor -ts_adjoint_monitor
446 output_file: output/adr_ex5adj_mf_2.out
447
448 TEST*/
449