Actual source code: ex5opt_ic.c
petsc-3.12.0 2019-09-29
1: static char help[] = "Demonstrates adjoint sensitivity analysis for Reaction-Diffusion Equations.\n";
3: /*
4: See ex5.c for details on the equation.
5: This code demonestrates the TSAdjoint/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.
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: */
18: #include <petscsys.h>
19: #include <petscdm.h>
20: #include <petscdmda.h>
21: #include <petscts.h>
22: #include <petsctao.h>
24: typedef struct {
25: PetscScalar u,v;
26: } Field;
28: typedef struct {
29: PetscReal D1,D2,gamma,kappa;
30: TS ts;
31: Vec U;
32: } AppCtx;
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*);
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;
51: VecDuplicate(U,&Uob);
52: PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
53: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
54: VecLoad(Uob,viewer);
55: PetscViewerDestroy(&viewer);
56: VecAYPX(Uob,-1.,U);
57: VecScale(Uob,2.0);
58: VecAXPY(lambda,1.,Uob);
59: VecDestroy(&Uob);
60: return(0);
61: }
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: {
70: PetscViewerCreate(PetscObjectComm((PetscObject)da),viewer);
71: PetscViewerSetType(*viewer,PETSCVIEWERBINARY);
72: PetscViewerFileSetMode(*viewer,FILE_MODE_WRITE);
73: PetscViewerFileSetName(*viewer,filename);
74: return(0);
75: }
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;
87: TSGetDM(ts,&da);
88: TSSolve(ts,U);
89: PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
90: OutputBIN(da,filename,&viewer);
91: VecView(U,viewer);
92: PetscViewerDestroy(&viewer);
93: return(0);
94: }
96: PetscErrorCode InitialConditions(DM da,Vec U)
97: {
99: PetscInt i,j,xs,ys,xm,ym,Mx,My;
100: Field **u;
101: PetscReal hx,hy,x,y;
104: 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);
106: hx = 2.5/(PetscReal)Mx;
107: hy = 2.5/(PetscReal)My;
109: DMDAVecGetArray(da,U,&u);
110: /* Get local grid boundaries */
111: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
121: u[j][i].u = 1.0 - 2.0*u[j][i].v;
122: }
123: }
125: DMDAVecRestoreArray(da,U,&u);
126: return(0);
127: }
129: PetscErrorCode PerturbedInitialConditions(DM da,Vec U)
130: {
132: PetscInt i,j,xs,ys,xm,ym,Mx,My;
133: Field **u;
134: PetscReal hx,hy,x,y;
137: 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);
139: hx = 2.5/(PetscReal)Mx;
140: hy = 2.5/(PetscReal)My;
142: DMDAVecGetArray(da,U,&u);
143: /* Get local grid boundaries */
144: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
154: u[j][i].u = 1.0 - 2.0*u[j][i].v;
155: }
156: }
158: DMDAVecRestoreArray(da,U,&u);
159: return(0);
160: }
162: PetscErrorCode PerturbedInitialConditions2(DM da,Vec U)
163: {
165: PetscInt i,j,xs,ys,xm,ym,Mx,My;
166: Field **u;
167: PetscReal hx,hy,x,y;
170: 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);
172: hx = 2.5/(PetscReal)Mx;
173: hy = 2.5/(PetscReal)My;
175: DMDAVecGetArray(da,U,&u);
176: /* Get local grid boundaries */
177: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
187: u[j][i].u = 1.0 - 2.0*u[j][i].v;
188: }
189: }
191: DMDAVecRestoreArray(da,U,&u);
192: return(0);
193: }
195: PetscErrorCode PerturbedInitialConditions3(DM da,Vec U)
196: {
198: PetscInt i,j,xs,ys,xm,ym,Mx,My;
199: Field **u;
200: PetscReal hx,hy,x,y;
203: 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);
205: hx = 2.5/(PetscReal)Mx;
206: hy = 2.5/(PetscReal)My;
208: DMDAVecGetArray(da,U,&u);
209: /* Get local grid boundaries */
210: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
220: u[j][i].u = 1.0 - 2.0*u[j][i].v;
221: }
222: }
224: DMDAVecRestoreArray(da,U,&u);
225: return(0);
226: }
228: int main(int argc,char **argv)
229: {
231: DM da;
232: AppCtx appctx;
233: PetscBool forwardonly = PETSC_FALSE,implicitform = PETSC_FALSE;
234: PetscInt perturbic = 1;
236: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
237: PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);
238: PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);
239: PetscOptionsGetInt(NULL,NULL,"-perturbic",&perturbic,NULL);
241: appctx.D1 = 8.0e-5;
242: appctx.D2 = 4.0e-5;
243: appctx.gamma = .024;
244: appctx.kappa = .06;
246: /* Create distributed array (DMDA) to manage parallel grid and vectors */
247: DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,64,64,PETSC_DECIDE,PETSC_DECIDE,2,1,NULL,NULL,&da);
248: DMSetFromOptions(da);
249: DMSetUp(da);
250: DMDASetFieldName(da,0,"u");
251: DMDASetFieldName(da,1,"v");
253: /* Extract global vectors from DMDA; then duplicate for remaining
254: vectors that are the same types */
255: DMCreateGlobalVector(da,&appctx.U);
257: /* Create timestepping solver context */
258: TSCreate(PETSC_COMM_WORLD,&appctx.ts);
259: TSSetType(appctx.ts,TSCN);
260: TSSetDM(appctx.ts,da);
261: TSSetProblemType(appctx.ts,TS_NONLINEAR);
262: TSSetEquationType(appctx.ts,TS_EQ_ODE_EXPLICIT); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
263: if (!implicitform) {
264: TSSetRHSFunction(appctx.ts,NULL,RHSFunction,&appctx);
265: TSSetRHSJacobian(appctx.ts,NULL,NULL,RHSJacobian,&appctx);
266: } else {
267: TSSetIFunction(appctx.ts,NULL,IFunction,&appctx);
268: TSSetIJacobian(appctx.ts,NULL,NULL,IJacobian,&appctx);
269: }
271: /* Set initial conditions */
272: InitialConditions(da,appctx.U);
273: TSSetSolution(appctx.ts,appctx.U);
275: /* Set solver options */
276: TSSetMaxTime(appctx.ts,2000.0);
277: TSSetTimeStep(appctx.ts,0.5);
278: TSSetExactFinalTime(appctx.ts,TS_EXACTFINALTIME_MATCHSTEP);
279: TSSetFromOptions(appctx.ts);
281: GenerateOBs(appctx.ts,appctx.U,&appctx);
283: if (!forwardonly) {
284: Tao tao;
285: Vec P;
286: Vec lambda[1];
287: PetscLogStage opt_stage;
289: PetscLogStageRegister("Optimization",&opt_stage);
290: PetscLogStagePush(opt_stage);
291: if (perturbic == 1) {
292: PerturbedInitialConditions(da,appctx.U);
293: } else if (perturbic == 2) {
294: PerturbedInitialConditions2(da,appctx.U);
295: } else if (perturbic == 3) {
296: PerturbedInitialConditions3(da,appctx.U);
297: }
299: VecDuplicate(appctx.U,&lambda[0]);
300: TSSetCostGradients(appctx.ts,1,lambda,NULL);
302: /* Have the TS save its trajectory needed by TSAdjointSolve() */
303: TSSetSaveTrajectory(appctx.ts);
305: /* Create TAO solver and set desired solution method */
306: TaoCreate(PETSC_COMM_WORLD,&tao);
307: TaoSetType(tao,TAOBLMVM);
309: /* Set initial guess for TAO */
310: VecDuplicate(appctx.U,&P);
311: VecCopy(appctx.U,P);
312: TaoSetInitialVector(tao,P);
314: /* Set routine for function and gradient evaluation */
315: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);
317: /* Check for any TAO command line options */
318: TaoSetFromOptions(tao);
320: TaoSolve(tao);
321: TaoDestroy(&tao);
322: VecDestroy(&lambda[0]);
323: VecDestroy(&P);
324: PetscLogStagePop();
325: }
327: /* Free work space. All PETSc objects should be destroyed when they
328: are no longer needed. */
329: VecDestroy(&appctx.U);
330: TSDestroy(&appctx.ts);
331: DMDestroy(&da);
332: PetscFinalize();
333: return ierr;
334: }
336: /* ------------------------ TS callbacks ---------------------------- */
338: /*
339: RHSFunction - Evaluates nonlinear function, F(x).
341: Input Parameters:
342: . ts - the TS context
343: . X - input vector
344: . ptr - optional user-defined context, as set by TSSetRHSFunction()
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;
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;
361: TSGetDM(ts,&da);
362: DMGetLocalVector(da,&localU);
363: 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);
364: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
365: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
372: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
374: /* Get pointers to vector data */
375: DMDAVecGetArrayRead(da,localU,&u);
376: DMDAVecGetArray(da,F,&f);
378: /* Get local grid boundaries */
379: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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: PetscLogFlops(16*xm*ym);
396: DMDAVecRestoreArrayRead(da,localU,&u);
397: DMDAVecRestoreArray(da,F,&f);
398: DMRestoreLocalVector(da,&localU);
399: return(0);
400: }
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;
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];
416: TSGetDM(ts,&da);
417: DMGetLocalVector(da,&localU);
418: 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);
420: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
421: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
428: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
430: /* Get pointers to vector data */
431: DMDAVecGetArrayRead(da,localU,&u);
433: /* Get local grid boundaries */
434: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
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;
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);*/
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;
471: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;
480: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
481: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
482: }
483: }
484: PetscLogFlops(19*xm*ym);
486: DMDAVecRestoreArrayRead(da,localU,&u);
487: DMRestoreLocalVector(da,&localU);
488: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
489: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
490: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
491: return(0);
492: }
494: /*
495: IFunction - Evaluates implicit nonlinear function, xdot - F(x).
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()
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;
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;
518: TSGetDM(ts,&da);
519: DMGetLocalVector(da,&localU);
520: 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);
521: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
522: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
529: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
531: /* Get pointers to vector data */
532: DMDAVecGetArrayRead(da,localU,&u);
533: DMDAVecGetArray(da,F,&f);
534: DMDAVecGetArrayRead(da,Udot,&udot);
536: /* Get local grid boundaries */
537: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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: PetscLogFlops(16*xm*ym);
554: DMDAVecRestoreArrayRead(da,localU,&u);
555: DMDAVecRestoreArray(da,F,&f);
556: DMDAVecRestoreArrayRead(da,Udot,&udot);
557: DMRestoreLocalVector(da,&localU);
558: return(0);
559: }
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;
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];
575: TSGetDM(ts,&da);
576: DMGetLocalVector(da,&localU);
577: 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);
579: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
580: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
587: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
589: /* Get pointers to vector data */
590: DMDAVecGetArrayRead(da,localU,&u);
592: /* Get local grid boundaries */
593: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
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;
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;
629: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;
638: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
639: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
640: }
641: }
642: PetscLogFlops(19*xm*ym);
644: DMDAVecRestoreArrayRead(da,localU,&u);
645: DMRestoreLocalVector(da,&localU);
646: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
647: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
648: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
649: return(0);
650: }
652: /* ------------------------ TAO callbacks ---------------------------- */
654: /*
655: FormFunctionAndGradient - Evaluates the function and corresponding gradient.
657: Input Parameters:
658: tao - the Tao context
659: P - the input vector
660: ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
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;
678: TSSetTime(appctx->ts,0.0);
679: TSGetTimeStep(appctx->ts,×tep);
680: if (timestep<0) {
681: TSSetTimeStep(appctx->ts,-timestep);
682: }
683: TSSetStepNumber(appctx->ts,0);
684: TSSetFromOptions(appctx->ts);
686: VecDuplicate(P,&SDiff);
687: VecCopy(P,appctx->U);
688: TSGetDM(appctx->ts,&da);
689: *f = 0;
691: TSSolve(appctx->ts,appctx->U);
692: PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
693: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
694: VecLoad(SDiff,viewer);
695: PetscViewerDestroy(&viewer);
696: VecAYPX(SDiff,-1.,appctx->U);
697: VecDot(SDiff,SDiff,&soberr);
698: *f += soberr;
700: TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);
701: VecSet(lambda[0],0.0);
702: InitializeLambda(da,lambda[0],appctx->U,appctx);
703: TSAdjointSolve(appctx->ts);
705: VecCopy(lambda[0],G);
707: VecDestroy(&SDiff);
708: return(0);
709: }
711: /*TEST
713: build:
714: requires: !complex !single
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
720: TEST*/