Actual source code: ex5opt_ic.c

petsc-3.13.4 2020-08-01
Report Typos and Errors
  1: static char help[] = "Demonstrates adjoint sensitivity analysis for Reaction-Diffusion Equations.\n";

  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.

 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:  */

 18:  #include <petscsys.h>
 19:  #include <petscdm.h>
 20:  #include <petscdmda.h>
 21:  #include <petscts.h>
 22:  #include <petsctao.h>

 24: typedef struct {
 25:   PetscScalar u,v;
 26: } Field;

 28: typedef struct {
 29:   PetscReal D1,D2,gamma,kappa;
 30:   TS        ts;
 31:   Vec       U;
 32: } AppCtx;

 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*);

 41: /*
 42:    Set terminal condition for the adjoint variable
 43:  */
 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;

 51:   VecDuplicate(U,&Uob);
 52:   PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
 53:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 54:   VecLoad(Uob,viewer);
 55:   PetscViewerDestroy(&viewer);
 56:   VecAYPX(Uob,-1.,U);
 57:   VecScale(Uob,2.0);
 58:   VecAXPY(lambda,1.,Uob);
 59:   VecDestroy(&Uob);
 60:   return(0);
 61: }

 63: /*
 64:    Set up a viewer that dumps data into a binary file
 65:  */
 66: PetscErrorCode OutputBIN(DM da, const char *filename, PetscViewer *viewer)
 67: {

 70:   PetscViewerCreate(PetscObjectComm((PetscObject)da),viewer);
 71:   PetscViewerSetType(*viewer,PETSCVIEWERBINARY);
 72:   PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);
 73:   PetscViewerFileSetName(*viewer,filename);
 74:   return(0);
 75: }

 77: /*
 78:    Generate a reference solution and save it to a binary file
 79:  */
 80: PetscErrorCode GenerateOBs(TS ts,Vec U,AppCtx *appctx)
 81: {
 82:   char           filename[PETSC_MAX_PATH_LEN] = "";
 83:   PetscViewer    viewer;
 84:   DM             da;

 87:   TSGetDM(ts,&da);
 88:   TSSolve(ts,U);
 89:   PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
 90:   OutputBIN(da,filename,&viewer);
 91:   VecView(U,viewer);
 92:   PetscViewerDestroy(&viewer);
 93:   return(0);
 94: }

 96: PetscErrorCode InitialConditions(DM da,Vec U)
 97: {
 99:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
100:   Field          **u;
101:   PetscReal      hx,hy,x,y;

104:   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);

106:   hx = 2.5/(PetscReal)Mx;
107:   hy = 2.5/(PetscReal)My;

109:   DMDAVecGetArray(da,U,&u);
110:   /* Get local grid boundaries */
111:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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;

121:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
122:     }
123:   }

125:   DMDAVecRestoreArray(da,U,&u);
126:   return(0);
127: }

129: PetscErrorCode PerturbedInitialConditions(DM da,Vec U)
130: {
132:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
133:   Field          **u;
134:   PetscReal      hx,hy,x,y;

137:   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);

139:   hx = 2.5/(PetscReal)Mx;
140:   hy = 2.5/(PetscReal)My;

142:   DMDAVecGetArray(da,U,&u);
143:   /* Get local grid boundaries */
144:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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;

154:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
155:     }
156:   }

158:   DMDAVecRestoreArray(da,U,&u);
159:   return(0);
160: }

162: PetscErrorCode PerturbedInitialConditions2(DM da,Vec U)
163: {
165:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
166:   Field          **u;
167:   PetscReal      hx,hy,x,y;

170:   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);

172:   hx = 2.5/(PetscReal)Mx;
173:   hy = 2.5/(PetscReal)My;

175:   DMDAVecGetArray(da,U,&u);
176:   /* Get local grid boundaries */
177:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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;

187:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
188:     }
189:   }

191:   DMDAVecRestoreArray(da,U,&u);
192:   return(0);
193: }

