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