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