Actual source code: ex5opt_ic.c

  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 discrepancy 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 "reaction_diffusion.h"
 19: #include <petscdm.h>
 20: #include <petscdmda.h>
 21: #include <petsctao.h>

 23: /* User-defined routines */
 24: extern PetscErrorCode FormFunctionAndGradient(Tao, Vec, PetscReal *, Vec, void *);

 26: /*
 27:    Set terminal condition for the adjoint variable
 28:  */
 29: PetscErrorCode InitializeLambda(DM da, Vec lambda, Vec U, AppCtx *appctx)
 30: {
 31:   char        filename[PETSC_MAX_PATH_LEN] = "";
 32:   PetscViewer viewer;
 33:   Vec         Uob;

 35:   VecDuplicate(U, &Uob);
 36:   PetscSNPrintf(filename, sizeof filename, "ex5opt.ob");
 37:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);
 38:   VecLoad(Uob, viewer);
 39:   PetscViewerDestroy(&viewer);
 40:   VecAYPX(Uob, -1., U);
 41:   VecScale(Uob, 2.0);
 42:   VecAXPY(lambda, 1., Uob);
 43:   VecDestroy(&Uob);
 44:   return 0;
 45: }

 47: /*
 48:    Set up a viewer that dumps data into a binary file
 49:  */
 50: PetscErrorCode OutputBIN(DM da, const char *filename, PetscViewer *viewer)
 51: {
 52:   PetscViewerCreate(PetscObjectComm((PetscObject)da), viewer);
 53:   PetscViewerSetType(*viewer, PETSCVIEWERBINARY);
 54:   PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE);
 55:   PetscViewerFileSetName(*viewer, filename);
 56:   return 0;
 57: }

 59: /*
 60:    Generate a reference solution and save it to a binary file
 61:  */
 62: PetscErrorCode GenerateOBs(TS ts, Vec U, AppCtx *appctx)
 63: {
 64:   char        filename[PETSC_MAX_PATH_LEN] = "";
 65:   PetscViewer viewer;
 66:   DM          da;

 68:   TSGetDM(ts, &da);
 69:   TSSolve(ts, U);
 70:   PetscSNPrintf(filename, sizeof filename, "ex5opt.ob");
 71:   OutputBIN(da, filename, &viewer);
 72:   VecView(U, viewer);
 73:   PetscViewerDestroy(&viewer);
 74:   return 0;
 75: }

 77: PetscErrorCode InitialConditions(DM da, Vec U)
 78: {
 79:   PetscInt  i, j, xs, ys, xm, ym, Mx, My;
 80:   Field   **u;
 81:   PetscReal hx, hy, x, y;

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

 85:   hx = 2.5 / (PetscReal)Mx;
 86:   hy = 2.5 / (PetscReal)My;

 88:   DMDAVecGetArray(da, U, &u);
 89:   /* Get local grid boundaries */
 90:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

 92:   /* Compute function over the locally owned part of the grid */
 93:   for (j = ys; j < ys + ym; j++) {
 94:     y = j * hy;
 95:     for (i = xs; i < xs + xm; i++) {
 96:       x = i * hx;
 97:       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);
 98:       else u[j][i].v = 0.0;

100:       u[j][i].u = 1.0 - 2.0 * u[j][i].v;
101:     }
102:   }

104:   DMDAVecRestoreArray(da, U, &u);
105:   return 0;
106: }

108: PetscErrorCode PerturbedInitialConditions(DM da, Vec U)
109: {
110:   PetscInt  i, j, xs, ys, xm, ym, Mx, My;
111:   Field   **u;
112:   PetscReal hx, hy, x, y;

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

116:   hx = 2.5 / (PetscReal)Mx;
117:   hy = 2.5 / (PetscReal)My;

119:   DMDAVecGetArray(da, U, &u);
120:   /* Get local grid boundaries */
121:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

123:   /* Compute function over the locally owned part of the grid */
124:   for (j = ys; j < ys + ym; j++) {
125:     y = j * hy;
126:     for (i = xs; i < xs + xm; i++) {
127:       x = i * hx;
128:       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);
129:       else u[j][i].v = 0.0;

131:       u[j][i].u = 1.0 - 2.0 * u[j][i].v;
132:     }
133:   }

135:   DMDAVecRestoreArray(da, U, &u);
136:   return 0;
137: }