195: PetscErrorCode PerturbedInitialConditions3(DM da,Vec U)
196: {
198:   PetscInt       i,j,xs,ys,xm,ym,Mx,My;
199:   Field          **u;
200:   PetscReal      hx,hy,x,y;

203:   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);

205:   hx = 2.5/(PetscReal)Mx;
206:   hy = 2.5/(PetscReal)My;

208:   DMDAVecGetArray(da,U,&u);
209:   /* Get local grid boundaries */
210:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

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;

220:       u[j][i].u = 1.0 - 2.0*u[j][i].v;
221:     }
222:   }

224:   DMDAVecRestoreArray(da,U,&u);
225:   return(0);
226: }

228: int main(int argc,char **argv)
229: {
231:   DM             da;
232:   AppCtx         appctx;
233:   PetscBool      forwardonly = PETSC_FALSE,implicitform = PETSC_FALSE;
234:   PetscInt       perturbic = 1;

236:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
237:   PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);
238:   PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);
239:   PetscOptionsGetInt(NULL,NULL,"-perturbic",&perturbic,NULL);

241:   appctx.D1    = 8.0e-5;
242:   appctx.D2    = 4.0e-5;
243:   appctx.gamma = .024;
244:   appctx.kappa = .06;

246:   /* Create distributed array (DMDA) to manage parallel grid and vectors */
247:   DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,64,64,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
248:   DMSetFromOptions(da);
249:   DMSetUp(da);
250:   DMDASetFieldName(da,0,"u");
251:   DMDASetFieldName(da,1,"v");

253:   /* Extract global vectors from DMDA; then duplicate for remaining
254:      vectors that are the same types */
255:   DMCreateGlobalVector(da,&appctx.U);

257:   /* Create timestepping solver context */
258:   TSCreate(PETSC_COMM_WORLD,&appctx.ts);
259:   TSSetType(appctx.ts,TSCN);
260:   TSSetDM(appctx.ts,da);
261:   TSSetProblemType(appctx.ts,TS_NONLINEAR);
262:   TSSetEquationType(appctx.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
263:   if (!implicitform) {
264:     TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
265:     TSSetRHSJacobian(appctx.ts,NULL,NULL,RHSJacobian,&appctx);
266:   } else {
267:     TSSetIFunction(appctx.ts,NULL,IFunction,&appctx);
268:     TSSetIJacobian(appctx.ts,NULL,NULL,IJacobian,&appctx);
269:   }

271:   /* Set initial conditions */
272:   InitialConditions(da,appctx.U);
273:   TSSetSolution(appctx.ts,appctx.U);

275:   /* Set solver options */
276:   TSSetMaxTime(appctx.ts,2000.0);
277:   TSSetTimeStep(appctx.ts,0.5);
278:   TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
279:   TSSetFromOptions(appctx.ts);

281:   GenerateOBs(appctx.ts,appctx.U,&appctx);

283:   if (!forwardonly) {
284:     Tao tao;
285:     Vec P;
286:     Vec lambda[1];
287:     PetscLogStage opt_stage;

289:     PetscLogStageRegister("Optimization",&opt_stage);
290:     PetscLogStagePush(opt_stage);
291:     if (perturbic == 1) {
292:       PerturbedInitialConditions(da,appctx.U);
293:     } else if (perturbic == 2) {
294:       PerturbedInitialConditions2(da,appctx.U);
295:     } else if (perturbic == 3) {
296:       PerturbedInitialConditions3(da,appctx.U);
297:     }

299:     VecDuplicate(appctx.U,&lambda[0]);
300:     TSSetCostGradients(appctx.ts,1,lambda,NULL);

302:     /* Have the TS save its trajectory needed by TSAdjointSolve() */
303:     TSSetSaveTrajectory(appctx.ts);

305:     /* Create TAO solver and set desired solution method */
306:     TaoCreate(PETSC_COMM_WORLD,&tao);
307:     TaoSetType(tao,TAOBLMVM);

309:     /* Set initial guess for TAO */
310:     VecDuplicate(appctx.U,&P);
311:     VecCopy(appctx.U,P);
312:     TaoSetInitialVector(tao,P);

314:     /* Set routine for function and gradient evaluation */
315:     TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);

317:     /* Check for any TAO command line options */
318:     TaoSetFromOptions(tao);

320:     TaoSolve(tao);
321:     TaoDestroy(&tao);
322:     VecDestroy(&lambda[0]);
323:     VecDestroy(&P);
324:     PetscLogStagePop();
325:   }

