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 #if defined(PETSC_USE_LOG) 288 PetscLogStage opt_stage; 289 #endif 290 291 ierr = PetscLogStageRegister("Optimization",&opt_stage);CHKERRQ(ierr); 292 ierr = PetscLogStagePush(opt_stage);CHKERRQ(ierr); 293 if (perturbic == 1) { 294 ierr = PerturbedInitialConditions(da,appctx.U);CHKERRQ(ierr); 295 } else if (perturbic == 2) { 296 ierr = PerturbedInitialConditions2(da,appctx.U);CHKERRQ(ierr); 297 } else if (perturbic == 3) { 298 ierr = PerturbedInitialConditions3(da,appctx.U);CHKERRQ(ierr); 299 } 300 301 ierr = VecDuplicate(appctx.U,&lambda[0]);CHKERRQ(ierr); 302 ierr = TSSetCostGradients(appctx.ts,1,lambda,NULL);CHKERRQ(ierr); 303 304 /* Have the TS save its trajectory needed by TSAdjointSolve() */ 305 ierr = TSSetSaveTrajectory(appctx.ts);CHKERRQ(ierr); 306 307 /* Create TAO solver and set desired solution method */ 308 ierr = TaoCreate(PETSC_COMM_WORLD,&tao);CHKERRQ(ierr); 309 ierr = TaoSetType(tao,TAOBLMVM);CHKERRQ(ierr); 310 311 /* Set initial guess for TAO */ 312 ierr = VecDuplicate(appctx.U,&P);CHKERRQ(ierr); 313 ierr = VecCopy(appctx.U,P);CHKERRQ(ierr); 314 ierr = TaoSetInitialVector(tao,P);CHKERRQ(ierr); 315 316 /* Set routine for function and gradient evaluation */ 317 ierr = TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);CHKERRQ(ierr); 318 319 /* Check for any TAO command line options */ 320 ierr = TaoSetFromOptions(tao);CHKERRQ(ierr); 321 322 ierr = TaoSolve(tao);CHKERRQ(ierr); 323 ierr = TaoDestroy(&tao);CHKERRQ(ierr); 324 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr); 325 ierr = VecDestroy(&P);CHKERRQ(ierr); 326 ierr = PetscLogStagePop();CHKERRQ(ierr); 327 } 328 329 /* Free work space. All PETSc objects should be destroyed when they 330 are no longer needed. */ 331 ierr = VecDestroy(&appctx.U);CHKERRQ(ierr); 332 ierr = TSDestroy(&appctx.ts);CHKERRQ(ierr); 333 ierr = DMDestroy(&da);CHKERRQ(ierr); 334 ierr = PetscFinalize(); 335 return ierr; 336 } 337 338 /* ------------------------ TS callbacks ---------------------------- */ 339 340 /* 341 RHSFunction - Evaluates nonlinear function, F(x). 342 343 Input Parameters: 344 . ts - the TS context 345 . X - input vector 346 . ptr - optional user-defined context, as set by TSSetRHSFunction() 347 348 Output Parameter: 349 . F - function vector 350 */ 351 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr) 352 { 353 AppCtx *appctx = (AppCtx*)ptr; 354 DM da; 355 PetscErrorCode ierr; 356 PetscInt i,j,Mx,My,xs,ys,xm,ym; 357 PetscReal hx,hy,sx,sy; 358 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 359 Field **u,**f; 360 Vec localU; 361 362 PetscFunctionBegin; 363 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 364 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 365 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); 366 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 367 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 368 369 /* Scatter ghost points to local vector,using the 2-step process 370 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 371 By placing code between these two statements, computations can be 372 done while messages are in transition. */ 373 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 374 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 375 376 /* Get pointers to vector data */ 377 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 378 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 379 380 /* Get local grid boundaries */ 381 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 382 383 /* Compute function over the locally owned part of the grid */ 384 for (j=ys; j<ys+ym; j++) { 385 for (i=xs; i<xs+xm; i++) { 386 uc = u[j][i].u; 387 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 388 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 389 vc = u[j][i].v; 390 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 391 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 392 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); 393 f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; 394 } 395 } 396 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 397 398 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 399 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 400 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx) 405 { 406 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 407 DM da; 408 PetscErrorCode ierr; 409 PetscInt i,j,Mx,My,xs,ys,xm,ym; 410 PetscReal hx,hy,sx,sy; 411 PetscScalar uc,vc; 412 Field **u; 413 Vec localU; 414 MatStencil stencil[6],rowstencil; 415 PetscScalar entries[6]; 416 417 PetscFunctionBegin; 418 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 419 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 420 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); 421 422 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 423 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 424 425 /* Scatter ghost points to local vector,using the 2-step process 426 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 427 By placing code between these two statements, computations can be 428 done while messages are in transition. */ 429 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 430 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 431 432 /* Get pointers to vector data */ 433 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 434 435 /* Get local grid boundaries */ 436 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 437 438 stencil[0].k = 0; 439 stencil[1].k = 0; 440 stencil[2].k = 0; 441 stencil[3].k = 0; 442 stencil[4].k = 0; 443 stencil[5].k = 0; 444 rowstencil.k = 0; 445 446 /* Compute function over the locally owned part of the grid */ 447 for (j=ys; j<ys+ym; j++) { 448 stencil[0].j = j-1; 449 stencil[1].j = j+1; 450 stencil[2].j = j; 451 stencil[3].j = j; 452 stencil[4].j = j; 453 stencil[5].j = j; 454 rowstencil.k = 0; rowstencil.j = j; 455 for (i=xs; i<xs+xm; i++) { 456 uc = u[j][i].u; 457 vc = u[j][i].v; 458 459 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 460 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 461 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 462 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 463 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/ 464 465 stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy; 466 stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy; 467 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx; 468 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx; 469 stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma; 470 stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc; 471 rowstencil.i = i; rowstencil.c = 0; 472 473 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 474 stencil[0].c = 1; entries[0] = appctx->D2*sy; 475 stencil[1].c = 1; entries[1] = appctx->D2*sy; 476 stencil[2].c = 1; entries[2] = appctx->D2*sx; 477 stencil[3].c = 1; entries[3] = appctx->D2*sx; 478 stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa; 479 stencil[5].c = 0; entries[5] = vc*vc; 480 rowstencil.c = 1; 481 482 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 483 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 484 } 485 } 486 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 487 488 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 489 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 490 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 491 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 492 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 493 PetscFunctionReturn(0); 494 } 495 496 /* 497 IFunction - Evaluates implicit nonlinear function, xdot - F(x). 498 499 Input Parameters: 500 . ts - the TS context 501 . U - input vector 502 . Udot - input vector 503 . ptr - optional user-defined context, as set by TSSetRHSFunction() 504 505 Output Parameter: 506 . F - function vector 507 */ 508 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr) 509 { 510 AppCtx *appctx = (AppCtx*)ptr; 511 DM da; 512 PetscErrorCode ierr; 513 PetscInt i,j,Mx,My,xs,ys,xm,ym; 514 PetscReal hx,hy,sx,sy; 515 PetscScalar uc,uxx,uyy,vc,vxx,vyy; 516 Field **u,**f,**udot; 517 Vec localU; 518 519 PetscFunctionBegin; 520 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 521 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 522 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); 523 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 524 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 525 526 /* Scatter ghost points to local vector,using the 2-step process 527 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 528 By placing code between these two statements, computations can be 529 done while messages are in transition. */ 530 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 531 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 532 533 /* Get pointers to vector data */ 534 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 535 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 536 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr); 537 538 /* Get local grid boundaries */ 539 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 540 541 /* Compute function over the locally owned part of the grid */ 542 for (j=ys; j<ys+ym; j++) { 543 for (i=xs; i<xs+xm; i++) { 544 uc = u[j][i].u; 545 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 546 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 547 vc = u[j][i].v; 548 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 549 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 550 f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc)); 551 f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc); 552 } 553 } 554 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr); 555 556 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 557 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 558 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr); 559 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 560 PetscFunctionReturn(0); 561 } 562 563 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx) 564 { 565 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */ 566 DM da; 567 PetscErrorCode ierr; 568 PetscInt i,j,Mx,My,xs,ys,xm,ym; 569 PetscReal hx,hy,sx,sy; 570 PetscScalar uc,vc; 571 Field **u; 572 Vec localU; 573 MatStencil stencil[6],rowstencil; 574 PetscScalar entries[6]; 575 576 PetscFunctionBegin; 577 ierr = TSGetDM(ts,&da);CHKERRQ(ierr); 578 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr); 579 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); 580 581 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx); 582 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy); 583 584 /* Scatter ghost points to local vector,using the 2-step process 585 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 586 By placing code between these two statements, computations can be 587 done while messages are in transition.*/ 588 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 589 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr); 590 591 /* Get pointers to vector data */ 592 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr); 593 594 /* Get local grid boundaries */ 595 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 596 597 stencil[0].k = 0; 598 stencil[1].k = 0; 599 stencil[2].k = 0; 600 stencil[3].k = 0; 601 stencil[4].k = 0; 602 stencil[5].k = 0; 603 rowstencil.k = 0; 604 605 /* Compute function over the locally owned part of the grid */ 606 for (j=ys; j<ys+ym; j++) { 607 stencil[0].j = j-1; 608 stencil[1].j = j+1; 609 stencil[2].j = j; 610 stencil[3].j = j; 611 stencil[4].j = j; 612 stencil[5].j = j; 613 rowstencil.k = 0; rowstencil.j = j; 614 for (i=xs; i<xs+xm; i++) { 615 uc = u[j][i].u; 616 vc = u[j][i].v; 617 618 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx; 619 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy; 620 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx; 621 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy; 622 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc); */ 623 stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy; 624 stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy; 625 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx; 626 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx; 627 stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a; 628 stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc; 629 rowstencil.i = i; rowstencil.c = 0; 630 631 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 632 stencil[0].c = 1; entries[0] = -appctx->D2*sy; 633 stencil[1].c = 1; entries[1] = -appctx->D2*sy; 634 stencil[2].c = 1; entries[2] = -appctx->D2*sx; 635 stencil[3].c = 1; entries[3] = -appctx->D2*sx; 636 stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a; 637 stencil[5].c = 0; entries[5] = -vc*vc; 638 rowstencil.c = 1; 639 640 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr); 641 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */ 642 } 643 } 644 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr); 645 646 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr); 647 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr); 648 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 649 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 650 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 651 PetscFunctionReturn(0); 652 } 653 654 /* ------------------------ TAO callbacks ---------------------------- */ 655 656 /* 657 FormFunctionAndGradient - Evaluates the function and corresponding gradient. 658 659 Input Parameters: 660 tao - the Tao context 661 P - the input vector 662 ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine() 663 664 Output Parameters: 665 f - the newly evaluated function 666 G - the newly evaluated gradient 667 */ 668 PetscErrorCode FormFunctionAndGradient(Tao tao,Vec P,PetscReal *f,Vec G,void *ctx) 669 { 670 AppCtx *appctx = (AppCtx*)ctx; 671 PetscReal soberr,timestep; 672 Vec *lambda; 673 Vec SDiff; 674 DM da; 675 char filename[PETSC_MAX_PATH_LEN]=""; 676 PetscViewer viewer; 677 PetscErrorCode ierr; 678 679 PetscFunctionBeginUser; 680 ierr = TSSetTime(appctx->ts,0.0);CHKERRQ(ierr); 681 ierr = TSGetTimeStep(appctx->ts,×tep);CHKERRQ(ierr); 682 if (timestep<0) { 683 ierr = TSSetTimeStep(appctx->ts,-timestep);CHKERRQ(ierr); 684 } 685 ierr = TSSetStepNumber(appctx->ts,0);CHKERRQ(ierr); 686 ierr = TSSetFromOptions(appctx->ts);CHKERRQ(ierr); 687 688 ierr = VecDuplicate(P,&SDiff);CHKERRQ(ierr); 689 ierr = VecCopy(P,appctx->U);CHKERRQ(ierr); 690 ierr = TSGetDM(appctx->ts,&da);CHKERRQ(ierr); 691 *f = 0; 692 693 ierr = TSSolve(appctx->ts,appctx->U);CHKERRQ(ierr); 694 ierr = PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");CHKERRQ(ierr); 695 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);CHKERRQ(ierr); 696 ierr = VecLoad(SDiff,viewer);CHKERRQ(ierr); 697 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 698 ierr = VecAYPX(SDiff,-1.,appctx->U);CHKERRQ(ierr); 699 ierr = VecDot(SDiff,SDiff,&soberr);CHKERRQ(ierr); 700 *f += soberr; 701 702 ierr = TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);CHKERRQ(ierr); 703 ierr = VecSet(lambda[0],0.0);CHKERRQ(ierr); 704 ierr = InitializeLambda(da,lambda[0],appctx->U,appctx);CHKERRQ(ierr); 705 ierr = TSAdjointSolve(appctx->ts);CHKERRQ(ierr); 706 707 ierr = VecCopy(lambda[0],G);CHKERRQ(ierr); 708 709 ierr = VecDestroy(&SDiff);CHKERRQ(ierr); 710 PetscFunctionReturn(0); 711 } 712 713 /*TEST 714 715 build: 716 requires: !complex !single 717 718 test: 719 args: -ts_max_steps 5 -ts_type rk -ts_rk_type 3 -ts_trajectory_type memory -tao_monitor -tao_view -tao_gatol 1e-6 720 output_file: output/ex5opt_ic_1.out 721 722 TEST*/ 723