139: PetscErrorCode PerturbedInitialConditions2(DM da, Vec U)
140: {
141:   PetscInt  i, j, xs, ys, xm, ym, Mx, My;
142:   Field   **u;
143:   PetscReal hx, hy, x, y;

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

147:   hx = 2.5 / (PetscReal)Mx;
148:   hy = 2.5 / (PetscReal)My;

150:   DMDAVecGetArray(da, U, &u);
151:   /* Get local grid boundaries */
152:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

154:   /* Compute function over the locally owned part of the grid */
155:   for (j = ys; j < ys + ym; j++) {
156:     y = j * hy;
157:     for (i = xs; i < xs + xm; i++) {
158:       x = i * hx;
159:       if (PetscApproximateGTE(x, 1.0) && PetscApproximateLTE(x, 1.5) && PetscApproximateGTE(y, 1.0) && PetscApproximateLTE(y, 1.5))
160:         u[j][i].v = PetscPowReal(PetscSinReal(4.0 * PETSC_PI * (x - 0.5)), 2.0) * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * y), 2.0) / 4.0;
161:       else u[j][i].v = 0.0;

163:       u[j][i].u = 1.0 - 2.0 * u[j][i].v;
164:     }
165:   }

167:   DMDAVecRestoreArray(da, U, &u);
168:   return 0;
169: }

171: PetscErrorCode PerturbedInitialConditions3(DM da, Vec U)
172: {
173:   PetscInt  i, j, xs, ys, xm, ym, Mx, My;
174:   Field   **u;
175:   PetscReal hx, hy, x, y;

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

179:   hx = 2.5 / (PetscReal)Mx;
180:   hy = 2.5 / (PetscReal)My;

182:   DMDAVecGetArray(da, U, &u);
183:   /* Get local grid boundaries */
184:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

186:   /* Compute function over the locally owned part of the grid */
187:   for (j = ys; j < ys + ym; j++) {
188:     y = j * hy;
189:     for (i = xs; i < xs + xm; i++) {
190:       x = i * hx;
191:       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);
192:       else u[j][i].v = 0.0;

194:       u[j][i].u = 1.0 - 2.0 * u[j][i].v;
195:     }
196:   }

198:   DMDAVecRestoreArray(da, U, &u);
199:   return 0;
200: }

202: int main(int argc, char **argv)
203: {
204:   DM        da;
205:   AppCtx    appctx;
206:   PetscBool forwardonly = PETSC_FALSE, implicitform = PETSC_FALSE;
207:   PetscInt  perturbic = 1;

210:   PetscInitialize(&argc, &argv, (char *)0, help);
211:   PetscOptionsGetBool(NULL, NULL, "-forwardonly", &forwardonly, NULL);
212:   PetscOptionsGetBool(NULL, NULL, "-implicitform", &implicitform, NULL);
213:   PetscOptionsGetInt(NULL, NULL, "-perturbic", &perturbic, NULL);

215:   appctx.D1    = 8.0e-5;
216:   appctx.D2    = 4.0e-5;
217:   appctx.gamma = .024;
218:   appctx.kappa = .06;
219:   appctx.aijpc = PETSC_FALSE;

221:   /* Create distributed array (DMDA) to manage parallel grid and vectors */
222:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 64, 64, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da);
223:   DMSetFromOptions(da);
224:   DMSetUp(da);
225:   DMDASetFieldName(da, 0, "u");
226:   DMDASetFieldName(da, 1, "v");

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

