Actual source code: ex5adj.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 interface to a system of time-dependent partial differential equations.
6: It computes the sensitivity of a component in the final solution, which locates in the center of the 2D domain, w.r.t. the initial conditions.
7: The user does not need to provide any additional functions. The required functions in the original simulation are reused in the adjoint run.
9: Runtime options:
10: -forwardonly - run the forward simulation without adjoint
11: -implicitform - provide IFunction and IJacobian to TS, if not set, RHSFunction and RHSJacobian will be used
12: -aijpc - set the preconditioner matrix to be aij (the Jacobian matrix can be of a different type such as ELL)
13: */
14: #include "reaction_diffusion.h"
15: #include <petscdm.h>
16: #include <petscdmda.h>
18: PetscErrorCode InitialConditions(DM, Vec);
20: PetscErrorCode InitializeLambda(DM da, Vec lambda, PetscReal x, PetscReal y)
21: {
22: PetscInt i, j, Mx, My, xs, ys, xm, ym;
23: Field **l;
25: 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);
26: /* locate the global i index for x and j index for y */
27: i = (PetscInt)(x * (Mx - 1));
28: j = (PetscInt)(y * (My - 1));
29: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
31: if (xs <= i && i < xs + xm && ys <= j && j < ys + ym) {
32: /* the i,j vertex is on this process */
33: DMDAVecGetArray(da, lambda, &l);
34: l[j][i].u = 1.0;
35: l[j][i].v = 1.0;
36: DMDAVecRestoreArray(da, lambda, &l);
37: }
38: return 0;
39: }
41: int main(int argc, char **argv)
42: {
43: TS ts; /* ODE integrator */
44: Vec x; /* solution */
45: DM da;
46: AppCtx appctx;
47: Vec lambda[1];
48: PetscBool forwardonly = PETSC_FALSE, implicitform = PETSC_TRUE;
51: PetscInitialize(&argc, &argv, (char *)0, help);
52: PetscOptionsGetBool(NULL, NULL, "-forwardonly", &forwardonly, NULL);
53: PetscOptionsGetBool(NULL, NULL, "-implicitform", &implicitform, NULL);
54: appctx.aijpc = PETSC_FALSE;
55: PetscOptionsGetBool(NULL, NULL, "-aijpc", &appctx.aijpc, NULL);
57: appctx.D1 = 8.0e-5;
58: appctx.D2 = 4.0e-5;
59: appctx.gamma = .024;
60: appctx.kappa = .06;
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Create distributed array (DMDA) to manage parallel grid and vectors
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 64, 64, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, NULL, &da);
66: DMSetFromOptions(da);
67: DMSetUp(da);
68: DMDASetFieldName(da, 0, "u");
69: DMDASetFieldName(da, 1, "v");
71: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
72: Extract global vectors from DMDA; then duplicate for remaining
73: vectors that are the same types
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75: DMCreateGlobalVector(da, &x);
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Create timestepping solver context
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: TSCreate(PETSC_COMM_WORLD, &ts);
81: TSSetDM(ts, da);
82: TSSetProblemType(ts, TS_NONLINEAR);
83: TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
84: if (!implicitform) {
85: TSSetType(ts, TSRK);
86: TSSetRHSFunction(ts, NULL, RHSFunction, &appctx);
87: TSSetRHSJacobian(ts, NULL, NULL, RHSJacobian, &appctx);
88: } else {
89: TSSetType(ts, TSCN);
90: TSSetIFunction(ts, NULL, IFunction, &appctx);
91: if (appctx.aijpc) {
92: Mat A, B;
94: DMSetMatType(da, MATSELL);
95: DMCreateMatrix(da, &A);
96: MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &B);
97: /* FIXME do we need to change viewer to display matrix in natural ordering as DMCreateMatrix_DA does? */
98: TSSetIJacobian(ts, A, B, IJacobian, &appctx);
99: MatDestroy(&A);
100: MatDestroy(&B);
101: } else {
102: TSSetIJacobian(ts, NULL, NULL, IJacobian, &appctx);
103: }
104: }
106: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107: Set initial conditions
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: InitialConditions(da, x);
110: TSSetSolution(ts, x);
112: /*
113: Have the TS save its trajectory so that TSAdjointSolve() may be used
114: */
115: if (!forwardonly) TSSetSaveTrajectory(ts);
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Set solver options
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: TSSetMaxTime(ts, 200.0);
121: TSSetTimeStep(ts, 0.5);
122: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
123: TSSetFromOptions(ts);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Solve ODE system
127: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128: TSSolve(ts, x);
129: if (!forwardonly) {
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Start the Adjoint model
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: VecDuplicate(x, &lambda[0]);
134: /* Reset initial conditions for the adjoint integration */
135: InitializeLambda(da, lambda[0], 0.5, 0.5);
136: TSSetCostGradients(ts, 1, lambda, NULL);
137: TSAdjointSolve(ts);
138: VecDestroy(&lambda[0]);
139: }
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Free work space. All PETSc objects should be destroyed when they
142: are no longer needed.
143: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: VecDestroy(&x);
145: TSDestroy(&ts);
146: DMDestroy(&da);
147: PetscFinalize();
148: return 0;
149: }
151: /* ------------------------------------------------------------------- */
152: PetscErrorCode InitialConditions(DM da, Vec U)
153: {
154: PetscInt i, j, xs, ys, xm, ym, Mx, My;
155: Field **u;
156: PetscReal hx, hy, x, y;
158: 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);
160: hx = 2.5 / (PetscReal)Mx;
161: hy = 2.5 / (PetscReal)My;
163: /*
164: Get pointers to vector data
165: */
166: DMDAVecGetArray(da, U, &u);
168: /*
169: Get local grid boundaries
170: */
171: DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
173: /*
174: Compute function over the locally owned part of the grid
175: */
176: for (j = ys; j < ys + ym; j++) {
177: y = j * hy;
178: for (i = xs; i < xs + xm; i++) {
179: x = i * hx;
180: if (PetscApproximateGTE(x, 1.0) && PetscApproximateLTE(x, 1.5) && PetscApproximateGTE(y, 1.0) && PetscApproximateLTE(y, 1.5))
181: u[j][i].v = PetscPowReal(PetscSinReal(4.0 * PETSC_PI * x), 2.0) * PetscPowReal(PetscSinReal(4.0 * PETSC_PI * y), 2.0) / 4.0;
182: else u[j][i].v = 0.0;
184: u[j][i].u = 1.0 - 2.0 * u[j][i].v;
185: }
186: }
188: /*
189: Restore vectors
190: */
191: DMDAVecRestoreArray(da, U, &u);
192: return 0;
193: }
195: /*TEST
197: build:
198: depends: reaction_diffusion.c
199: requires: !complex !single
201: test:
202: args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -da_grid_x 20 -da_grid_y 20
203: output_file: output/ex5adj_1.out
205: test:
206: suffix: 2
207: nsize: 2
208: args: -ts_max_steps 10 -ts_dt 10 -ts_monitor -ts_adjoint_monitor -ksp_monitor_short -da_grid_x 20 -da_grid_y 20 -ts_trajectory_dirname Test-dir -ts_trajectory_file_template test-%06D.cp
210: test:
211: suffix: 3
212: nsize: 2
213: args: -ts_max_steps 10 -ts_dt 10 -ts_adjoint_monitor_draw_sensi
215: test:
216: suffix: 4
217: nsize: 2
218: args: -ts_max_steps 10 -ts_dt 10 -ts_monitor -ts_adjoint_monitor -ksp_monitor_short -da_grid_x 20 -da_grid_y 20 -snes_fd_color
219: output_file: output/ex5adj_2.out
221: test:
222: suffix: 5
223: nsize: 2
224: args: -ts_max_steps 10 -implicitform 0 -ts_type rk -ts_rk_type 4 -ts_monitor -ts_adjoint_monitor -da_grid_x 20 -da_grid_y 20 -snes_fd_color
225: output_file: output/ex5adj_1.out
227: test:
228: suffix: knl
229: args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ts_trajectory_solution_only 0 -malloc_hbw -ts_trajectory_use_dram 1
230: output_file: output/ex5adj_3.out
231: requires: knl
233: test:
234: suffix: sell
235: nsize: 4
236: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type none
237: output_file: output/ex5adj_sell_1.out
239: test:
240: suffix: aijsell
241: nsize: 4
242: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell -pc_type none
243: output_file: output/ex5adj_sell_1.out
245: test:
246: suffix: sell2
247: nsize: 4
248: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type sor
249: output_file: output/ex5adj_sell_2.out
251: test:
252: suffix: aijsell2
253: nsize: 4
254: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type sor
255: output_file: output/ex5adj_sell_2.out
257: test:
258: suffix: sell3
259: nsize: 4
260: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type bjacobi -mg_levels_pc_type bjacobi
261: output_file: output/ex5adj_sell_3.out
263: test:
264: suffix: sell4
265: nsize: 4
266: args: -forwardonly -implicitform -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type bjacobi -mg_levels_pc_type bjacobi
267: output_file: output/ex5adj_sell_4.out
269: test:
270: suffix: sell5
271: nsize: 4
272: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -aijpc
273: output_file: output/ex5adj_sell_5.out
275: test:
276: suffix: aijsell5
277: nsize: 4
278: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell
279: output_file: output/ex5adj_sell_5.out
281: test:
282: suffix: sell6
283: args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ts_trajectory_solution_only 0 -dm_mat_type sell -pc_type jacobi
284: output_file: output/ex5adj_sell_6.out
286: test:
287: suffix: gamg1
288: args: -pc_type gamg -pc_gamg_esteig_ksp_type gmres -pc_gamg_esteig_ksp_max_it 10 -pc_mg_levels 2 -ksp_monitor_short -ts_max_steps 2 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ksp_rtol 1e-2 -pc_gamg_use_sa_esteig 0
289: output_file: output/ex5adj_gamg.out
291: test:
292: suffix: gamg2
293: args: -pc_type gamg -pc_gamg_esteig_ksp_type gmres -pc_gamg_esteig_ksp_max_it 10 -pc_mg_levels 2 -ksp_monitor_short -ts_max_steps 2 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ksp_use_explicittranspose -ksp_rtol 1e-2 -pc_gamg_use_sa_esteig 0
294: output_file: output/ex5adj_gamg.out
295: TEST*/