Actual source code: ex5opt_ic.c
petsc-3.14.5 2021-03-03
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: #if defined(PETSC_USE_LOG)
288: PetscLogStage opt_stage;
289: #endif
291: PetscLogStageRegister("Optimization",&opt_stage);
292: PetscLogStagePush(opt_stage);
293: if (perturbic == 1) {
294: PerturbedInitialConditions(da,appctx.U);
295: } else if (perturbic == 2) {
296: PerturbedInitialConditions2(da,appctx.U);
297: } else if (perturbic == 3) {
298: PerturbedInitialConditions3(da,appctx.U);
299: }
301: VecDuplicate(appctx.U,&lambda[0]);
302: TSSetCostGradients(appctx.ts,1,lambda,NULL);
304: /* Have the TS save its trajectory needed by TSAdjointSolve() */
305: TSSetSaveTrajectory(appctx.ts);
307: /* Create TAO solver and set desired solution method */
308: TaoCreate(PETSC_COMM_WORLD,&tao);
309: TaoSetType(tao,TAOBLMVM);
311: /* Set initial guess for TAO */
312: VecDuplicate(appctx.U,&P);
313: VecCopy(appctx.U,P);
314: TaoSetInitialVector(tao,P);
316: /* Set routine for function and gradient evaluation */
317: TaoSetObjectiveAndGradientRoutine(tao,FormFunctionAndGradient,&appctx);
319: /* Check for any TAO command line options */
320: TaoSetFromOptions(tao);
322: TaoSolve(tao);
323: TaoDestroy(&tao);
324: VecDestroy(&lambda[0]);
325: VecDestroy(&P);
326: PetscLogStagePop();
327: }
329: /* Free work space. All PETSc objects should be destroyed when they
330: are no longer needed. */
331: VecDestroy(&appctx.U);
332: TSDestroy(&appctx.ts);
333: DMDestroy(&da);
334: PetscFinalize();
335: return ierr;
336: }
338: /* ------------------------ TS callbacks ---------------------------- */
340: /*
341: RHSFunction - Evaluates nonlinear function, F(x).
343: Input Parameters:
344: . ts - the TS context
345: . X - input vector
346: . ptr - optional user-defined context, as set by TSSetRHSFunction()
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;
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;
363: TSGetDM(ts,&da);
364: DMGetLocalVector(da,&localU);
365: 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);
366: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
367: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
374: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
376: /* Get pointers to vector data */
377: DMDAVecGetArrayRead(da,localU,&u);
378: DMDAVecGetArray(da,F,&f);
380: /* Get local grid boundaries */
381: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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: PetscLogFlops(16.0*xm*ym);
398: DMDAVecRestoreArrayRead(da,localU,&u);
399: DMDAVecRestoreArray(da,F,&f);
400: DMRestoreLocalVector(da,&localU);
401: return(0);
402: }
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;
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];
418: TSGetDM(ts,&da);
419: DMGetLocalVector(da,&localU);
420: 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);
422: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
423: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
430: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
432: /* Get pointers to vector data */
433: DMDAVecGetArrayRead(da,localU,&u);
435: /* Get local grid boundaries */
436: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
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;
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);*/
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;
473: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;
482: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
483: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
484: }
485: }
486: PetscLogFlops(19.0*xm*ym);
488: DMDAVecRestoreArrayRead(da,localU,&u);
489: DMRestoreLocalVector(da,&localU);
490: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
491: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
492: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
493: return(0);
494: }
496: /*
497: IFunction - Evaluates implicit nonlinear function, xdot - F(x).
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()
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;
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;
520: TSGetDM(ts,&da);
521: DMGetLocalVector(da,&localU);
522: 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);
523: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
524: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
531: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
533: /* Get pointers to vector data */
534: DMDAVecGetArrayRead(da,localU,&u);
535: DMDAVecGetArray(da,F,&f);
536: DMDAVecGetArrayRead(da,Udot,&udot);
538: /* Get local grid boundaries */
539: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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: PetscLogFlops(16.0*xm*ym);
556: DMDAVecRestoreArrayRead(da,localU,&u);
557: DMDAVecRestoreArray(da,F,&f);
558: DMDAVecRestoreArrayRead(da,Udot,&udot);
559: DMRestoreLocalVector(da,&localU);
560: return(0);
561: }
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;
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];
577: TSGetDM(ts,&da);
578: DMGetLocalVector(da,&localU);
579: 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);
581: hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
582: hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
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: DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);
589: DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);
591: /* Get pointers to vector data */
592: DMDAVecGetArrayRead(da,localU,&u);
594: /* Get local grid boundaries */
595: DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);
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;
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;
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;
631: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
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;
640: MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);
641: /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
642: }
643: }
644: PetscLogFlops(19.0*xm*ym);
646: DMDAVecRestoreArrayRead(da,localU,&u);
647: DMRestoreLocalVector(da,&localU);
648: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
649: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
650: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
651: return(0);
652: }
654: /* ------------------------ TAO callbacks ---------------------------- */
656: /*
657: FormFunctionAndGradient - Evaluates the function and corresponding gradient.
659: Input Parameters:
660: tao - the Tao context
661: P - the input vector
662: ctx - optional user-defined context, as set by TaoSetObjectiveAndGradientRoutine()
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;
680: TSSetTime(appctx->ts,0.0);
681: TSGetTimeStep(appctx->ts,×tep);
682: if (timestep<0) {
683: TSSetTimeStep(appctx->ts,-timestep);
684: }
685: TSSetStepNumber(appctx->ts,0);
686: TSSetFromOptions(appctx->ts);
688: VecDuplicate(P,&SDiff);
689: VecCopy(P,appctx->U);
690: TSGetDM(appctx->ts,&da);
691: *f = 0;
693: TSSolve(appctx->ts,appctx->U);
694: PetscSNPrintf(filename,sizeof filename,"ex5opt.ob");
695: PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
696: VecLoad(SDiff,viewer);
697: PetscViewerDestroy(&viewer);
698: VecAYPX(SDiff,-1.,appctx->U);
699: VecDot(SDiff,SDiff,&soberr);
700: *f += soberr;
702: TSGetCostGradients(appctx->ts,NULL,&lambda,NULL);
703: VecSet(lambda[0],0.0);
704: InitializeLambda(da,lambda[0],appctx->U,appctx);
705: TSAdjointSolve(appctx->ts);
707: VecCopy(lambda[0],G);
709: VecDestroy(&SDiff);
710: return(0);
711: }
713: /*TEST
715: build:
716: requires: !complex !single
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
722: TEST*/