Actual source code: ex5opt_ic.c

petsc-3.14.3 2021-01-09
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: #if defined(PETSC_USE_LOG)
288:     PetscLogStage opt_stage;
289: #endif

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

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

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

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

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

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

319:     /* Check for any TAO command line options */
320:     TaoSetFromOptions(tao);

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

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

338: /* ------------------------ TS callbacks ---------------------------- */

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

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

348:    Output Parameter:
349: .  F - function vector
350:  */
351: PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
352: {
353:   AppCtx         *appctx = (AppCtx*)ptr;
354:   DM             da;
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;

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

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:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
374:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

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

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

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:   PetscLogFlops(16.0*xm*ym);

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

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;
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];

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

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

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:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
430:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

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

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

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;

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;

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

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;

473:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;

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

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

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

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()

505:    Output Parameter:
506: .  F - function vector
507:  */
508: PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
509: {
510:   AppCtx         *appctx = (AppCtx*)ptr;
511:   DM             da;
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;

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

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:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
531:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

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

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

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:   PetscLogFlops(16.0*xm*ym);

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

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;
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];

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

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

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:   DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
589:   DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);

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

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

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;

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;

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;

631:       MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;

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

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

654: /* ------------------------ TAO callbacks ---------------------------- */

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

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

664:    Output Parameters:
665:    f   - the newly evaluated function
666:    G   - the newly evaluated gradient
667: */
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;

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

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

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

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

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

709:   VecDestroy(&SDiff);
710:   return(0);
711: }

713: /*TEST

715:    build:
716:       requires: !complex !single

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

722: TEST*/