1 static char help[] = "Demonstrates adjoint sensitivity analysis for Reaction-Diffusion Equations.\n";
2 
3 /*
4    See ex5.c for details on the equation.
5    This code demonestrates the TSAdjoint/TAO interface to solve an inverse initial value problem built on a system of
6    time-dependent partial differential equations.
7    In this problem, the initial value for the PDE is unknown, but the output (the final solution of the PDE) is known.
8    We want to determine the initial value that can produce the given output.
9    We formulate the problem as a nonlinear optimization problem that minimizes the discrepency between the simulated
10    result and given reference solution, calculate the gradient of the objective function with the discrete adjoint
11    solver, and solve the optimization problem with TAO.
12 
13    Runtime options:
14      -forwardonly  - run only the forward simulation
15      -implicitform - provide IFunction and IJacobian to TS, if not set, RHSFunction and RHSJacobian will be used
16  */
17 
18 #include <petscsys.h>
19 #include <petscdm.h>
20 #include <petscdmda.h>
21 #include <petscts.h>
22 #include <petsctao.h>
23 
24 typedef struct {
25   PetscScalar u,v;
26 } Field;
27 
28 typedef struct {
29   PetscReal D1,D2,gamma,kappa;
30   TS        ts;
31   Vec       U;
32 } AppCtx;
33 
34 /* User-defined routines */
35 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*),InitialConditions(DM,Vec);
36 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
37 extern PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
38 extern PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
39 extern PetscErrorCode FormFunctionAndGradient(Tao,Vec,PetscReal*,Vec,void*);
40 
41 /*
42    Set terminal condition for the adjoint variable
43  */
InitializeLambda(DM da,Vec lambda,Vec U,AppCtx * appctx)44 PetscErrorCode InitializeLambda(DM da,Vec lambda,Vec U,AppCtx *appctx)
45 {
46   char           filename[PETSC_MAX_PATH_LEN]="";
47   PetscViewer    viewer;
48   Vec            Uob;
49   PetscErrorCode ierr;
50 
51   ierr = VecDuplicate(U,&Uob);CHKERRQ(ierr);
52   ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr);
53   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
54   ierr = VecLoad(Uob,viewer);CHKERRQ(ierr);
55   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
56   ierr = VecAYPX(Uob,-1.,U);CHKERRQ(ierr);
57   ierr = VecScale(Uob,2.0);CHKERRQ(ierr);
58   ierr = VecAXPY(lambda,1.,Uob);CHKERRQ(ierr);
59   ierr = VecDestroy(&Uob);CHKERRQ(ierr);
60   PetscFunctionReturn(0);
61 }
62 
63 /*
64    Set up a viewer that dumps data into a binary file
65  */
OutputBIN(DM da,const char * filename,PetscViewer * viewer)66 PetscErrorCode OutputBIN(DM da, const char *filename, PetscViewer *viewer)
67 {
68   PetscErrorCode ierr;
69 
70   ierr = PetscViewerCreate(PetscObjectComm((PetscObject)da),viewer);CHKERRQ(ierr);
71   ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr);
72   ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr);
73   ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr);
74   PetscFunctionReturn(0);
75 }
76 
77 /*
78    Generate a reference solution and save it to a binary file
79  */
GenerateOBs(TS ts,Vec U,AppCtx * appctx)80 PetscErrorCode GenerateOBs(TS ts,Vec U,AppCtx *appctx)
81 {
82   char           filename[PETSC_MAX_PATH_LEN] = "";
83   PetscViewer    viewer;
84   DM             da;
85   PetscErrorCode ierr;
86 
87   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
88   ierr = TSSolve(ts,U);CHKERRQ(ierr);
89   ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr);
90   ierr = OutputBIN(da,filename,&viewer);CHKERRQ(ierr);
91   ierr = VecView(U,viewer);CHKERRQ(ierr);
92   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
93   PetscFunctionReturn(0);
94 }
95 
InitialConditions(DM da,Vec U)96 PetscErrorCode InitialConditions(DM da,Vec U)
97 {
98   PetscErrorCode ierr;
99   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
100   Field          **u;
101   PetscReal      hx,hy,x,y;
102 
103   PetscFunctionBegin;
104   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);
105 
106   hx = 2.5/(PetscReal)Mx;
107   hy = 2.5/(PetscReal)My;
108 
109   ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
110   /* Get local grid boundaries */
111   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
112 
113   /* Compute function over the locally owned part of the grid */
114   for (j=ys; j<ys+ym; j++) {
115     y = j*hy;
116     for (i=xs; i<xs+xm; i++) {
117       x = i*hx;
118       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);
119       else u[j][i].v = 0.0;
120 
121       u[j][i].u = 1.0 - 2.0*u[j][i].v;
122     }
123   }
124 
125   ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
126   PetscFunctionReturn(0);
127 }
128 
PerturbedInitialConditions(DM da,Vec U)129 PetscErrorCode PerturbedInitialConditions(DM da,Vec U)
130 {
131   PetscErrorCode ierr;
132   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
133   Field          **u;
134   PetscReal      hx,hy,x,y;
135 
136   PetscFunctionBegin;
137   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);
138 
139   hx = 2.5/(PetscReal)Mx;
140   hy = 2.5/(PetscReal)My;
141 
142   ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
143   /* Get local grid boundaries */
144   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
145 
146   /* Compute function over the locally owned part of the grid */
147   for (j=ys; j<ys+ym; j++) {
148     y = j*hy;
149     for (i=xs; i<xs+xm; i++) {
150       x = i*hx;
151       if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*0.5*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*0.5*y),2.0);
152       else u[j][i].v = 0.0;
153 
154       u[j][i].u = 1.0 - 2.0*u[j][i].v;
155     }
156   }
157 
158   ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
159   PetscFunctionReturn(0);
160 }
161 
PerturbedInitialConditions2(DM da,Vec U)162 PetscErrorCode PerturbedInitialConditions2(DM da,Vec U)
163 {
164   PetscErrorCode ierr;
165   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
166   Field          **u;
167   PetscReal      hx,hy,x,y;
168 
169   PetscFunctionBegin;
170   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);
171 
172   hx = 2.5/(PetscReal)Mx;
173   hy = 2.5/(PetscReal)My;
174 
175   ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
176   /* Get local grid boundaries */
177   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
178 
179   /* Compute function over the locally owned part of the grid */
180   for (j=ys; j<ys+ym; j++) {
181     y = j*hy;
182     for (i=xs; i<xs+xm; i++) {
183       x = i*hx;
184       if ((1.0 <= x-0.5) && (x-0.5 <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*(x-0.5)),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0);
185       else u[j][i].v = 0.0;
186 
187       u[j][i].u = 1.0 - 2.0*u[j][i].v;
188     }
189   }
190 
191   ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
192   PetscFunctionReturn(0);
193 }
194 
PerturbedInitialConditions3(DM da,Vec U)195 PetscErrorCode PerturbedInitialConditions3(DM da,Vec U)
196 {
197   PetscErrorCode ierr;
198   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
199   Field          **u;
200   PetscReal      hx,hy,x,y;
201 
202   PetscFunctionBegin;
203   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);
204 
205   hx = 2.5/(PetscReal)Mx;
206   hy = 2.5/(PetscReal)My;
207 
208   ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
209   /* Get local grid boundaries */
210   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
211 
212   /* Compute function over the locally owned part of the grid */
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 ((0.5 <= x) && (x <= 2.0) && (0.5 <= y) && (y <= 2.0)) 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   ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
225   PetscFunctionReturn(0);
226 }
227 
main(int argc,char ** argv)228 int main(int argc,char **argv)
229 {
230   PetscErrorCode ierr;
231   DM             da;
232   AppCtx         appctx;
233   PetscBool      forwardonly = PETSC_FALSE,implicitform = PETSC_FALSE;
234   PetscInt       perturbic = 1;
235 
236   ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
237   ierr = PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);CHKERRQ(ierr);
238   ierr = PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);CHKERRQ(ierr);
239   ierr = PetscOptionsGetInt(NULL,NULL,"-perturbic",&perturbic,NULL);CHKERRQ(ierr);
240 
241   appctx.D1    = 8.0e-5;
242   appctx.D2    = 4.0e-5;
243   appctx.gamma = .024;
244   appctx.kappa = .06;
245 
246   /* Create distributed array (DMDA) to manage parallel grid and vectors */
247   ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,64,64,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);CHKERRQ(ierr);
248   ierr = DMSetFromOptions(da);CHKERRQ(ierr);
249   ierr = DMSetUp(da);CHKERRQ(ierr);
250   ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
251   ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);
252 
253   /* Extract global vectors from DMDA; then duplicate for remaining
254      vectors that are the same types */
255   ierr = DMCreateGlobalVector(da,&appctx.U);CHKERRQ(ierr);
256 
257   /* Create timestepping solver context */
258   ierr = TSCreate(PETSC_COMM_WORLD,&appctx.ts);CHKERRQ(ierr);
259   ierr = TSSetType(appctx.ts,TSCN);CHKERRQ(ierr);
260   ierr = TSSetDM(appctx.ts,da);CHKERRQ(ierr);
261   ierr = TSSetProblemType(appctx.ts,TS_NONLINEAR);CHKERRQ(ierr);
262   ierr = TSSetEquationType(appctx.ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
263   if (!implicitform) {
264     ierr = TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);CHKERRQ(ierr);
265     ierr = TSSetRHSJacobian(appctx.ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
266   } else {
267     ierr = TSSetIFunction(appctx.ts,NULL,IFunction,&appctx);CHKERRQ(ierr);
268     ierr = TSSetIJacobian(appctx.ts,NULL,NULL,IJacobian,&appctx);CHKERRQ(ierr);
269   }
270 
271   /* Set initial conditions */
272   ierr = InitialConditions(da,appctx.U);CHKERRQ(ierr);
273   ierr = TSSetSolution(appctx.ts,appctx.U);CHKERRQ(ierr);
274 
275   /* Set solver options */
276   ierr = TSSetMaxTime(appctx.ts,2000.0);CHKERRQ(ierr);
277   ierr = TSSetTimeStep(appctx.ts,0.5);CHKERRQ(ierr);
278   ierr = TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
279   ierr = TSSetFromOptions(appctx.ts);CHKERRQ(ierr);
280 
281   ierr = GenerateOBs(appctx.ts,appctx.U,&appctx);CHKERRQ(ierr);
282 
283   if (!forwardonly) {
284     Tao           tao;
285     Vec           P;
286     Vec           lambda[1];
287 #if defined(PETSC_USE_LOG)
288     PetscLogStage opt_stage;
289 #endif
290 
291     ierr = PetscLogStageRegister("Optimization",&opt_stage);CHKERRQ(ierr);
292     ierr = PetscLogStagePush(opt_stage);CHKERRQ(ierr);
293     if (perturbic == 1) {
294       ierr = PerturbedInitialConditions(da,appctx.U);CHKERRQ(ierr);
295     } else if (perturbic == 2) {
296       ierr = PerturbedInitialConditions2(da,appctx.U);CHKERRQ(ierr);
297     } else if (perturbic == 3) {
298       ierr = PerturbedInitialConditions3(da,appctx.U);CHKERRQ(ierr);
299     }
300 
301     ierr = VecDuplicate(appctx.U,&lambda[0]);CHKERRQ(ierr);
302     ierr = TSSetCostGradients(appctx.ts,1,lambda,NULL);CHKERRQ(ierr);
303 
304     /* Have the TS save its trajectory needed by TSAdjointSolve() */
305     ierr = TSSetSaveTrajectory(appctx.ts);CHKERRQ(ierr);
306 
307     /* Create TAO solver and set desired solution method */
308     ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr);
309     ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr);
310 
311     /* Set initial guess for TAO */
312     ierr = VecDuplicate(appctx.U,&P);CHKERRQ(ierr);
313     ierr = VecCopy(appctx.U,P);CHKERRQ(ierr);
314     ierr = TaoSetInitialVector(tao,P);CHKERRQ(ierr);
315 
316     /* Set routine for function and gradient evaluation */
317     ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);CHKERRQ(ierr);
318 
319     /* Check for any TAO command line options */
320     ierr = TaoSetFromOptions(tao);CHKERRQ(ierr);
321 
322     ierr = TaoSolve(tao);CHKERRQ(ierr);
323     ierr = TaoDestroy(&tao);CHKERRQ(ierr);
324     ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr);
325     ierr = VecDestroy(&P);CHKERRQ(ierr);
326     ierr = PetscLogStagePop();CHKERRQ(ierr);
327   }
328 
329   /* Free work space.  All PETSc objects should be destroyed when they
330      are no longer needed. */
331   ierr = VecDestroy(&appctx.U);CHKERRQ(ierr);
332   ierr = TSDestroy(&appctx.ts);CHKERRQ(ierr);
333   ierr = DMDestroy(&da);CHKERRQ(ierr);
334   ierr = PetscFinalize();
335   return ierr;
336 }
337 
338 /* ------------------------ TS callbacks ---------------------------- */
339 
340 /*
341    RHSFunction - Evaluates nonlinear function, F(x).
342 
343    Input Parameters:
344 .  ts - the TS context
345 .  X - input vector
346 .  ptr - optional user-defined context, as set by TSSetRHSFunction()
347 
348    Output Parameter:
349 .  F - function vector
350  */
RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void * ptr)351 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
352 {
353   AppCtx         *appctx = (AppCtx*)ptr;
354   DM             da;
355   PetscErrorCode ierr;
356   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
357   PetscReal      hx,hy,sx,sy;
358   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
359   Field          **u,**f;
360   Vec            localU;
361 
362   PetscFunctionBegin;
363   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
364   ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
365   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);
366   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
367   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
368 
369   /* Scatter ghost points to local vector,using the 2-step process
370      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
371      By placing code between these two statements, computations can be
372      done while messages are in transition. */
373   ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
374   ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
375 
376   /* Get pointers to vector data */
377   ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
378   ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
379 
380   /* Get local grid boundaries */
381   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
382 
383   /* Compute function over the locally owned part of the grid */
384   for (j=ys; j<ys+ym; j++) {
385     for (i=xs; i<xs+xm; i++) {
386       uc        = u[j][i].u;
387       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
388       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
389       vc        = u[j][i].v;
390       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
391       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
392       f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);
393       f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc;
394     }
395   }
396   ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
397 
398   ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
399   ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
400   ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
401   PetscFunctionReturn(0);
402 }
403 
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void * ctx)404 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
405 {
406   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
407   DM             da;
408   PetscErrorCode ierr;
409   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
410   PetscReal      hx,hy,sx,sy;
411   PetscScalar    uc,vc;
412   Field          **u;
413   Vec            localU;
414   MatStencil     stencil[6],rowstencil;
415   PetscScalar    entries[6];
416 
417   PetscFunctionBegin;
418   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
419   ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
420   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);
421 
422   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
423   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
424 
425   /* Scatter ghost points to local vector,using the 2-step process
426      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
427      By placing code between these two statements, computations can be
428      done while messages are in transition. */
429   ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
430   ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
431 
432   /* Get pointers to vector data */
433   ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
434 
435   /* Get local grid boundaries */
436   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
437 
438   stencil[0].k = 0;
439   stencil[1].k = 0;
440   stencil[2].k = 0;
441   stencil[3].k = 0;
442   stencil[4].k = 0;
443   stencil[5].k = 0;
444   rowstencil.k = 0;
445 
446   /* Compute function over the locally owned part of the grid */
447   for (j=ys; j<ys+ym; j++) {
448     stencil[0].j = j-1;
449     stencil[1].j = j+1;
450     stencil[2].j = j;
451     stencil[3].j = j;
452     stencil[4].j = j;
453     stencil[5].j = j;
454     rowstencil.k = 0; rowstencil.j = j;
455     for (i=xs; i<xs+xm; i++) {
456       uc = u[j][i].u;
457       vc = u[j][i].v;
458 
459       /* uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
460          uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
461          vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
462          vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
463          f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/
464 
465       stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
466       stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
467       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
468       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
469       stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
470       stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
471       rowstencil.i = i; rowstencil.c = 0;
472 
473       ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
474       stencil[0].c = 1; entries[0] = appctx->D2*sy;
475       stencil[1].c = 1; entries[1] = appctx->D2*sy;
476       stencil[2].c = 1; entries[2] = appctx->D2*sx;
477       stencil[3].c = 1; entries[3] = appctx->D2*sx;
478       stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
479       stencil[5].c = 0; entries[5] = vc*vc;
480       rowstencil.c = 1;
481 
482       ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
483       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
484     }
485   }
486   ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr);
487 
488   ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
489   ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
490   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
491   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
492   ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
493   PetscFunctionReturn(0);
494 }
495 
496 /*
497    IFunction - Evaluates implicit nonlinear function, xdot - F(x).
498 
499    Input Parameters:
500 .  ts - the TS context
501 .  U - input vector
502 .  Udot - input vector
503 .  ptr - optional user-defined context, as set by TSSetRHSFunction()
504 
505    Output Parameter:
506 .  F - function vector
507  */
IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)508 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
509 {
510   AppCtx         *appctx = (AppCtx*)ptr;
511   DM             da;
512   PetscErrorCode ierr;
513   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
514   PetscReal      hx,hy,sx,sy;
515   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
516   Field          **u,**f,**udot;
517   Vec            localU;
518 
519   PetscFunctionBegin;
520   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
521   ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
522   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);
523   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
524   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
525 
526   /* Scatter ghost points to local vector,using the 2-step process
527      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
528      By placing code between these two statements, computations can be
529      done while messages are in transition. */
530   ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
531   ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
532 
533   /* Get pointers to vector data */
534   ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
535   ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
536   ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr);
537 
538   /* Get local grid boundaries */
539   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
540 
541   /* Compute function over the locally owned part of the grid */
542   for (j=ys; j<ys+ym; j++) {
543     for (i=xs; i<xs+xm; i++) {
544       uc        = u[j][i].u;
545       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
546       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
547       vc        = u[j][i].v;
548       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
549       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
550       f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc));
551       f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc);
552     }
553   }
554   ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
555 
556   ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
557   ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
558   ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr);
559   ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
560   PetscFunctionReturn(0);
561 }
562 
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void * ctx)563 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx)
564 {
565   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
566   DM             da;
567   PetscErrorCode ierr;
568   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
569   PetscReal      hx,hy,sx,sy;
570   PetscScalar    uc,vc;
571   Field          **u;
572   Vec            localU;
573   MatStencil     stencil[6],rowstencil;
574   PetscScalar    entries[6];
575 
576   PetscFunctionBegin;
577   ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
578   ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
579   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);
580 
581   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
582   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
583 
584   /* Scatter ghost points to local vector,using the 2-step process
585      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
586      By placing code between these two statements, computations can be
587      done while messages are in transition.*/
588   ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
589   ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
590 
591   /* Get pointers to vector data */
592   ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
593 
594   /* Get local grid boundaries */
595   ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
596 
597   stencil[0].k = 0;
598   stencil[1].k = 0;
599   stencil[2].k = 0;
600   stencil[3].k = 0;
601   stencil[4].k = 0;
602   stencil[5].k = 0;
603   rowstencil.k = 0;
604 
605   /* Compute function over the locally owned part of the grid */
606   for (j=ys; j<ys+ym; j++) {
607     stencil[0].j = j-1;
608     stencil[1].j = j+1;
609     stencil[2].j = j;
610     stencil[3].j = j;
611     stencil[4].j = j;
612     stencil[5].j = j;
613     rowstencil.k = 0; rowstencil.j = j;
614     for (i=xs; i<xs+xm; i++) {
615       uc = u[j][i].u;
616       vc = u[j][i].v;
617 
618       /* uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
619          uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
620          vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
621          vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
622          f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); */
623       stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy;
624       stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy;
625       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx;
626       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx;
627       stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a;
628       stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc;
629       rowstencil.i = i; rowstencil.c = 0;
630 
631       ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
632       stencil[0].c = 1; entries[0] = -appctx->D2*sy;
633       stencil[1].c = 1; entries[1] = -appctx->D2*sy;
634       stencil[2].c = 1; entries[2] = -appctx->D2*sx;
635       stencil[3].c = 1; entries[3] = -appctx->D2*sx;
636       stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a;
637       stencil[5].c = 0; entries[5] = -vc*vc;
638       rowstencil.c = 1;
639 
640       ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
641       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
642     }
643   }
644   ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr);
645 
646   ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
647   ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
648   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
649   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
650   ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
651   PetscFunctionReturn(0);
652 }
653 
654 /* ------------------------ TAO callbacks ---------------------------- */
655 
656 /*
657    FormFunctionAndGradient - Evaluates the function and corresponding gradient.
658 
659    Input Parameters:
660    tao - the Tao context
661    P   - the input vector
662    ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
663 
664    Output Parameters:
665    f   - the newly evaluated function
666    G   - the newly evaluated gradient
667 */
FormFunctionAndGradient(Tao tao,Vec P,PetscReal * f,Vec G,void * ctx)668 PetscErrorCode FormFunctionAndGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
669 {
670   AppCtx         *appctx = (AppCtx*)ctx;
671   PetscReal      soberr,timestep;
672   Vec            *lambda;
673   Vec            SDiff;
674   DM             da;
675   char           filename[PETSC_MAX_PATH_LEN]="";
676   PetscViewer    viewer;
677   PetscErrorCode ierr;
678 
679   PetscFunctionBeginUser;
680   ierr = TSSetTime(appctx->ts,0.0);CHKERRQ(ierr);
681   ierr = TSGetTimeStep(appctx->ts,&timestep);CHKERRQ(ierr);
682   if (timestep<0) {
683     ierr = TSSetTimeStep(appctx->ts,-timestep);CHKERRQ(ierr);
684   }
685   ierr = TSSetStepNumber(appctx->ts,0);CHKERRQ(ierr);
686   ierr = TSSetFromOptions(appctx->ts);CHKERRQ(ierr);
687 
688   ierr = VecDuplicate(P,&SDiff);CHKERRQ(ierr);
689   ierr = VecCopy(P,appctx->U);CHKERRQ(ierr);
690   ierr = TSGetDM(appctx->ts,&da);CHKERRQ(ierr);
691   *f = 0;
692 
693   ierr = TSSolve(appctx->ts,appctx->U);CHKERRQ(ierr);
694   ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr);
695   ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr);
696   ierr = VecLoad(SDiff,viewer);CHKERRQ(ierr);
697   ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
698   ierr = VecAYPX(SDiff,-1.,appctx->U);CHKERRQ(ierr);
699   ierr = VecDot(SDiff,SDiff,&soberr);CHKERRQ(ierr);
700   *f += soberr;
701 
702   ierr = TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);CHKERRQ(ierr);
703   ierr = VecSet(lambda[0],0.0);CHKERRQ(ierr);
704   ierr = InitializeLambda(da,lambda[0],appctx->U,appctx);CHKERRQ(ierr);
705   ierr = TSAdjointSolve(appctx->ts);CHKERRQ(ierr);
706 
707   ierr = VecCopy(lambda[0],G);CHKERRQ(ierr);
708 
709   ierr = VecDestroy(&SDiff);CHKERRQ(ierr);
710   PetscFunctionReturn(0);
711 }
712 
713 /*TEST
714 
715    build:
716       requires: !complex !single
717 
718    test:
719       args: -ts_max_steps 5 -ts_type rk -ts_rk_type 3 -ts_trajectory_type memory -tao_monitor -tao_view -tao_gatol 1e-6
720       output_file: output/ex5opt_ic_1.out
721 
722 TEST*/
723