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