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/TAO interface to solve an inverse initial value problem built on a system of 6 time-dependent partial differential equations. 7 In this problem, the initial value for the PDE is unknown, but the output (the final solution of the PDE) is known. 8 We want to determine the initial value that can produce the given output. 9 We formulate the problem as a nonlinear optimization problem that minimizes the discrepency between the simulated 10 result and given reference solution, calculate the gradient of the objective function with the discrete adjoint 11 solver, and solve the optimization problem with TAO. 12 13 Runtime options: 14 -forwardonly - run only the forward simulation 15 -implicitform - provide IFunction and IJacobian to TS, if not set, RHSFunction and RHSJacobian will be used 16 */ 17 18 #include <petscsys.h> 19 #include <petscdm.h> 20 #include <petscdmda.h> 21 #include <petscts.h> 22 #include <petsctao.h> 23 24 typedef struct { 25 PetscScalar u,v; 26 } Field; 27 28 typedef struct { 29 PetscReal D1,D2,gamma,kappa; 30 TS ts; 31 Vec U; 32 } AppCtx; 33 34 /* User-defined routines */ 35 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*),InitialConditions(DM,Vec); 36 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*); 37 extern PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*); 38 extern PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*); 39 extern PetscErrorCode FormFunctionAndGradient(Tao,Vec,PetscReal*,Vec,void*); 40 41 /* 42 Set terminal condition for the adjoint variable 43 */ 44 PetscErrorCode InitializeLambda(DM da,Vec lambda,Vec U,AppCtx *appctx) 45 { 46 char filename[PETSC_MAX_PATH_LEN]=""; 47 PetscViewer viewer; 48 Vec Uob; 49 PetscErrorCode ierr; 50 51 ierr = VecDuplicate(U,&Uob);CHKERRQ(ierr); 52 ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr); 53 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 54 ierr = VecLoad(Uob,viewer);CHKERRQ(ierr); 55 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 56 ierr = VecAYPX(Uob,-1.,U);CHKERRQ(ierr); 57 ierr = VecScale(Uob,2.0);CHKERRQ(ierr); 58 ierr = VecAXPY(lambda,1.,Uob);CHKERRQ(ierr); 59 ierr = VecDestroy(&Uob);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 /* 64 Set up a viewer that dumps data into a binary file 65 */ 66 PetscErrorCode OutputBIN(DM da, const char *filename, PetscViewer *viewer) 67 { 68 PetscErrorCode ierr; 69 70 ierr = PetscViewerCreate(PetscObjectComm((PetscObject)da),viewer);CHKERRQ(ierr); 71 ierr = PetscViewerSetType(*viewer,PETSCVIEWERBINARY);CHKERRQ(ierr); 72 ierr = PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);CHKERRQ(ierr); 73 ierr = PetscViewerFileSetName(*viewer,filename);CHKERRQ(ierr); 74 PetscFunctionReturn(0); 75 } 76 77 /* 78 Generate a reference solution and save it to a binary file 79 */ 80 PetscErrorCode GenerateOBs(TS ts,Vec U,AppCtx *appctx) 81 { 82 char filename[PETSC_MAX_PATH_LEN] = ""; 83 PetscViewer viewer; 84 DM da; 85 PetscErrorCode ierr; 86 87 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 88 ierr = TSSolve(ts,U);CHKERRQ(ierr); 89 ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr); 90 ierr = OutputBIN(da,filename,&viewer);CHKERRQ(ierr); 91 ierr = VecView(U,viewer);CHKERRQ(ierr); 92 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 93 PetscFunctionReturn(0); 94 } 95 96 PetscErrorCode InitialConditions(DM da,Vec U) 97 { 98 PetscErrorCode ierr; 99 PetscInt i,j,xs,ys,xm,ym,Mx,My; 100 Field **u; 101 PetscReal hx,hy,x,y; 102 103 PetscFunctionBegin; 104 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); 105 106 hx = 2.5/(PetscReal)Mx; 107 hy = 2.5/(PetscReal)My; 108 109 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 110 /* Get local grid boundaries */ 111 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 112 113 /* Compute function over the locally owned part of the grid */ 114 for (j=ys; j<ys+ym; j++) { 115 y = j*hy; 116 for (i=xs; i<xs+xm; i++) { 117 x = i*hx; 118 if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0); 119 else u[j][i].v = 0.0; 120 121 u[j][i].u = 1.0 - 2.0*u[j][i].v; 122 } 123 } 124 125 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode PerturbedInitialConditions(DM da,Vec U) 130 { 131 PetscErrorCode ierr; 132 PetscInt i,j,xs,ys,xm,ym,Mx,My; 133 Field **u; 134 PetscReal hx,hy,x,y; 135 136 PetscFunctionBegin; 137 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); 138 139 hx = 2.5/(PetscReal)Mx; 140 hy = 2.5/(PetscReal)My; 141 142 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 143 /* Get local grid boundaries */ 144 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 145 146 /* Compute function over the locally owned part of the grid */ 147 for (j=ys; j<ys+ym; j++) { 148 y = j*hy; 149 for (i=xs; i<xs+xm; i++) { 150 x = i*hx; 151 if ((1.0 <= x) && (x <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*0.5*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*0.5*y),2.0); 152 else u[j][i].v = 0.0; 153 154 u[j][i].u = 1.0 - 2.0*u[j][i].v; 155 } 156 } 157 158 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 PetscErrorCode PerturbedInitialConditions2(DM da,Vec U) 163 { 164 PetscErrorCode ierr; 165 PetscInt i,j,xs,ys,xm,ym,Mx,My; 166 Field **u; 167 PetscReal hx,hy,x,y; 168 169 PetscFunctionBegin; 170 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); 171 172 hx = 2.5/(PetscReal)Mx; 173 hy = 2.5/(PetscReal)My; 174 175 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 176 /* Get local grid boundaries */ 177 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 178 179 /* Compute function over the locally owned part of the grid */ 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 ((1.0 <= x-0.5) && (x-0.5 <= 1.5) && (1.0 <= y) && (y <= 1.5)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*(x-0.5)),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.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 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 192 PetscFunctionReturn(0); 193 } 194 195 PetscErrorCode PerturbedInitialConditions3(DM da,Vec U) 196 { 197 PetscErrorCode ierr; 198 PetscInt i,j,xs,ys,xm,ym,Mx,My; 199 Field **u; 200 PetscReal hx,hy,x,y; 201 202 PetscFunctionBegin; 203 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); 204 205 hx = 2.5/(PetscReal)Mx; 206 hy = 2.5/(PetscReal)My; 207 208 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr); 209 /* Get local grid boundaries */ 210 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 211 212 /* Compute function over the locally owned part of the grid */ 213 for (j=ys; j<ys+ym; j++) { 214 y = j*hy; 215 for (i=xs; i<xs+xm; i++) { 216 x = i*hx; 217 if ((0.5 <= x) && (x <= 2.0) && (0.5 <= y) && (y <= 2.0)) u[j][i].v = .25*PetscPowReal(PetscSinReal(4.0*PETSC_PI*x),2.0)*PetscPowReal(PetscSinReal(4.0*PETSC_PI*y),2.0); 218 else u[j][i].v = 0.0; 219 220 u[j][i].u = 1.0 - 2.0*u[j][i].v; 221 } 222 } 223 224 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr); 225 PetscFunctionReturn(0); 226 } 227 228 int main(int argc,char **argv) 229 { 230 PetscErrorCode ierr; 231 DM da; 232 AppCtx appctx; 233 PetscBool forwardonly = PETSC_FALSE,implicitform = PETSC_FALSE; 234 PetscInt perturbic = 1; 235 236 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 237 ierr = PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);CHKERRQ(ierr); 238 ierr = PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);CHKERRQ(ierr); 239 ierr = PetscOptionsGetInt(NULL,NULL,"-perturbic",&perturbic,NULL);CHKERRQ(ierr); 240 241 appctx.D1 = 8.0e-5; 242 appctx.D2 = 4.0e-5; 243 appctx.gamma = .024; 244 appctx.kappa = .06; 245 246 /* Create distributed array (DMDA) to manage parallel grid and vectors */ 247 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); 248 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 249 ierr = DMSetUp(da);CHKERRQ(ierr); 250 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr); 251 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr); 252 253 /* Extract global vectors from DMDA; then duplicate for remaining 254 vectors that are the same types */ 255 ierr = DMCreateGlobalVector(da,&appctx.U);CHKERRQ(ierr); 256 257 /* Create timestepping solver context */ 258 ierr = TSCreate(PETSC_COMM_WORLD,&appctx.ts);CHKERRQ(ierr); 259 ierr = TSSetType(appctx.ts,TSCN);CHKERRQ(ierr); 260 ierr = TSSetDM(appctx.ts,da);CHKERRQ(ierr); 261 ierr = TSSetProblemType(appctx.ts,TS_NONLINEAR);CHKERRQ(ierr); 262 ierr = TSSetEquationType(appctx.ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */ 263 if (!implicitform) { 264 ierr = TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);CHKERRQ(ierr); 265 ierr = TSSetRHSJacobian(appctx.ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr); 266 } else { 267 ierr = TSSetIFunction(appctx.ts,NULL,IFunction,&appctx);CHKERRQ(ierr); 268 ierr = TSSetIJacobian(appctx.ts,NULL,NULL,IJacobian,&appctx);CHKERRQ(ierr); 269 } 270 271 /* Set initial conditions */ 272 ierr = InitialConditions(da,appctx.U);CHKERRQ(ierr); 273 ierr = TSSetSolution(appctx.ts,appctx.U);CHKERRQ(ierr); 274 275 /* Set solver options */ 276 ierr = TSSetMaxTime(appctx.ts,2000.0);CHKERRQ(ierr); 277 ierr = TSSetTimeStep(appctx.ts,0.5);CHKERRQ(ierr); 278 ierr = TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr); 279 ierr = TSSetFromOptions(appctx.ts);CHKERRQ(ierr); 280 281 ierr = GenerateOBs(appctx.ts,appctx.U,&appctx);CHKERRQ(ierr); 282 283 if (!forwardonly) { 284 Tao tao; 285 Vec P; 286 Vec lambda[1]; 287 PetscLogStage opt_stage; 288 289 ierr = PetscLogStageRegister("Optimization",&opt_stage);CHKERRQ(ierr); 290 ierr = PetscLogStagePush(opt_stage);CHKERRQ(ierr); 291 if (perturbic == 1) { 292 ierr = PerturbedInitialConditions(da,appctx.U);CHKERRQ(ierr); 293 } else if (perturbic == 2) { 294 ierr = PerturbedInitialConditions2(da,appctx.U);CHKERRQ(ierr); 295 } else if (perturbic == 3) { 296 ierr = PerturbedInitialConditions3(da,appctx.U);CHKERRQ(ierr); 297 } 298 299 ierr = VecDuplicate(appctx.U,&lambda[0]);CHKERRQ(ierr); 300 ierr = TSSetCostGradients(appctx.ts,1,lambda,NULL);CHKERRQ(ierr); 301 302 /* Have the TS save its trajectory needed by TSAdjointSolve() */ 303 ierr = TSSetSaveTrajectory(appctx.ts);CHKERRQ(ierr); 304 305 /* Create TAO solver and set desired solution method */ 306 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 307 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr); 308 309 /* Set initial guess for TAO */ 310 ierr = VecDuplicate(appctx.U,&P);CHKERRQ(ierr); 311 ierr = VecCopy(appctx.U,P);CHKERRQ(ierr); 312 ierr = TaoSetInitialVector(tao,P);CHKERRQ(ierr); 313 314 /* Set routine for function and gradient evaluation */ 315 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);CHKERRQ(ierr); 316 317 /* Check for any TAO command line options */ 318 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 319 320 ierr = TaoSolve(tao);CHKERRQ(ierr); 321 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 322 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr); 323 ierr = VecDestroy(&P);CHKERRQ(ierr); 324 ierr = PetscLogStagePop();CHKERRQ(ierr); 325 } 326 327 /* Free work space. All PETSc objects should be destroyed when they 328 are no longer needed. */ 329 ierr = VecDestroy(&appctx.U);CHKERRQ(ierr); 330 ierr = TSDestroy(&appctx.ts);CHKERRQ(ierr); 331 ierr = DMDestroy(&da);CHKERRQ(ierr); 332 ierr = PetscFinalize(); 333 return ierr; 334 } 335 336 /* ------------------------ TS callbacks ---------------------------- */ 337 338 /* 339 RHSFunction - Evaluates nonlinear function, F(x). 340 341 Input Parameters: 342 . ts - the TS context 343 . X - input vector 344 . ptr - optional user-defined context, as set by TSSetRHSFunction() 345 346 Output Parameter: 347 . F - function vector 348 */ 349 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr) 350 { 351 AppCtx *appctx = (AppCtx*)ptr; 352 DM da; 353 PetscErrorCode ierr; 354 PetscInt i,j,Mx,My,xs,ys,xm,ym; 355 PetscReal hx,hy,sx,sy; 356 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 357 Field **u,**f; 358 Vec localU; 359 360 PetscFunctionBegin; 361 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 362 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 363 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); 364 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 365 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 366 367 /* Scatter ghost points to local vector,using the 2-step process 368 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 369 By placing code between these two statements, computations can be 370 done while messages are in transition. */ 371 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 372 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 373 374 /* Get pointers to vector data */ 375 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 376 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 377 378 /* Get local grid boundaries */ 379 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 380 381 /* Compute function over the locally owned part of the grid */ 382 for (j=ys; j<ys+ym; j++) { 383 for (i=xs; i<xs+xm; i++) { 384 uc = u[j][i].u; 385 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 386 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 387 vc = u[j][i].v; 388 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 389 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 390 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); 391 f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; 392 } 393 } 394 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 395 396 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 397 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 398 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 399 PetscFunctionReturn(0); 400 } 401 402 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx) 403 { 404 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 405 DM da; 406 PetscErrorCode ierr; 407 PetscInt i,j,Mx,My,xs,ys,xm,ym; 408 PetscReal hx,hy,sx,sy; 409 PetscScalar uc,vc; 410 Field **u; 411 Vec localU; 412 MatStencil stencil[6],rowstencil; 413 PetscScalar entries[6]; 414 415 PetscFunctionBegin; 416 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 417 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 418 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); 419 420 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 421 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 422 423 /* Scatter ghost points to local vector,using the 2-step process 424 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 425 By placing code between these two statements, computations can be 426 done while messages are in transition. */ 427 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 428 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 429 430 /* Get pointers to vector data */ 431 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 432 433 /* Get local grid boundaries */ 434 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 435 436 stencil[0].k = 0; 437 stencil[1].k = 0; 438 stencil[2].k = 0; 439 stencil[3].k = 0; 440 stencil[4].k = 0; 441 stencil[5].k = 0; 442 rowstencil.k = 0; 443 444 /* Compute function over the locally owned part of the grid */ 445 for (j=ys; j<ys+ym; j++) { 446 stencil[0].j = j-1; 447 stencil[1].j = j+1; 448 stencil[2].j = j; 449 stencil[3].j = j; 450 stencil[4].j = j; 451 stencil[5].j = j; 452 rowstencil.k = 0; rowstencil.j = j; 453 for (i=xs; i<xs+xm; i++) { 454 uc = u[j][i].u; 455 vc = u[j][i].v; 456 457 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 458 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 459 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 460 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 461 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ 462 463 stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy; 464 stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy; 465 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx; 466 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx; 467 stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma; 468 stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc; 469 rowstencil.i = i; rowstencil.c = 0; 470 471 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 472 stencil[0].c = 1; entries[0] = appctx->D2*sy; 473 stencil[1].c = 1; entries[1] = appctx->D2*sy; 474 stencil[2].c = 1; entries[2] = appctx->D2*sx; 475 stencil[3].c = 1; entries[3] = appctx->D2*sx; 476 stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa; 477 stencil[5].c = 0; entries[5] = vc*vc; 478 rowstencil.c = 1; 479 480 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 481 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 482 } 483 } 484 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 485 486 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 487 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 488 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 489 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 490 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 491 PetscFunctionReturn(0); 492 } 493 494 /* 495 IFunction - Evaluates implicit nonlinear function, xdot - F(x). 496 497 Input Parameters: 498 . ts - the TS context 499 . U - input vector 500 . Udot - input vector 501 . ptr - optional user-defined context, as set by TSSetRHSFunction() 502 503 Output Parameter: 504 . F - function vector 505 */ 506 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 507 { 508 AppCtx *appctx = (AppCtx*)ptr; 509 DM da; 510 PetscErrorCode ierr; 511 PetscInt i,j,Mx,My,xs,ys,xm,ym; 512 PetscReal hx,hy,sx,sy; 513 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 514 Field **u,**f,**udot; 515 Vec localU; 516 517 PetscFunctionBegin; 518 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 519 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 520 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); 521 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 522 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 523 524 /* Scatter ghost points to local vector,using the 2-step process 525 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 526 By placing code between these two statements, computations can be 527 done while messages are in transition. */ 528 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 529 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 530 531 /* Get pointers to vector data */ 532 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 533 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 534 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 535 536 /* Get local grid boundaries */ 537 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 538 539 /* Compute function over the locally owned part of the grid */ 540 for (j=ys; j<ys+ym; j++) { 541 for (i=xs; i<xs+xm; i++) { 542 uc = u[j][i].u; 543 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 544 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 545 vc = u[j][i].v; 546 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 547 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 548 f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc)); 549 f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc); 550 } 551 } 552 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 553 554 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 555 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 556 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 557 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 558 PetscFunctionReturn(0); 559 } 560 561 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx) 562 { 563 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 564 DM da; 565 PetscErrorCode ierr; 566 PetscInt i,j,Mx,My,xs,ys,xm,ym; 567 PetscReal hx,hy,sx,sy; 568 PetscScalar uc,vc; 569 Field **u; 570 Vec localU; 571 MatStencil stencil[6],rowstencil; 572 PetscScalar entries[6]; 573 574 PetscFunctionBegin; 575 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 576 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 577 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); 578 579 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 580 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 581 582 /* Scatter ghost points to local vector,using the 2-step process 583 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 584 By placing code between these two statements, computations can be 585 done while messages are in transition.*/ 586 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 587 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 588 589 /* Get pointers to vector data */ 590 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 591 592 /* Get local grid boundaries */ 593 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 594 595 stencil[0].k = 0; 596 stencil[1].k = 0; 597 stencil[2].k = 0; 598 stencil[3].k = 0; 599 stencil[4].k = 0; 600 stencil[5].k = 0; 601 rowstencil.k = 0; 602 603 /* Compute function over the locally owned part of the grid */ 604 for (j=ys; j<ys+ym; j++) { 605 stencil[0].j = j-1; 606 stencil[1].j = j+1; 607 stencil[2].j = j; 608 stencil[3].j = j; 609 stencil[4].j = j; 610 stencil[5].j = j; 611 rowstencil.k = 0; rowstencil.j = j; 612 for (i=xs; i<xs+xm; i++) { 613 uc = u[j][i].u; 614 vc = u[j][i].v; 615 616 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 617 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 618 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 619 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 620 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); */ 621 stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy; 622 stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy; 623 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx; 624 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx; 625 stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a; 626 stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc; 627 rowstencil.i = i; rowstencil.c = 0; 628 629 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 630 stencil[0].c = 1; entries[0] = -appctx->D2*sy; 631 stencil[1].c = 1; entries[1] = -appctx->D2*sy; 632 stencil[2].c = 1; entries[2] = -appctx->D2*sx; 633 stencil[3].c = 1; entries[3] = -appctx->D2*sx; 634 stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a; 635 stencil[5].c = 0; entries[5] = -vc*vc; 636 rowstencil.c = 1; 637 638 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 639 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 640 } 641 } 642 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 643 644 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 645 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 646 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 647 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 648 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 649 PetscFunctionReturn(0); 650 } 651 652 /* ------------------------ TAO callbacks ---------------------------- */ 653 654 /* 655 FormFunctionAndGradient - Evaluates the function and corresponding gradient. 656 657 Input Parameters: 658 tao - the Tao context 659 P - the input vector 660 ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine() 661 662 Output Parameters: 663 f - the newly evaluated function 664 G - the newly evaluated gradient 665 */ 666 PetscErrorCode FormFunctionAndGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx) 667 { 668 AppCtx *appctx = (AppCtx*)ctx; 669 PetscReal soberr,timestep; 670 Vec *lambda; 671 Vec SDiff; 672 DM da; 673 char filename[PETSC_MAX_PATH_LEN]=""; 674 PetscViewer viewer; 675 PetscErrorCode ierr; 676 677 PetscFunctionBeginUser; 678 ierr = TSSetTime(appctx->ts,0.0);CHKERRQ(ierr); 679 ierr = TSGetTimeStep(appctx->ts,×tep);CHKERRQ(ierr); 680 if (timestep<0) { 681 ierr = TSSetTimeStep(appctx->ts,-timestep);CHKERRQ(ierr); 682 } 683 ierr = TSSetStepNumber(appctx->ts,0);CHKERRQ(ierr); 684 ierr = TSSetFromOptions(appctx->ts);CHKERRQ(ierr); 685 686 ierr = VecDuplicate(P,&SDiff);CHKERRQ(ierr); 687 ierr = VecCopy(P,appctx->U);CHKERRQ(ierr); 688 ierr = TSGetDM(appctx->ts,&da);CHKERRQ(ierr); 689 *f = 0; 690 691 ierr = TSSolve(appctx->ts,appctx->U);CHKERRQ(ierr); 692 ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr); 693 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 694 ierr = VecLoad(SDiff,viewer);CHKERRQ(ierr); 695 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 696 ierr = VecAYPX(SDiff,-1.,appctx->U);CHKERRQ(ierr); 697 ierr = VecDot(SDiff,SDiff,&soberr);CHKERRQ(ierr); 698 *f += soberr; 699 700 ierr = TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);CHKERRQ(ierr); 701 ierr = VecSet(lambda[0],0.0);CHKERRQ(ierr); 702 ierr = InitializeLambda(da,lambda[0],appctx->U,appctx);CHKERRQ(ierr); 703 ierr = TSAdjointSolve(appctx->ts);CHKERRQ(ierr); 704 705 ierr = VecCopy(lambda[0],G);CHKERRQ(ierr); 706 707 ierr = VecDestroy(&SDiff);CHKERRQ(ierr); 708 PetscFunctionReturn(0); 709 } 710 711 /*TEST 712 713 build: 714 requires: !complex !single 715 716 test: 717 args: -ts_max_steps 5 -ts_type rk -ts_rk_type 3 -ts_trajectory_type memory -tao_monitor -tao_view -tao_gatol 1e-6 718 output_file: output/ex5opt_ic_1.out 719 720 TEST*/ 721