1 /*
2 The Problem:
3 Solve the convection-diffusion equation:
4
5 u_t+a*(u_x+u_y)=epsilon*(u_xx+u_yy)
6 u=0 at x=0, y=0
7 u_x=0 at x=1
8 u_y=0 at y=1
9 u = exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0))) at t=0
10
11 This program tests the routine of computing the Jacobian by the
12 finite difference method as well as PETSc.
13
14 */
15
16 static char help[] = "Solve the convection-diffusion equation. \n\n";
17
18 #include <petscts.h>
19
20 typedef struct
21 {
22 PetscInt m; /* the number of mesh points in x-direction */
23 PetscInt n; /* the number of mesh points in y-direction */
24 PetscReal dx; /* the grid space in x-direction */
25 PetscReal dy; /* the grid space in y-direction */
26 PetscReal a; /* the convection coefficient */
27 PetscReal epsilon; /* the diffusion coefficient */
28 PetscReal tfinal;
29 } Data;
30
31 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
32 extern PetscErrorCode Initial(Vec,void*);
33 extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
34 extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
35 extern PetscErrorCode PostStep(TS);
36
main(int argc,char ** argv)37 int main(int argc,char **argv)
38 {
39 PetscErrorCode ierr;
40 PetscInt time_steps=100,iout,NOUT=1;
41 Vec global;
42 PetscReal dt,ftime,ftime_original;
43 TS ts;
44 PetscViewer viewfile;
45 Mat J = 0;
46 Vec x;
47 Data data;
48 PetscInt mn;
49 PetscBool flg;
50 MatColoring mc;
51 ISColoring iscoloring;
52 MatFDColoring matfdcoloring = 0;
53 PetscBool fd_jacobian_coloring = PETSC_FALSE;
54 SNES snes;
55 KSP ksp;
56 PC pc;
57
58 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
59
60 /* set data */
61 data.m = 9;
62 data.n = 9;
63 data.a = 1.0;
64 data.epsilon = 0.1;
65 data.dx = 1.0/(data.m+1.0);
66 data.dy = 1.0/(data.n+1.0);
67 mn = (data.m)*(data.n);
68 ierr = PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);CHKERRQ(ierr);
69
70 /* set initial conditions */
71 ierr = VecCreate(PETSC_COMM_WORLD,&global);CHKERRQ(ierr);
72 ierr = VecSetSizes(global,PETSC_DECIDE,mn);CHKERRQ(ierr);
73 ierr = VecSetFromOptions(global);CHKERRQ(ierr);
74 ierr = Initial(global,&data);CHKERRQ(ierr);
75 ierr = VecDuplicate(global,&x);CHKERRQ(ierr);
76
77 /* create timestep context */
78 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
79 ierr = TSMonitorSet(ts,Monitor,&data,NULL);CHKERRQ(ierr);
80 ierr = TSSetType(ts,TSEULER);CHKERRQ(ierr);
81 dt = 0.1;
82 ftime_original = data.tfinal = 1.0;
83
84 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
85 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr);
86 ierr = TSSetMaxTime(ts,ftime_original);CHKERRQ(ierr);
87 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
88 ierr = TSSetSolution(ts,global);CHKERRQ(ierr);
89
90 /* set user provided RHSFunction and RHSJacobian */
91 ierr = TSSetRHSFunction(ts,NULL,RHSFunction,&data);CHKERRQ(ierr);
92 ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);
93 ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);CHKERRQ(ierr);
94 ierr = MatSetFromOptions(J);CHKERRQ(ierr);
95 ierr = MatSeqAIJSetPreallocation(J,5,NULL);CHKERRQ(ierr);
96 ierr = MatMPIAIJSetPreallocation(J,5,NULL,5,NULL);CHKERRQ(ierr);
97
98 ierr = PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);CHKERRQ(ierr);
99 if (!flg) {
100 ierr = TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);CHKERRQ(ierr);
101 } else {
102 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
103 ierr = PetscOptionsHasName(NULL,NULL,"-fd_color",&fd_jacobian_coloring);CHKERRQ(ierr);
104 if (fd_jacobian_coloring) { /* Use finite differences with coloring */
105 /* Get data structure of J */
106 PetscBool pc_diagonal;
107 ierr = PetscOptionsHasName(NULL,NULL,"-pc_diagonal",&pc_diagonal);CHKERRQ(ierr);
108 if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
109 PetscInt rstart,rend,i;
110 PetscScalar zero=0.0;
111 ierr = MatGetOwnershipRange(J,&rstart,&rend);CHKERRQ(ierr);
112 for (i=rstart; i<rend; i++) {
113 ierr = MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);CHKERRQ(ierr);
114 }
115 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
116 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
117 } else {
118 /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
119 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
120 ierr = TSSetUp(ts);CHKERRQ(ierr);
121 ierr = SNESComputeJacobianDefault(snes,x,J,J,ts);CHKERRQ(ierr);
122 }
123
124 /* create coloring context */
125 ierr = MatColoringCreate(J,&mc);CHKERRQ(ierr);
126 ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr);
127 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr);
128 ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr);
129 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr);
130 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr);
131 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);CHKERRQ(ierr);
132 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
133 ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr);
134 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr);
135 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
136 } else { /* Use finite differences (slow) */
137 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr);
138 }
139 }
140
141 /* Pick up a Petsc preconditioner */
142 /* one can always set method or preconditioner during the run time */
143 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
144 ierr = SNESGetKSP(snes,&ksp);CHKERRQ(ierr);
145 ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);
146 ierr = PCSetType(pc,PCJACOBI);CHKERRQ(ierr);
147 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
148
149 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
150 ierr = TSSetUp(ts);CHKERRQ(ierr);
151
152 /* Test TSSetPostStep() */
153 ierr = PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);CHKERRQ(ierr);
154 if (flg) {
155 ierr = TSSetPostStep(ts,PostStep);CHKERRQ(ierr);
156 }
157
158 ierr = PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);CHKERRQ(ierr);
159 for (iout=1; iout<=NOUT; iout++) {
160 ierr = TSSetMaxSteps(ts,time_steps);CHKERRQ(ierr);
161 ierr = TSSetMaxTime(ts,iout*ftime_original/NOUT);CHKERRQ(ierr);
162 ierr = TSSolve(ts,global);CHKERRQ(ierr);
163 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
164 ierr = TSSetTime(ts,ftime);CHKERRQ(ierr);
165 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
166 }
167 /* Interpolate solution at tfinal */
168 ierr = TSGetSolution(ts,&global);CHKERRQ(ierr);
169 ierr = TSInterpolate(ts,ftime_original,global);CHKERRQ(ierr);
170
171 ierr = PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);CHKERRQ(ierr);
172 if (flg) { /* print solution into a MATLAB file */
173 ierr = PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);CHKERRQ(ierr);
174 ierr = PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);CHKERRQ(ierr);
175 ierr = VecView(global,viewfile);CHKERRQ(ierr);
176 ierr = PetscViewerPopFormat(viewfile);CHKERRQ(ierr);
177 ierr = PetscViewerDestroy(&viewfile);CHKERRQ(ierr);
178 }
179
180 /* free the memories */
181 ierr = TSDestroy(&ts);CHKERRQ(ierr);
182 ierr = VecDestroy(&global);CHKERRQ(ierr);
183 ierr = VecDestroy(&x);CHKERRQ(ierr);
184 ierr = MatDestroy(&J);CHKERRQ(ierr);
185 if (fd_jacobian_coloring) {ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);}
186 ierr = PetscFinalize();
187 return ierr;
188 }
189
190 /* -------------------------------------------------------------------*/
191 /* the initial function */
f_ini(PetscReal x,PetscReal y)192 PetscReal f_ini(PetscReal x,PetscReal y)
193 {
194 PetscReal f;
195
196 f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
197 return f;
198 }
199
Initial(Vec global,void * ctx)200 PetscErrorCode Initial(Vec global,void *ctx)
201 {
202 Data *data = (Data*)ctx;
203 PetscInt m,row,col;
204 PetscReal x,y,dx,dy;
205 PetscScalar *localptr;
206 PetscInt i,mybase,myend,locsize;
207 PetscErrorCode ierr;
208
209 PetscFunctionBeginUser;
210 /* make the local copies of parameters */
211 m = data->m;
212 dx = data->dx;
213 dy = data->dy;
214
215 /* determine starting point of each processor */
216 ierr = VecGetOwnershipRange(global,&mybase,&myend);CHKERRQ(ierr);
217 ierr = VecGetLocalSize(global,&locsize);CHKERRQ(ierr);
218
219 /* Initialize the array */
220 ierr = VecGetArrayWrite(global,&localptr);CHKERRQ(ierr);
221
222 for (i=0; i<locsize; i++) {
223 row = 1+(mybase+i)-((mybase+i)/m)*m;
224 col = (mybase+i)/m+1;
225 x = dx*row;
226 y = dy*col;
227 localptr[i] = f_ini(x,y);
228 }
229
230 ierr = VecRestoreArrayWrite(global,&localptr);CHKERRQ(ierr);
231 PetscFunctionReturn(0);
232 }
233
Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void * ctx)234 PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
235 {
236 VecScatter scatter;
237 IS from,to;
238 PetscInt i,n,*idx,nsteps,maxsteps;
239 Vec tmp_vec;
240 PetscErrorCode ierr;
241 const PetscScalar *tmp;
242
243 PetscFunctionBeginUser;
244 ierr = TSGetStepNumber(ts,&nsteps);CHKERRQ(ierr);
245 /* display output at selected time steps */
246 ierr = TSGetMaxSteps(ts, &maxsteps);CHKERRQ(ierr);
247 if (nsteps % 10 != 0) PetscFunctionReturn(0);
248
249 /* Get the size of the vector */
250 ierr = VecGetSize(global,&n);CHKERRQ(ierr);
251
252 /* Set the index sets */
253 ierr = PetscMalloc1(n,&idx);CHKERRQ(ierr);
254 for (i=0; i<n; i++) idx[i]=i;
255
256 /* Create local sequential vectors */
257 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);CHKERRQ(ierr);
258
259 /* Create scatter context */
260 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
261 ierr = ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr);
262 ierr = VecScatterCreate(global,from,tmp_vec,to,&scatter);CHKERRQ(ierr);
263 ierr = VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
264 ierr = VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
265
266 ierr = VecGetArrayRead(tmp_vec,&tmp);CHKERRQ(ierr);
267 ierr = PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));CHKERRQ(ierr);
268 ierr = VecRestoreArrayRead(tmp_vec,&tmp);CHKERRQ(ierr);
269
270 ierr = PetscFree(idx);CHKERRQ(ierr);
271 ierr = ISDestroy(&from);CHKERRQ(ierr);
272 ierr = ISDestroy(&to);CHKERRQ(ierr);
273 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr);
274 ierr = VecDestroy(&tmp_vec);CHKERRQ(ierr);
275 PetscFunctionReturn(0);
276 }
277
RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void * ptr)278 PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr)
279 {
280 Data *data = (Data*)ptr;
281 PetscScalar v[5];
282 PetscInt idx[5],i,j,row;
283 PetscErrorCode ierr;
284 PetscInt m,n,mn;
285 PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
286
287 PetscFunctionBeginUser;
288 m = data->m;
289 n = data->n;
290 mn = m*n;
291 dx = data->dx;
292 dy = data->dy;
293 a = data->a;
294 epsilon = data->epsilon;
295
296 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
297 xl = 0.5*a/dx+epsilon/(dx*dx);
298 xr = -0.5*a/dx+epsilon/(dx*dx);
299 yl = 0.5*a/dy+epsilon/(dy*dy);
300 yr = -0.5*a/dy+epsilon/(dy*dy);
301
302 row = 0;
303 v[0] = xc; v[1] = xr; v[2] = yr;
304 idx[0] = 0; idx[1] = 2; idx[2] = m;
305 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);
306
307 row = m-1;
308 v[0] = 2.0*xl; v[1] = xc; v[2] = yr;
309 idx[0] = m-2; idx[1] = m-1; idx[2] = m-1+m;
310 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);
311
312 for (i=1; i<m-1; i++) {
313 row = i;
314 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = yr;
315 idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m;
316 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr);
317 }
318
319 for (j=1; j<n-1; j++) {
320 row = j*m;
321 v[0] = xc; v[1] = xr; v[2] = yl; v[3] = yr;
322 idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m;
323 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr);
324
325 row = j*m+m-1;
326 v[0] = xc; v[1] = 2.0*xl; v[2] = yl; v[3] = yr;
327 idx[0] = j*m+m-1; idx[1] = j*m+m-1-1; idx[2] = j*m+m-1-m; idx[3] = j*m+m-1+m;
328 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr);
329
330 for (i=1; i<m-1; i++) {
331 row = j*m+i;
332 v[0] = xc; v[1] = xl; v[2] = xr; v[3] = yl; v[4]=yr;
333 idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m;
334 idx[4] = j*m+i+m;
335 ierr = MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);CHKERRQ(ierr);
336 }
337 }
338
339 row = mn-m;
340 v[0] = xc; v[1] = xr; v[2] = 2.0*yl;
341 idx[0] = mn-m; idx[1] = mn-m+1; idx[2] = mn-m-m;
342 ierr = MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);
343
344 row = mn-1;
345 v[0] = xc; v[1] = 2.0*xl; v[2] = 2.0*yl;
346 idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m;
347 ierr = MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);CHKERRQ(ierr);
348
349 for (i=1; i<m-1; i++) {
350 row = mn-m+i;
351 v[0] = xl; v[1] = xc; v[2] = xr; v[3] = 2.0*yl;
352 idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m;
353 ierr = MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);CHKERRQ(ierr);
354 }
355
356 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
357 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
358
359 PetscFunctionReturn(0);
360 }
361
362 /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void * ctx)363 PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
364 {
365 Data *data = (Data*)ctx;
366 PetscInt m,n,mn;
367 PetscReal dx,dy;
368 PetscReal xc,xl,xr,yl,yr;
369 PetscReal a,epsilon;
370 PetscScalar *outptr;
371 const PetscScalar *inptr;
372 PetscInt i,j,len;
373 PetscErrorCode ierr;
374 IS from,to;
375 PetscInt *idx;
376 VecScatter scatter;
377 Vec tmp_in,tmp_out;
378
379 PetscFunctionBeginUser;
380 m = data->m;
381 n = data->n;
382 mn = m*n;
383 dx = data->dx;
384 dy = data->dy;
385 a = data->a;
386 epsilon = data->epsilon;
387
388 xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
389 xl = 0.5*a/dx+epsilon/(dx*dx);
390 xr = -0.5*a/dx+epsilon/(dx*dx);
391 yl = 0.5*a/dy+epsilon/(dy*dy);
392 yr = -0.5*a/dy+epsilon/(dy*dy);
393
394 /* Get the length of parallel vector */
395 ierr = VecGetSize(globalin,&len);CHKERRQ(ierr);
396
397 /* Set the index sets */
398 ierr = PetscMalloc1(len,&idx);CHKERRQ(ierr);
399 for (i=0; i<len; i++) idx[i]=i;
400
401 /* Create local sequential vectors */
402 ierr = VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);CHKERRQ(ierr);
403 ierr = VecDuplicate(tmp_in,&tmp_out);CHKERRQ(ierr);
404
405 /* Create scatter context */
406 ierr = ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);CHKERRQ(ierr);
407 ierr = ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);CHKERRQ(ierr);
408 ierr = VecScatterCreate(globalin,from,tmp_in,to,&scatter);CHKERRQ(ierr);
409 ierr = VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
410 ierr = VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
411 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr);
412
413 /*Extract income array - include ghost points */
414 ierr = VecGetArrayRead(tmp_in,&inptr);CHKERRQ(ierr);
415
416 /* Extract outcome array*/
417 ierr = VecGetArrayWrite(tmp_out,&outptr);CHKERRQ(ierr);
418
419 outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
420 outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
421 for (i=1; i<m-1; i++) {
422 outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m];
423 }
424
425 for (j=1; j<n-1; j++) {
426 outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m];
427 outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+ yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
428 for (i=1; i<m-1; i++) {
429 outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]+yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
430 }
431 }
432
433 outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
434 outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
435 for (i=1; i<m-1; i++) {
436 outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]+2*yl*inptr[mn-m+i-m];
437 }
438
439 ierr = VecRestoreArrayRead(tmp_in,&inptr);CHKERRQ(ierr);
440 ierr = VecRestoreArrayWrite(tmp_out,&outptr);CHKERRQ(ierr);
441
442 ierr = VecScatterCreate(tmp_out,from,globalout,to,&scatter);CHKERRQ(ierr);
443 ierr = VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
444 ierr = VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
445
446 /* Destroy idx aand scatter */
447 ierr = VecDestroy(&tmp_in);CHKERRQ(ierr);
448 ierr = VecDestroy(&tmp_out);CHKERRQ(ierr);
449 ierr = ISDestroy(&from);CHKERRQ(ierr);
450 ierr = ISDestroy(&to);CHKERRQ(ierr);
451 ierr = VecScatterDestroy(&scatter);CHKERRQ(ierr);
452
453 ierr = PetscFree(idx);CHKERRQ(ierr);
454 PetscFunctionReturn(0);
455 }
456
PostStep(TS ts)457 PetscErrorCode PostStep(TS ts)
458 {
459 PetscErrorCode ierr;
460 PetscReal t;
461
462 PetscFunctionBeginUser;
463 ierr = TSGetTime(ts,&t);CHKERRQ(ierr);
464 ierr = PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t);CHKERRQ(ierr);
465 PetscFunctionReturn(0);
466 }
467
468 /*TEST
469
470 test:
471 args: -ts_fd -ts_type beuler
472 output_file: output/ex4.out
473
474 test:
475 suffix: 2
476 args: -ts_fd -ts_type beuler
477 nsize: 2
478 output_file: output/ex4.out
479
480 test:
481 suffix: 3
482 args: -ts_fd -ts_type cn
483
484 test:
485 suffix: 4
486 args: -ts_fd -ts_type cn
487 output_file: output/ex4_3.out
488 nsize: 2
489
490 test:
491 suffix: 5
492 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
493 output_file: output/ex4.out
494
495 test:
496 suffix: 6
497 args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
498 output_file: output/ex4.out
499 nsize: 2
500
501 test:
502 suffix: 7
503 requires: !single
504 args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1
505
506 test:
507 suffix: 8
508 requires: !single
509 args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view
510
511 TEST*/
512
513