327:   /* Free work space.  All PETSc objects should be destroyed when they
328:      are no longer needed. */
329:   VecDestroy(&appctx.U);
330:   TSDestroy(&appctx.ts);
331:   DMDestroy(&da);
332:   PetscFinalize();
333:   return ierr;
334: }

336: /* ------------------------ TS callbacks ---------------------------- */

338: /*
339:    RHSFunction - Evaluates nonlinear function, F(x).

341:    Input Parameters:
342: .  ts - the TS context
343: .  X - input vector
344: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

346:    Output Parameter:
347: .  F - function vector
348:  */
349: PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
350: {
351:   AppCtx         *appctx = (AppCtx*)ptr;
352:   DM             da;
354:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
355:   PetscReal      hx,hy,sx,sy;
356:   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
357:   Field          **u,**f;
358:   Vec            localU;

361:   TSGetDM(ts,&da);
362:   DMGetLocalVector(da,&localU);
363:   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);
364:   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
365:   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);

367:   /* Scatter ghost points to local vector,using the 2-step process
368:      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
369:      By placing code between these two statements, computations can be
370:      done while messages are in transition. */
371:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
372:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

374:   /* Get pointers to vector data */
375:   DMDAVecGetArrayRead(da,localU,&u);
376:   DMDAVecGetArray(da,F,&f);

378:   /* Get local grid boundaries */
379:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

381:   /* Compute function over the locally owned part of the grid */
382:   for (j=ys; j<ys+ym; j++) {
383:     for (i=xs; i<xs+xm; i++) {
384:       uc        = u[j][i].u;
385:       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
386:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
387:       vc        = u[j][i].v;
388:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
389:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
390:       f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);
391:       f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc;
392:     }
393:   }
394:   PetscLogFlops(16.0*xm*ym);

396:   DMDAVecRestoreArrayRead(da,localU,&u);
397:   DMDAVecRestoreArray(da,F,&f);
398:   DMRestoreLocalVector(da,&localU);
399:   return(0);
400: }

402: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
403: {
404:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined Section 1.5 Writing Application Codes with PETSc context */
405:   DM             da;
407:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
408:   PetscReal      hx,hy,sx,sy;
409:   PetscScalar    uc,vc;
410:   Field          **u;
411:   Vec            localU;
412:   MatStencil     stencil[6],rowstencil;
413:   PetscScalar    entries[6];

416:   TSGetDM(ts,&da);
417:   DMGetLocalVector(da,&localU);
418:   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);

420:   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
421:   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);

423:   /* Scatter ghost points to local vector,using the 2-step process
424:      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
425:      By placing code between these two statements, computations can be
426:      done while messages are in transition. */
427:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
428:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

430:   /* Get pointers to vector data */
431:   DMDAVecGetArrayRead(da,localU,&u);

433:   /* Get local grid boundaries */
434:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

436:   stencil[0].k = 0;
437:   stencil[1].k = 0;
438:   stencil[2].k = 0;
439:   stencil[3].k = 0;
440:   stencil[4].k = 0;
441:   stencil[5].k = 0;
442:   rowstencil.k = 0;

