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