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