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;
50: PetscInitialize(&argc,&argv,(char*)0,help);
51: PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);
52: PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);
53: appctx.aijpc = PETSC_FALSE;
54: PetscOptionsGetBool(NULL,NULL,"-aijpc",&appctx.aijpc,NULL);
56: appctx.D1 = 8.0e-5;
57: appctx.D2 = 4.0e-5;
58: appctx.gamma = .024;
59: appctx.kappa = .06;
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Create distributed array (DMDA) to manage parallel grid and vectors
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,64,64,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
65: DMSetFromOptions(da);
66: DMSetUp(da);
67: DMDASetFieldName(da,0,"u");
68: DMDASetFieldName(da,1,"v");
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Extract global vectors from DMDA; then duplicate for remaining
72: vectors that are the same types
73: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
74: DMCreateGlobalVector(da,&x);
76: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77: Create timestepping solver context
78: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: TSCreate(PETSC_COMM_WORLD,&ts);
80: TSSetDM(ts,da);
81: TSSetProblemType(ts,TS_NONLINEAR);
82: TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
83: if (!implicitform) {
84: TSSetType(ts,TSRK);
85: TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);
86: TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);
87: } else {
88: TSSetType(ts,TSCN);
89: TSSetIFunction(ts,NULL,IFunction,&appctx);
90: if (appctx.aijpc) {
91: Mat A,B;
93: DMSetMatType(da,MATSELL);
94: DMCreateMatrix(da,&A);
95: MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);
96: /* FIXME do we need to change viewer to display matrix in natural ordering as DMCreateMatrix_DA does? */
97: TSSetIJacobian(ts,A,B,IJacobian,&appctx);
98: MatDestroy(&A);
99: MatDestroy(&B);
100: } else {
101: TSSetIJacobian(ts,NULL,NULL,IJacobian,&appctx);
102: }
103: }
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Set initial conditions
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: InitialConditions(da,x);
109: TSSetSolution(ts,x);
111: /*
112: Have the TS save its trajectory so that TSAdjointSolve() may be used
113: */
114: if (!forwardonly) TSSetSaveTrajectory(ts);
116: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: Set solver options
118: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: TSSetMaxTime(ts,200.0);
120: TSSetTimeStep(ts,0.5);
121: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);
122: TSSetFromOptions(ts);
124: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: Solve ODE system
126: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127: TSSolve(ts,x);
128: if (!forwardonly) {
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Start the Adjoint model
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: VecDuplicate(x,&lambda[0]);
133: /* Reset initial conditions for the adjoint integration */
134: InitializeLambda(da,lambda[0],0.5,0.5);
135: TSSetCostGradients(ts,1,lambda,NULL);
136: TSAdjointSolve(ts);
137: VecDestroy(&lambda[0]);
138: }
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Free work space. All PETSc objects should be destroyed when they
141: are no longer needed.
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: VecDestroy(&x);
144: TSDestroy(&ts);
145: DMDestroy(&da);
146: PetscFinalize();
147: return 0;
148: }
150: /* ------------------------------------------------------------------- */
151: PetscErrorCode InitialConditions(DM da,Vec U)
152: {
153: PetscInt i,j,xs,ys,xm,ym,Mx,My;
154: Field **u;
155: PetscReal hx,hy,x,y;
157: 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);
159: hx = 2.5/(PetscReal)Mx;
160: hy = 2.5/(PetscReal)My;
162: /*
163: Get pointers to vector data
164: */
165: DMDAVecGetArray(da,U,&u);
167: /*
168: Get local grid boundaries
169: */
170: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
172: /*
173: Compute function over the locally owned part of the grid
174: */
175: for (j=ys; j<ys+ym; j++) {
176: y = j*hy;
177: for (i=xs; i<xs+xm; i++) {
178: x = i*hx;
179: if (PetscApproximateGTE(x,1.0) && PetscApproximateLTE(x,1.5) && PetscApproximateGTE(y,1.0) && PetscApproximateLTE(y,1.5)) u[j][i].v = PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0)/4.0;
180: else u[j][i].v = 0.0;
182: u[j][i].u = 1.0 - 2.0*u[j][i].v;
183: }
184: }
186: /*
187: Restore vectors
188: */
189: DMDAVecRestoreArray(da,U,&u);
190: return 0;
191: }
193: /*TEST
195: build:
196: depends: reaction_diffusion.c
197: requires: !complex !single
199: test:
200: args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -da_grid_x 20 -da_grid_y 20
201: output_file: output/ex5adj_1.out
203: test:
204: suffix: 2
205: nsize: 2
206: 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
208: test:
209: suffix: 3
210: nsize: 2
211: args: -ts_max_steps 10 -ts_dt 10 -ts_adjoint_monitor_draw_sensi
213: test:
214: suffix: 4
215: nsize: 2
216: 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
217: output_file: output/ex5adj_2.out
219: test:
220: suffix: 5
221: nsize: 2
222: 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
223: output_file: output/ex5adj_1.out
225: test:
226: suffix: knl
227: 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
228: output_file: output/ex5adj_3.out
229: requires: knl
231: test:
232: suffix: sell
233: nsize: 4
234: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type none
235: output_file: output/ex5adj_sell_1.out
237: test:
238: suffix: aijsell
239: nsize: 4
240: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell -pc_type none
241: output_file: output/ex5adj_sell_1.out
243: test:
244: suffix: sell2
245: nsize: 4
246: 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
247: output_file: output/ex5adj_sell_2.out
249: test:
250: suffix: aijsell2
251: nsize: 4
252: 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
253: output_file: output/ex5adj_sell_2.out
255: test:
256: suffix: sell3
257: nsize: 4
258: 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
259: output_file: output/ex5adj_sell_3.out
261: test:
262: suffix: sell4
263: nsize: 4
264: 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
265: output_file: output/ex5adj_sell_4.out
267: test:
268: suffix: sell5
269: nsize: 4
270: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -aijpc
271: output_file: output/ex5adj_sell_5.out
273: test:
274: suffix: aijsell5
275: nsize: 4
276: args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell
277: output_file: output/ex5adj_sell_5.out
279: test:
280: suffix: sell6
281: 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
282: output_file: output/ex5adj_sell_6.out
284: test:
285: suffix: gamg1
286: 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
287: output_file: output/ex5adj_gamg.out
289: test:
290: suffix: gamg2
291: 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
292: output_file: output/ex5adj_gamg.out
293: TEST*/