232:   /* Create timestepping solver context */
233:   TSCreate(PETSC_COMM_WORLD, &appctx.ts);
234:   TSSetType(appctx.ts, TSCN);
235:   TSSetDM(appctx.ts, da);
236:   TSSetProblemType(appctx.ts, TS_NONLINEAR);
237:   TSSetEquationType(appctx.ts, TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
238:   if (!implicitform) {
239:     TSSetRHSFunction(appctx.ts, NULL, RHSFunction, &appctx);
240:     TSSetRHSJacobian(appctx.ts, NULL, NULL, RHSJacobian, &appctx);
241:   } else {
242:     TSSetIFunction(appctx.ts, NULL, IFunction, &appctx);
243:     TSSetIJacobian(appctx.ts, NULL, NULL, IJacobian, &appctx);
244:   }

246:   /* Set initial conditions */
247:   InitialConditions(da, appctx.U);
248:   TSSetSolution(appctx.ts, appctx.U);

250:   /* Set solver options */
251:   TSSetMaxTime(appctx.ts, 2000.0);
252:   TSSetTimeStep(appctx.ts, 0.5);
253:   TSSetExactFinalTime(appctx.ts, TS_EXACTFINALTIME_MATCHSTEP);
254:   TSSetFromOptions(appctx.ts);

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

258:   if (!forwardonly) {
259:     Tao tao;
260:     Vec P;
261:     Vec lambda[1];
262: #if defined(PETSC_USE_LOG)
263:     PetscLogStage opt_stage;
264: #endif

266:     PetscLogStageRegister("Optimization", &opt_stage);
267:     PetscLogStagePush(opt_stage);
268:     if (perturbic == 1) {
269:       PerturbedInitialConditions(da, appctx.U);
270:     } else if (perturbic == 2) {
271:       PerturbedInitialConditions2(da, appctx.U);
272:     } else if (perturbic == 3) {
273:       PerturbedInitialConditions3(da, appctx.U);
274:     }

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

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

282:     /* Create TAO solver and set desired solution method */
283:     TaoCreate(PETSC_COMM_WORLD, &tao);
284:     TaoSetType(tao, TAOBLMVM);

286:     /* Set initial guess for TAO */
287:     VecDuplicate(appctx.U, &P);
288:     VecCopy(appctx.U, P);
289:     TaoSetSolution(tao, P);

291:     /* Set routine for function and gradient evaluation */
292:     TaoSetObjectiveAndGradient(tao, NULL, FormFunctionAndGradient, &appctx);

294:     /* Check for any TAO command line options */
295:     TaoSetFromOptions(tao);

297:     TaoSolve(tao);
298:     TaoDestroy(&tao);
299:     VecDestroy(&lambda[0]);
300:     VecDestroy(&P);
301:     PetscLogStagePop();
302:   }

304:   /* Free work space.  All PETSc objects should be destroyed when they
305:      are no longer needed. */
306:   VecDestroy(&appctx.U);
307:   TSDestroy(&appctx.ts);
308:   DMDestroy(&da);
309:   PetscFinalize();
310:   return 0;
311: }

313: /* ------------------------ TAO callbacks ---------------------------- */

315: /*
316:    FormFunctionAndGradient - Evaluates the function and corresponding gradient.

318:    Input Parameters:
319:    tao - the Tao context
320:    P   - the input vector
321:    ctx - optional user-defined context, as set by TaoSetObjectiveAndGradient()

323:    Output Parameters:
324:    f   - the newly evaluated function
325:    G   - the newly evaluated gradient
326: */
327: PetscErrorCode FormFunctionAndGradient(Tao tao, Vec P, PetscReal *f, Vec G, void *ctx)
328: {
329:   AppCtx     *appctx = (AppCtx *)ctx;
330:   PetscReal   soberr, timestep;
331:   Vec        *lambda;
332:   Vec         SDiff;
333:   DM          da;
334:   char        filename[PETSC_MAX_PATH_LEN] = "";
335:   PetscViewer viewer;

338:   TSSetTime(appctx->ts, 0.0);
339:   TSGetTimeStep(appctx->ts, &timestep);
340:   if (timestep < 0) TSSetTimeStep(appctx->ts, -timestep);
341:   TSSetStepNumber(appctx->ts, 0);
342:   TSSetFromOptions(appctx->ts);

344:   VecDuplicate(P, &SDiff);
345:   VecCopy(P, appctx->U);
346:   TSGetDM(appctx->ts, &da);
347:   *f = 0;

349:   TSSolve(appctx->ts, appctx->U);
350:   PetscSNPrintf(filename, sizeof filename, "ex5opt.ob");
351:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);
352:   VecLoad(SDiff, viewer);
353:   PetscViewerDestroy(&viewer);
354:   VecAYPX(SDiff, -1., appctx->U);
355:   VecDot(SDiff, SDiff, &soberr);
356:   *f += soberr;

358:   TSGetCostGradients(appctx->ts, NULL, &lambda, NULL);
359:   VecSet(lambda[0], 0.0);
360:   InitializeLambda(da, lambda[0], appctx->U, appctx);
361:   TSAdjointSolve(appctx->ts);

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

365:   VecDestroy(&SDiff);
366:   return 0;
367: }

369: /*TEST

371:    build:
372:       depends: reaction_diffusion.c
373:       requires: !complex !single

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

379: TEST*/