444:   /* Compute function over the locally owned part of the grid */
445:   for (j=ys; j<ys+ym; j++) {
446:     stencil[0].j = j-1;
447:     stencil[1].j = j+1;
448:     stencil[2].j = j;
449:     stencil[3].j = j;
450:     stencil[4].j = j;
451:     stencil[5].j = j;
452:     rowstencil.k = 0; rowstencil.j = j;
453:     for (i=xs; i<xs+xm; i++) {
454:       uc = u[j][i].u;
455:       vc = u[j][i].v;

457:       /* uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
458:          uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
459:          vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
460:          vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
461:          f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/

463:       stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
464:       stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
465:       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
466:       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
467:       stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
468:       stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
469:       rowstencil.i = i; rowstencil.c = 0;

471:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
472:       stencil[0].c = 1; entries[0] = appctx->D2*sy;
473:       stencil[1].c = 1; entries[1] = appctx->D2*sy;
474:       stencil[2].c = 1; entries[2] = appctx->D2*sx;
475:       stencil[3].c = 1; entries[3] = appctx->D2*sx;
476:       stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
477:       stencil[5].c = 0; entries[5] = vc*vc;
478:       rowstencil.c = 1;

480:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
481:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
482:     }
483:   }
484:   PetscLogFlops(19.0*xm*ym);

486:   DMDAVecRestoreArrayRead(da,localU,&u);
487:   DMRestoreLocalVector(da,&localU);
488:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
489:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
490:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
491:   return(0);
492: }

494: /*
495:    IFunction - Evaluates implicit nonlinear function, xdot - F(x).

497:    Input Parameters:
498: .  ts - the TS context
499: .  U - input vector
500: .  Udot - input vector
501: .  ptr - optional user-defined context, as set by TSSetRHSFunction()

503:    Output Parameter:
504: .  F - function vector
505:  */
506: PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
507: {
508:   AppCtx         *appctx = (AppCtx*)ptr;
509:   DM             da;
511:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
512:   PetscReal      hx,hy,sx,sy;
513:   PetscScalar    uc,uxx,uyy,vc,vxx,vyy;
514:   Field          **u,**f,**udot;
515:   Vec            localU;

518:   TSGetDM(ts,&da);
519:   DMGetLocalVector(da,&localU);
520:   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);
521:   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
522:   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);

524:   /* Scatter ghost points to local vector,using the 2-step process
525:      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
526:      By placing code between these two statements, computations can be
527:      done while messages are in transition. */
528:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
529:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

531:   /* Get pointers to vector data */
532:   DMDAVecGetArrayRead(da,localU,&u);
533:   DMDAVecGetArray(da,F,&f);
534:   DMDAVecGetArrayRead(da,Udot,&udot);

536:   /* Get local grid boundaries */
537:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

539:   /* Compute function over the locally owned part of the grid */
540:   for (j=ys; j<ys+ym; j++) {
541:     for (i=xs; i<xs+xm; i++) {
542:       uc        = u[j][i].u;
543:       uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
544:       uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
545:       vc        = u[j][i].v;
546:       vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
547:       vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
548:       f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc) );
549:       f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc );
550:     }
551:   }
552:   PetscLogFlops(16.0*xm*ym);

554:   DMDAVecRestoreArrayRead(da,localU,&u);
555:   DMDAVecRestoreArray(da,F,&f);
556:   DMDAVecRestoreArrayRead(da,Udot,&udot);
557:   DMRestoreLocalVector(da,&localU);
558:   return(0);
559: }

561: PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx)
562: {
563:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined Section 1.5 Writing Application Codes with PETSc context */
564:   DM             da;
566:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
567:   PetscReal      hx,hy,sx,sy;
568:   PetscScalar    uc,vc;
569:   Field          **u;
570:   Vec            localU;
571:   MatStencil     stencil[6],rowstencil;
572:   PetscScalar    entries[6];

575:   TSGetDM(ts,&da);
576:   DMGetLocalVector(da,&localU);
577:   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);

579:   hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
580:   hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);

582:   /* Scatter ghost points to local vector,using the 2-step process
583:      DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
584:      By placing code between these two statements, computations can be
585:      done while messages are in transition.*/
586:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
587:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

589:   /* Get pointers to vector data */
590:   DMDAVecGetArrayRead(da,localU,&u);

592:   /* Get local grid boundaries */
593:   DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);

