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,×tep);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