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