595:   stencil[0].k = 0;
596:   stencil[1].k = 0;
597:   stencil[2].k = 0;
598:   stencil[3].k = 0;
599:   stencil[4].k = 0;
600:   stencil[5].k = 0;
601:   rowstencil.k = 0;

603:   /* Compute function over the locally owned part of the grid */
604:   for (j=ys; j<ys+ym; j++) {
605:     stencil[0].j = j-1;
606:     stencil[1].j = j+1;
607:     stencil[2].j = j;
608:     stencil[3].j = j;
609:     stencil[4].j = j;
610:     stencil[5].j = j;
611:     rowstencil.k = 0; rowstencil.j = j;
612:     for (i=xs; i<xs+xm; i++) {
613:       uc = u[j][i].u;
614:       vc = u[j][i].v;

616:       /* uxx       = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
617:          uyy       = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
618:          vxx       = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
619:          vyy       = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
620:          f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); */
621:       stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy;
622:       stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy;
623:       stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx;
624:       stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx;
625:       stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a;
626:       stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc;
627:       rowstencil.i = i; rowstencil.c = 0;

629:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
630:       stencil[0].c = 1; entries[0] = -appctx->D2*sy;
631:       stencil[1].c = 1; entries[1] = -appctx->D2*sy;
632:       stencil[2].c = 1; entries[2] = -appctx->D2*sx;
633:       stencil[3].c = 1; entries[3] = -appctx->D2*sx;
634:       stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a;
635:       stencil[5].c = 0; entries[5] = -vc*vc;
636:       rowstencil.c = 1;

638:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
639:       /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
640:     }
641:   }
642:   PetscLogFlops(19.0*xm*ym);

644:   DMDAVecRestoreArrayRead(da,localU,&u);
645:   DMRestoreLocalVector(da,&localU);
646:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
647:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
648:   MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
649:   return(0);
650: }

652: /* ------------------------ TAO callbacks ---------------------------- */

654: /*
655:    FormFunctionAndGradient - Evaluates the function and corresponding gradient.

657:    Input Parameters:
658:    tao - the Tao context
659:    P   - the input vector
660:    ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

662:    Output Parameters:
663:    f   - the newly evaluated function
664:    G   - the newly evaluated gradient
665: */
666: PetscErrorCode FormFunctionAndGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx)
667: {
668:   AppCtx         *appctx = (AppCtx*)ctx;
669:   PetscReal      soberr,timestep;
670:   Vec            *lambda;
671:   Vec            SDiff;
672:   DM             da;
673:   char           filename[PETSC_MAX_PATH_LEN]="";
674:   PetscViewer    viewer;

678:   TSSetTime(appctx->ts,0.0);
679:   TSGetTimeStep(appctx->ts,&timestep);
680:   if (timestep<0) {
681:     TSSetTimeStep(appctx->ts,-timestep);
682:   }
683:   TSSetStepNumber(appctx->ts,0);
684:   TSSetFromOptions(appctx->ts);

686:   VecDuplicate(P,&SDiff);
687:   VecCopy(P,appctx->U);
688:   TSGetDM(appctx->ts,&da);
689:   *f = 0;

691:   TSSolve(appctx->ts,appctx->U);
692:   PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
693:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
694:   VecLoad(SDiff,viewer);
695:   PetscViewerDestroy(&viewer);
696:   VecAYPX(SDiff,-1.,appctx->U);
697:   VecDot(SDiff,SDiff,&soberr);
698:   *f += soberr;

700:   TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);
701:   VecSet(lambda[0],0.0);
702:   InitializeLambda(da,lambda[0],appctx->U,appctx);
703:   TSAdjointSolve(appctx->ts);

705:   VecCopy(lambda[0],G);

707:   VecDestroy(&SDiff);
708:   return(0);
709: }

711: /*TEST

713:    build:
714:       requires: !complex !single

716:    test:
717:       args: -ts_max_steps 5 -ts_type rk -ts_rk_type 3 -ts_trajectory_type memory -tao_monitor -tao_view -tao_gatol 1e-6
718:       output_file: output/ex5opt_ic_1.out

720: TEST*/