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 interface to a system of time-dependent partial differential equations.
6 It computes the sensitivity of a component in the final solution, which locates in the center of the 2D domain, w.r.t. the initial conditions.
7 The user does not need to provide any additional functions. The required functions in the original simulation are reused in the adjoint run.
8
9 Runtime options:
10 -forwardonly - run the forward simulation without adjoint
11 -implicitform - provide IFunction and IJacobian to TS, if not set, RHSFunction and RHSJacobian will be used
12 -aijpc - set the preconditioner matrix to be aij (the Jacobian matrix can be of a different type such as ELL)
13 */
14
15 #include <petscsys.h>
16 #include <petscdm.h>
17 #include <petscdmda.h>
18 #include <petscts.h>
19
20 typedef struct {
21 PetscScalar u,v;
22 } Field;
23
24 typedef struct {
25 PetscReal D1,D2,gamma,kappa;
26 PetscBool aijpc;
27 } AppCtx;
28
29 /*
30 User-defined routines
31 */
32 PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
33 PetscErrorCode InitialConditions(DM,Vec);
34 PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
35 PetscErrorCode IFunction(TS,PetscReal,Vec,Vec,Vec,void*);
36 PetscErrorCode IJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);
37
InitializeLambda(DM da,Vec lambda,PetscReal x,PetscReal y)38 PetscErrorCode InitializeLambda(DM da,Vec lambda,PetscReal x,PetscReal y)
39 {
40 PetscInt i,j,Mx,My,xs,ys,xm,ym;
41 PetscErrorCode ierr;
42 Field **l;
43 PetscFunctionBegin;
44
45 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);
46 /* locate the global i index for x and j index for y */
47 i = (PetscInt)(x*(Mx-1));
48 j = (PetscInt)(y*(My-1));
49 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
50
51 if (xs <= i && i < xs+xm && ys <= j && j < ys+ym) {
52 /* the i,j vertex is on this process */
53 ierr = DMDAVecGetArray(da,lambda,&l);CHKERRQ(ierr);
54 l[j][i].u = 1.0;
55 l[j][i].v = 1.0;
56 ierr = DMDAVecRestoreArray(da,lambda,&l);CHKERRQ(ierr);
57 }
58 PetscFunctionReturn(0);
59 }
60
main(int argc,char ** argv)61 int main(int argc,char **argv)
62 {
63 TS ts; /* ODE integrator */
64 Vec x; /* solution */
65 PetscErrorCode ierr;
66 DM da;
67 AppCtx appctx;
68 Vec lambda[1];
69 PetscBool forwardonly=PETSC_FALSE,implicitform=PETSC_TRUE;
70
71 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
72 ierr = PetscOptionsGetBool(NULL,NULL,"-forwardonly",&forwardonly,NULL);CHKERRQ(ierr);
73 ierr = PetscOptionsGetBool(NULL,NULL,"-implicitform",&implicitform,NULL);CHKERRQ(ierr);
74 appctx.aijpc = PETSC_FALSE;
75 ierr = PetscOptionsGetBool(NULL,NULL,"-aijpc",&appctx.aijpc,NULL);CHKERRQ(ierr);
76
77 appctx.D1 = 8.0e-5;
78 appctx.D2 = 4.0e-5;
79 appctx.gamma = .024;
80 appctx.kappa = .06;
81
82 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 Create distributed array (DMDA) to manage parallel grid and vectors
84 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 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);
86 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
87 ierr = DMSetUp(da);CHKERRQ(ierr);
88 ierr = DMDASetFieldName(da,0,"u");CHKERRQ(ierr);
89 ierr = DMDASetFieldName(da,1,"v");CHKERRQ(ierr);
90
91 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92 Extract global vectors from DMDA; then duplicate for remaining
93 vectors that are the same types
94 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr);
96
97 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
98 Create timestepping solver context
99 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
100 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
101 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
102 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
103 ierr = TSSetEquationType(ts,TS_EQ_ODE_EXPLICIT);CHKERRQ(ierr); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
104 if (!implicitform) {
105 ierr = TSSetType(ts,TSRK);CHKERRQ(ierr);
106 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&appctx);CHKERRQ(ierr);
107 ierr = TSSetRHSJacobian(ts,NULL,NULL,RHSJacobian,&appctx);CHKERRQ(ierr);
108 } else {
109 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr);
110 ierr = TSSetIFunction(ts,NULL,IFunction,&appctx);CHKERRQ(ierr);
111 if (appctx.aijpc) {
112 Mat A,B;
113
114 ierr = DMSetMatType(da,MATSELL);CHKERRQ(ierr);
115 ierr = DMCreateMatrix(da,&A);CHKERRQ(ierr);
116 ierr = MatConvert(A,MATAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr);
117 /* FIXME do we need to change viewer to display matrix in natural ordering as DMCreateMatrix_DA does? */
118 ierr = TSSetIJacobian(ts,A,B,IJacobian,&appctx);CHKERRQ(ierr);
119 ierr = MatDestroy(&A);CHKERRQ(ierr);
120 ierr = MatDestroy(&B);CHKERRQ(ierr);
121 } else {
122 ierr = TSSetIJacobian(ts,NULL,NULL,IJacobian,&appctx);CHKERRQ(ierr);
123 }
124 }
125
126 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 Set initial conditions
128 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129 ierr = InitialConditions(da,x);CHKERRQ(ierr);
130 ierr = TSSetSolution(ts,x);CHKERRQ(ierr);
131
132 /*
133 Have the TS save its trajectory so that TSAdjointSolve() may be used
134 */
135 if (!forwardonly) { ierr = TSSetSaveTrajectory(ts);CHKERRQ(ierr); }
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Set solver options
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140 ierr = TSSetMaxTime(ts,200.0);CHKERRQ(ierr);
141 ierr = TSSetTimeStep(ts,0.5);CHKERRQ(ierr);
142 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_MATCHSTEP);CHKERRQ(ierr);
143 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
144
145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 Solve ODE system
147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148 ierr = TSSolve(ts,x);CHKERRQ(ierr);
149 if (!forwardonly) {
150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 Start the Adjoint model
152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153 ierr = VecDuplicate(x,&lambda[0]);CHKERRQ(ierr);
154 /* Reset initial conditions for the adjoint integration */
155 ierr = InitializeLambda(da,lambda[0],0.5,0.5);CHKERRQ(ierr);
156 ierr = TSSetCostGradients(ts,1,lambda,NULL);CHKERRQ(ierr);
157 ierr = TSAdjointSolve(ts);CHKERRQ(ierr);
158 ierr = VecDestroy(&lambda[0]);CHKERRQ(ierr);
159 }
160 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
161 Free work space. All PETSc objects should be destroyed when they
162 are no longer needed.
163 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
164 ierr = VecDestroy(&x);CHKERRQ(ierr);
165 ierr = TSDestroy(&ts);CHKERRQ(ierr);
166 ierr = DMDestroy(&da);CHKERRQ(ierr);
167 ierr = PetscFinalize();
168 return ierr;
169 }
170
171 /* ------------------------------------------------------------------- */
172 /*
173 RHSFunction - Evaluates nonlinear function, F(x).
174
175 Input Parameters:
176 . ts - the TS context
177 . X - input vector
178 . ptr - optional user-defined context, as set by TSSetRHSFunction()
179
180 Output Parameter:
181 . F - function vector
182 */
RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void * ptr)183 PetscErrorCode RHSFunction(TS ts,PetscReal ftime,Vec U,Vec F,void *ptr)
184 {
185 AppCtx *appctx = (AppCtx*)ptr;
186 DM da;
187 PetscErrorCode ierr;
188 PetscInt i,j,Mx,My,xs,ys,xm,ym;
189 PetscReal hx,hy,sx,sy;
190 PetscScalar uc,uxx,uyy,vc,vxx,vyy;
191 Field **u,**f;
192 Vec localU;
193
194 PetscFunctionBegin;
195 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
196 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
197 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);
198 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
199 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
200
201 /*
202 Scatter ghost points to local vector,using the 2-step process
203 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
204 By placing code between these two statements, computations can be
205 done while messages are in transition.
206 */
207 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
208 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
209
210 /*
211 Get pointers to vector data
212 */
213 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
214 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
215
216 /*
217 Get local grid boundaries
218 */
219 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
220
221 /*
222 Compute function over the locally owned part of the grid
223 */
224 for (j=ys; j<ys+ym; j++) {
225 for (i=xs; i<xs+xm; i++) {
226 uc = u[j][i].u;
227 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
228 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
229 vc = u[j][i].v;
230 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
231 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
232 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);
233 f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc;
234 }
235 }
236 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
237
238 /*
239 Restore vectors
240 */
241 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
242 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
243 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
244 PetscFunctionReturn(0);
245 }
246
247 /* ------------------------------------------------------------------- */
InitialConditions(DM da,Vec U)248 PetscErrorCode InitialConditions(DM da,Vec U)
249 {
250 PetscErrorCode ierr;
251 PetscInt i,j,xs,ys,xm,ym,Mx,My;
252 Field **u;
253 PetscReal hx,hy,x,y;
254
255 PetscFunctionBegin;
256 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);
257
258 hx = 2.5/(PetscReal)Mx;
259 hy = 2.5/(PetscReal)My;
260
261 /*
262 Get pointers to vector data
263 */
264 ierr = DMDAVecGetArray(da,U,&u);CHKERRQ(ierr);
265
266 /*
267 Get local grid boundaries
268 */
269 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
270
271 /*
272 Compute function over the locally owned part of the grid
273 */
274 for (j=ys; j<ys+ym; j++) {
275 y = j*hy;
276 for (i=xs; i<xs+xm; i++) {
277 x = i*hx;
278 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);
279 else u[j][i].v = 0.0;
280
281 u[j][i].u = 1.0 - 2.0*u[j][i].v;
282 }
283 }
284
285 /*
286 Restore vectors
287 */
288 ierr = DMDAVecRestoreArray(da,U,&u);CHKERRQ(ierr);
289 PetscFunctionReturn(0);
290 }
291
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void * ctx)292 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat BB,void *ctx)
293 {
294 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
295 DM da;
296 PetscErrorCode ierr;
297 PetscInt i,j,Mx,My,xs,ys,xm,ym;
298 PetscReal hx,hy,sx,sy;
299 PetscScalar uc,vc;
300 Field **u;
301 Vec localU;
302 MatStencil stencil[6],rowstencil;
303 PetscScalar entries[6];
304
305 PetscFunctionBegin;
306 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
307 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
308 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);
309
310 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
311 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
312
313 /*
314 Scatter ghost points to local vector,using the 2-step process
315 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
316 By placing code between these two statements, computations can be
317 done while messages are in transition.
318 */
319 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
320 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
321
322 /*
323 Get pointers to vector data
324 */
325 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
326
327 /*
328 Get local grid boundaries
329 */
330 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
331
332 stencil[0].k = 0;
333 stencil[1].k = 0;
334 stencil[2].k = 0;
335 stencil[3].k = 0;
336 stencil[4].k = 0;
337 stencil[5].k = 0;
338 rowstencil.k = 0;
339 /*
340 Compute function over the locally owned part of the grid
341 */
342 for (j=ys; j<ys+ym; j++) {
343
344 stencil[0].j = j-1;
345 stencil[1].j = j+1;
346 stencil[2].j = j;
347 stencil[3].j = j;
348 stencil[4].j = j;
349 stencil[5].j = j;
350 rowstencil.k = 0; rowstencil.j = j;
351 for (i=xs; i<xs+xm; i++) {
352 uc = u[j][i].u;
353 vc = u[j][i].v;
354
355 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
356 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
357
358 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
359 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
360 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/
361
362 stencil[0].i = i; stencil[0].c = 0; entries[0] = appctx->D1*sy;
363 stencil[1].i = i; stencil[1].c = 0; entries[1] = appctx->D1*sy;
364 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = appctx->D1*sx;
365 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = appctx->D1*sx;
366 stencil[4].i = i; stencil[4].c = 0; entries[4] = -2.0*appctx->D1*(sx + sy) - vc*vc - appctx->gamma;
367 stencil[5].i = i; stencil[5].c = 1; entries[5] = -2.0*uc*vc;
368 rowstencil.i = i; rowstencil.c = 0;
369
370 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
371 if (appctx->aijpc) {
372 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
373 }
374 stencil[0].c = 1; entries[0] = appctx->D2*sy;
375 stencil[1].c = 1; entries[1] = appctx->D2*sy;
376 stencil[2].c = 1; entries[2] = appctx->D2*sx;
377 stencil[3].c = 1; entries[3] = appctx->D2*sx;
378 stencil[4].c = 1; entries[4] = -2.0*appctx->D2*(sx + sy) + 2.0*uc*vc - appctx->gamma - appctx->kappa;
379 stencil[5].c = 0; entries[5] = vc*vc;
380 rowstencil.c = 1;
381
382 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
383 if (appctx->aijpc) {
384 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
385 }
386 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
387 }
388 }
389
390 /*
391 Restore vectors
392 */
393 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr);
394 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
395 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
396 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
397 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
398 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
399 if (appctx->aijpc) {
400 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
401 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
402 ierr = MatSetOption(BB,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
403 }
404
405 PetscFunctionReturn(0);
406 }
407
408 /*
409 IFunction - Evaluates implicit nonlinear function, xdot - F(x).
410
411 Input Parameters:
412 . ts - the TS context
413 . U - input vector
414 . Udot - input vector
415 . ptr - optional user-defined context, as set by TSSetRHSFunction()
416
417 Output Parameter:
418 . F - function vector
419 */
IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void * ptr)420 PetscErrorCode IFunction(TS ts,PetscReal ftime,Vec U,Vec Udot,Vec F,void *ptr)
421 {
422 AppCtx *appctx = (AppCtx*)ptr;
423 DM da;
424 PetscErrorCode ierr;
425 PetscInt i,j,Mx,My,xs,ys,xm,ym;
426 PetscReal hx,hy,sx,sy;
427 PetscScalar uc,uxx,uyy,vc,vxx,vyy;
428 Field **u,**f,**udot;
429 Vec localU;
430
431 PetscFunctionBegin;
432 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
433 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
434 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);
435 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
436 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
437
438 /*
439 Scatter ghost points to local vector,using the 2-step process
440 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
441 By placing code between these two statements, computations can be
442 done while messages are in transition.
443 */
444 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
445 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
446
447 /*
448 Get pointers to vector data
449 */
450 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
451 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr);
452 ierr = DMDAVecGetArrayRead(da,Udot,&udot);CHKERRQ(ierr);
453
454 /*
455 Get local grid boundaries
456 */
457 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
458
459 /*
460 Compute function over the locally owned part of the grid
461 */
462 for (j=ys; j<ys+ym; j++) {
463 for (i=xs; i<xs+xm; i++) {
464 uc = u[j][i].u;
465 uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
466 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
467 vc = u[j][i].v;
468 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
469 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
470 f[j][i].u = udot[j][i].u - ( appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc));
471 f[j][i].v = udot[j][i].v - ( appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc);
472 }
473 }
474 ierr = PetscLogFlops(16.0*xm*ym);CHKERRQ(ierr);
475
476 /*
477 Restore vectors
478 */
479 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
480 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr);
481 ierr = DMDAVecRestoreArrayRead(da,Udot,&udot);CHKERRQ(ierr);
482 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
483 PetscFunctionReturn(0);
484 }
485
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void * ctx)486 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat BB,void *ctx)
487 {
488 AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
489 DM da;
490 PetscErrorCode ierr;
491 PetscInt i,j,Mx,My,xs,ys,xm,ym;
492 PetscReal hx,hy,sx,sy;
493 PetscScalar uc,vc;
494 Field **u;
495 Vec localU;
496 MatStencil stencil[6],rowstencil;
497 PetscScalar entries[6];
498
499 PetscFunctionBegin;
500 ierr = TSGetDM(ts,&da);CHKERRQ(ierr);
501 ierr = DMGetLocalVector(da,&localU);CHKERRQ(ierr);
502 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);
503
504 hx = 2.50/(PetscReal)Mx; sx = 1.0/(hx*hx);
505 hy = 2.50/(PetscReal)My; sy = 1.0/(hy*hy);
506
507 /*
508 Scatter ghost points to local vector,using the 2-step process
509 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
510 By placing code between these two statements, computations can be
511 done while messages are in transition.
512 */
513 ierr = DMGlobalToLocalBegin(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
514 ierr = DMGlobalToLocalEnd(da,U,INSERT_VALUES,localU);CHKERRQ(ierr);
515
516 /*
517 Get pointers to vector data
518 */
519 ierr = DMDAVecGetArrayRead(da,localU,&u);CHKERRQ(ierr);
520
521 /*
522 Get local grid boundaries
523 */
524 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
525
526 stencil[0].k = 0;
527 stencil[1].k = 0;
528 stencil[2].k = 0;
529 stencil[3].k = 0;
530 stencil[4].k = 0;
531 stencil[5].k = 0;
532 rowstencil.k = 0;
533 /*
534 Compute function over the locally owned part of the grid
535 */
536 for (j=ys; j<ys+ym; j++) {
537
538 stencil[0].j = j-1;
539 stencil[1].j = j+1;
540 stencil[2].j = j;
541 stencil[3].j = j;
542 stencil[4].j = j;
543 stencil[5].j = j;
544 rowstencil.k = 0; rowstencil.j = j;
545 for (i=xs; i<xs+xm; i++) {
546 uc = u[j][i].u;
547 vc = u[j][i].v;
548
549 /* uxx = (-2.0*uc + u[j][i-1].u + u[j][i+1].u)*sx;
550 uyy = (-2.0*uc + u[j-1][i].u + u[j+1][i].u)*sy;
551
552 vxx = (-2.0*vc + u[j][i-1].v + u[j][i+1].v)*sx;
553 vyy = (-2.0*vc + u[j-1][i].v + u[j+1][i].v)*sy;
554 f[j][i].u = appctx->D1*(uxx + uyy) - uc*vc*vc + appctx->gamma*(1.0 - uc);*/
555
556 stencil[0].i = i; stencil[0].c = 0; entries[0] = -appctx->D1*sy;
557 stencil[1].i = i; stencil[1].c = 0; entries[1] = -appctx->D1*sy;
558 stencil[2].i = i-1; stencil[2].c = 0; entries[2] = -appctx->D1*sx;
559 stencil[3].i = i+1; stencil[3].c = 0; entries[3] = -appctx->D1*sx;
560 stencil[4].i = i; stencil[4].c = 0; entries[4] = 2.0*appctx->D1*(sx + sy) + vc*vc + appctx->gamma + a;
561 stencil[5].i = i; stencil[5].c = 1; entries[5] = 2.0*uc*vc;
562 rowstencil.i = i; rowstencil.c = 0;
563
564 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
565 if (appctx->aijpc) {
566 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
567 }
568 stencil[0].c = 1; entries[0] = -appctx->D2*sy;
569 stencil[1].c = 1; entries[1] = -appctx->D2*sy;
570 stencil[2].c = 1; entries[2] = -appctx->D2*sx;
571 stencil[3].c = 1; entries[3] = -appctx->D2*sx;
572 stencil[4].c = 1; entries[4] = 2.0*appctx->D2*(sx + sy) - 2.0*uc*vc + appctx->gamma + appctx->kappa + a;
573 stencil[5].c = 0; entries[5] = -vc*vc;
574 rowstencil.c = 1;
575
576 ierr = MatSetValuesStencil(A,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
577 if (appctx->aijpc) {
578 ierr = MatSetValuesStencil(BB,1,&rowstencil,6,stencil,entries,INSERT_VALUES);CHKERRQ(ierr);
579 }
580 /* f[j][i].v = appctx->D2*(vxx + vyy) + uc*vc*vc - (appctx->gamma + appctx->kappa)*vc; */
581 }
582 }
583
584 /*
585 Restore vectors
586 */
587 ierr = PetscLogFlops(19.0*xm*ym);CHKERRQ(ierr);
588 ierr = DMDAVecRestoreArrayRead(da,localU,&u);CHKERRQ(ierr);
589 ierr = DMRestoreLocalVector(da,&localU);CHKERRQ(ierr);
590 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
591 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
592 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
593 if (appctx->aijpc) {
594 ierr = MatAssemblyBegin(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
595 ierr = MatAssemblyEnd(BB,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
596 ierr = MatSetOption(BB,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);
597 }
598 PetscFunctionReturn(0);
599 }
600
601
602 /*TEST
603
604 build:
605 requires: !complex !single
606
607 test:
608 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -da_grid_x 20 -da_grid_y 20
609 output_file: output/ex5adj_1.out
610
611 test:
612 suffix: 2
613 nsize: 2
614 args: -ts_max_steps 10 -ts_dt 10 -ts_monitor -ts_adjoint_monitor -ksp_monitor_short -da_grid_x 20 -da_grid_y 20 -ts_trajectory_dirname Test-dir -ts_trajectory_file_template test-%06D.cp
615
616 test:
617 suffix: 3
618 nsize: 2
619 args: -ts_max_steps 10 -ts_dt 10 -ts_adjoint_monitor_draw_sensi
620
621 test:
622 suffix: 4
623 nsize: 2
624 args: -ts_max_steps 10 -ts_dt 10 -ts_monitor -ts_adjoint_monitor -ksp_monitor_short -da_grid_x 20 -da_grid_y 20 -snes_fd_color
625 output_file: output/ex5adj_2.out
626
627 test:
628 suffix: 5
629 nsize: 2
630 args: -ts_max_steps 10 -implicitform 0 -ts_type rk -ts_rk_type 4 -ts_monitor -ts_adjoint_monitor -da_grid_x 20 -da_grid_y 20 -snes_fd_color
631 output_file: output/ex5adj_1.out
632
633 test:
634 suffix: knl
635 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ts_trajectory_solution_only 0 -malloc_hbw -ts_trajectory_use_dram 1
636 output_file: output/ex5adj_3.out
637 requires: knl
638
639 test:
640 suffix: sell
641 nsize: 4
642 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type none
643 output_file: output/ex5adj_sell_1.out
644
645 test:
646 suffix: aijsell
647 nsize: 4
648 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell -pc_type none
649 output_file: output/ex5adj_sell_1.out
650
651 test:
652 suffix: sell2
653 nsize: 4
654 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type sor
655 output_file: output/ex5adj_sell_2.out
656
657 test:
658 suffix: aijsell2
659 nsize: 4
660 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type sor
661 output_file: output/ex5adj_sell_2.out
662
663 test:
664 suffix: sell3
665 nsize: 4
666 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type bjacobi -mg_levels_pc_type bjacobi
667 output_file: output/ex5adj_sell_3.out
668
669 test:
670 suffix: sell4
671 nsize: 4
672 args: -forwardonly -implicitform -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -pc_type mg -pc_mg_levels 2 -mg_coarse_pc_type bjacobi -mg_levels_pc_type bjacobi
673 output_file: output/ex5adj_sell_4.out
674
675 test:
676 suffix: sell5
677 nsize: 4
678 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type sell -aijpc
679 output_file: output/ex5adj_sell_5.out
680
681 test:
682 suffix: aijsell5
683 nsize: 4
684 args: -forwardonly -ts_max_steps 10 -ts_monitor -snes_monitor_short -dm_mat_type aijsell
685 output_file: output/ex5adj_sell_5.out
686
687 test:
688 suffix: sell6
689 args: -ts_max_steps 10 -ts_monitor -ts_adjoint_monitor -ts_trajectory_type memory -ts_trajectory_solution_only 0 -dm_mat_type sell -pc_type jacobi
690 output_file: output/ex5adj_sell_6.out
691
692 TEST*/
693