Actual source code: ex5opt_ic.c

petsc-3.10.3 2018-12-18
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:   if (!implicitform) {
263:     TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
264:     TSSetRHSJacobian(appctx.ts,NULL,NULL,RHSJacobian,&appctx);
265:   } else {
266:     TSSetIFunction(appctx.ts,NULL,IFunction,&appctx);
267:     TSSetIJacobian(appctx.ts,NULL,NULL,IJacobian,&appctx);
268:   }

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

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

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

282:   if (!forwardonly) {
283:     Tao tao;
284:     Vec P;
285:     Vec lambda[1];

287:     if (perturbic == 1) {
288:       PerturbedInitialConditions(da,appctx.U);
289:     } else if (perturbic == 2) {
290:       PerturbedInitialConditions2(da,appctx.U);
291:     } else if (perturbic == 3) {
292:       PerturbedInitialConditions3(da,appctx.U);
293:     }

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

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

301:     /* Create TAO solver and set desired solution method */
302:     TaoCreate(PETSC_COMM_WORLD,&tao);
303:     TaoSetType(tao,TAOBLMVM);

305:     /* Set initial guess for TAO */
306:     VecDuplicate(appctx.U,&P);
307:     VecCopy(appctx.U,P);
308:     TaoSetInitialVector(tao,P);

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

313:     /* Check for any TAO command line options */
314:     TaoSetFromOptions(tao);

316:     /*
317:     TaoGetKSP(tao,&ksp);
318:     if (ksp) {
319:       KSPGetPC(ksp,&pc);
320:       PCSetType(pc,PCNONE);
321:     }
322:     */

324:     TaoSolve(tao);
325:     TaoDestroy(&tao);
326:     VecDestroy(&lambda[0]);
327:     VecDestroy(&P);
328:   }

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

339: /* ------------------------ TS callbacks ---------------------------- */

341: /*
342:    RHSFunction - Evaluates nonlinear function, F(x).

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

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

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

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

377:   /* Get pointers to vector data */
378:   DMDAVecGetArrayRead(da,localU,&u);
379:   DMDAVecGetArray(da,F,&f);

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

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

399:   DMDAVecRestoreArrayRead(da,localU,&u);
400:   DMDAVecRestoreArray(da,F,&f);
401:   DMRestoreLocalVector(da,&localU);
402:   return(0);
403: }

405: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
406: {
407:   AppCtx         *appctx = (AppCtx*)ctx;     /* user-defined application context */
408:   DM             da;
410:   PetscInt       i,j,Mx,My,xs,ys,xm,ym;
411:   PetscReal      hx,hy,sx,sy;
412:   PetscScalar    uc,vc;
413:   Field          **u;
414:   Vec            localU;
415:   MatStencil     stencil[6],rowstencil;
416:   PetscScalar    entries[6];

419:   TSGetDM(ts,&da);
420:   DMGetLocalVector(da,&localU);
421:   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);

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

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

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

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

439:   stencil[0].k = 0;
440:   stencil[1].k = 0;
441:   stencil[2].k = 0;
442:   stencil[3].k = 0;
443:   stencil[4].k = 0;
444:   stencil[5].k = 0;
445:   rowstencil.k = 0;

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

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

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

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

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

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

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

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

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

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

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

534:   /* Get pointers to vector data */
535:   DMDAVecGetArrayRead(da,localU,&u);
536:   DMDAVecGetArray(da,F,&f);
537:   DMDAVecGetArrayRead(da,Udot,&udot);

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

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

557:   DMDAVecRestoreArrayRead(da,localU,&u);
558:   DMDAVecRestoreArray(da,F,&f);
559:   DMDAVecRestoreArrayRead(da,Udot,&udot);
560:   DMRestoreLocalVector(da,&localU);
561:   return(0);
562: }

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

578:   TSGetDM(ts,&da);
579:   DMGetLocalVector(da,&localU);
580:   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);

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

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

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

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

598:   stencil[0].k = 0;
599:   stencil[1].k = 0;
600:   stencil[2].k = 0;
601:   stencil[3].k = 0;
602:   stencil[4].k = 0;
603:   stencil[5].k = 0;
604:   rowstencil.k = 0;

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

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

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

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

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

655: /* ------------------------ TAO callbacks ---------------------------- */

657: /*
658:    FormFunctionAndGradient - Evaluates the function and corresponding gradient.

660:    Input Parameters:
661:    tao - the Tao context
662:    P   - the input vector
663:    ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()

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

682:   TSSetTime(appctx->ts,0.0);
683:   TSGetTimeStep(appctx->ts,&timestep);
684:   if (timestep<0) {
685:     TSSetTimeStep(appctx->ts,-timestep);
686:   }
687:   TSSetStepNumber(appctx->ts,0);
688:   TSSetFromOptions(appctx->ts);

690:   VecDuplicate(P,&SDiff);
691:   VecCopy(P,appctx->U);
692:   TSGetDM(appctx->ts,&da);
693:   *f = 0;

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

704:   TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);
705:   VecSet(lambda[0],0.0);
706:   InitializeLambda(da,lambda[0],appctx->U,appctx);
707:   TSAdjointSolve(appctx->ts);

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

711:   VecDestroy(&SDiff);
712:   TSGetTrajectory(appctx->ts,&tj);
713:   TSTrajectoryReset(tj);
714:   return(0);
715: }

717: /*TEST

719:    build:
720:       requires: !complex !single

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

726: TEST*/