Actual source code: ex5opt_ic.c
petsc-3.12.4 2020-02-04
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*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 application 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*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*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 application 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*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*/