Actual source code: ex4.c
petsc-3.14.6 2021-03-30
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: 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;
58: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
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: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
70: /* set initial conditions */
71: VecCreate(PETSC_COMM_WORLD,&global);
72: VecSetSizes(global,PETSC_DECIDE,mn);
73: VecSetFromOptions(global);
74: Initial(global,&data);
75: VecDuplicate(global,&x);
77: /* create timestep context */
78: TSCreate(PETSC_COMM_WORLD,&ts);
79: TSMonitorSet(ts,Monitor,&data,NULL);
80: TSSetType(ts,TSEULER);
81: dt = 0.1;
82: ftime_original = data.tfinal = 1.0;
84: TSSetTimeStep(ts,dt);
85: TSSetMaxSteps(ts,time_steps);
86: TSSetMaxTime(ts,ftime_original);
87: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
88: TSSetSolution(ts,global);
90: /* set user provided RHSFunction and RHSJacobian */
91: TSSetRHSFunction(ts,NULL,RHSFunction,&data);
92: MatCreate(PETSC_COMM_WORLD,&J);
93: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
94: MatSetFromOptions(J);
95: MatSeqAIJSetPreallocation(J,5,NULL);
96: MatMPIAIJSetPreallocation(J,5,NULL,5,NULL);
98: PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);
99: if (!flg) {
100: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
101: } else {
102: TSGetSNES(ts,&snes);
103: PetscOptionsHasName(NULL,NULL,"-fd_color",&fd_jacobian_coloring);
104: if (fd_jacobian_coloring) { /* Use finite differences with coloring */
105: /* Get data structure of J */
106: PetscBool pc_diagonal;
107: PetscOptionsHasName(NULL,NULL,"-pc_diagonal",&pc_diagonal);
108: if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
109: PetscInt rstart,rend,i;
110: PetscScalar zero=0.0;
111: MatGetOwnershipRange(J,&rstart,&rend);
112: for (i=rstart; i<rend; i++) {
113: MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
114: }
115: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
117: } else {
118: /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
119: TSSetType(ts,TSBEULER);
120: TSSetUp(ts);
121: SNESComputeJacobianDefault(snes,x,J,J,ts);
122: }
124: /* create coloring context */
125: MatColoringCreate(J,&mc);
126: MatColoringSetType(mc,MATCOLORINGSL);
127: MatColoringSetFromOptions(mc);
128: MatColoringApply(mc,&iscoloring);
129: MatColoringDestroy(&mc);
130: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
131: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
132: MatFDColoringSetFromOptions(matfdcoloring);
133: MatFDColoringSetUp(J,iscoloring,matfdcoloring);
134: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
135: ISColoringDestroy(&iscoloring);
136: } else { /* Use finite differences (slow) */
137: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);
138: }
139: }
141: /* Pick up a Petsc preconditioner */
142: /* one can always set method or preconditioner during the run time */
143: TSGetSNES(ts,&snes);
144: SNESGetKSP(snes,&ksp);
145: KSPGetPC(ksp,&pc);
146: PCSetType(pc,PCJACOBI);
147: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
149: TSSetFromOptions(ts);
150: TSSetUp(ts);
152: /* Test TSSetPostStep() */
153: PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);
154: if (flg) {
155: TSSetPostStep(ts,PostStep);
156: }
158: PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);
159: for (iout=1; iout<=NOUT; iout++) {
160: TSSetMaxSteps(ts,time_steps);
161: TSSetMaxTime(ts,iout*ftime_original/NOUT);
162: TSSolve(ts,global);
163: TSGetSolveTime(ts,&ftime);
164: TSSetTime(ts,ftime);
165: TSSetTimeStep(ts,dt);
166: }
167: /* Interpolate solution at tfinal */
168: TSGetSolution(ts,&global);
169: TSInterpolate(ts,ftime_original,global);
171: PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);
172: if (flg) { /* print solution into a MATLAB file */
173: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
174: PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
175: VecView(global,viewfile);
176: PetscViewerPopFormat(viewfile);
177: PetscViewerDestroy(&viewfile);
178: }
180: /* free the memories */
181: TSDestroy(&ts);
182: VecDestroy(&global);
183: VecDestroy(&x);
184: MatDestroy(&J);
185: if (fd_jacobian_coloring) {MatFDColoringDestroy(&matfdcoloring);}
186: PetscFinalize();
187: return ierr;
188: }
190: /* -------------------------------------------------------------------*/
191: /* the initial function */
192: PetscReal f_ini(PetscReal x,PetscReal y)
193: {
194: PetscReal f;
196: f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
197: return f;
198: }
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;
210: /* make the local copies of parameters */
211: m = data->m;
212: dx = data->dx;
213: dy = data->dy;
215: /* determine starting point of each processor */
216: VecGetOwnershipRange(global,&mybase,&myend);
217: VecGetLocalSize(global,&locsize);
219: /* Initialize the array */
220: VecGetArrayWrite(global,&localptr);
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: }
230: VecRestoreArrayWrite(global,&localptr);
231: return(0);
232: }
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;
244: TSGetStepNumber(ts,&nsteps);
245: /* display output at selected time steps */
246: TSGetMaxSteps(ts, &maxsteps);
247: if (nsteps % 10 != 0) return(0);
249: /* Get the size of the vector */
250: VecGetSize(global,&n);
252: /* Set the index sets */
253: PetscMalloc1(n,&idx);
254: for (i=0; i<n; i++) idx[i]=i;
256: /* Create local sequential vectors */
257: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
259: /* Create scatter context */
260: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
261: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
262: VecScatterCreate(global,from,tmp_vec,to,&scatter);
263: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
264: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
266: VecGetArrayRead(tmp_vec,&tmp);
267: PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));
268: VecRestoreArrayRead(tmp_vec,&tmp);
270: PetscFree(idx);
271: ISDestroy(&from);
272: ISDestroy(&to);
273: VecScatterDestroy(&scatter);
274: VecDestroy(&tmp_vec);
275: return(0);
276: }
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;
284: PetscInt m,n,mn;
285: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
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;
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);
302: row = 0;
303: v[0] = xc; v[1] = xr; v[2] = yr;
304: idx[0] = 0; idx[1] = 2; idx[2] = m;
305: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
317: }
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: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
336: }
337: }
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: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
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: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
354: }
356: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
357: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
359: return(0);
360: }
362: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
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;
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;
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);
394: /* Get the length of parallel vector */
395: VecGetSize(globalin,&len);
397: /* Set the index sets */
398: PetscMalloc1(len,&idx);
399: for (i=0; i<len; i++) idx[i]=i;
401: /* Create local sequential vectors */
402: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
403: VecDuplicate(tmp_in,&tmp_out);
405: /* Create scatter context */
406: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
407: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
408: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
409: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
410: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
411: VecScatterDestroy(&scatter);
413: /*Extract income array - include ghost points */
414: VecGetArrayRead(tmp_in,&inptr);
416: /* Extract outcome array*/
417: VecGetArrayWrite(tmp_out,&outptr);
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: }
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: }
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: }
439: VecRestoreArrayRead(tmp_in,&inptr);
440: VecRestoreArrayWrite(tmp_out,&outptr);
442: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
443: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
444: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
446: /* Destroy idx aand scatter */
447: VecDestroy(&tmp_in);
448: VecDestroy(&tmp_out);
449: ISDestroy(&from);
450: ISDestroy(&to);
451: VecScatterDestroy(&scatter);
453: PetscFree(idx);
454: return(0);
455: }
457: PetscErrorCode PostStep(TS ts)
458: {
460: PetscReal t;
463: TSGetTime(ts,&t);
464: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t);
465: return(0);
466: }
468: /*TEST
470: test:
471: args: -ts_fd -ts_type beuler
472: output_file: output/ex4.out
474: test:
475: suffix: 2
476: args: -ts_fd -ts_type beuler
477: nsize: 2
478: output_file: output/ex4.out
480: test:
481: suffix: 3
482: args: -ts_fd -ts_type cn
484: test:
485: suffix: 4
486: args: -ts_fd -ts_type cn
487: output_file: output/ex4_3.out
488: nsize: 2
490: test:
491: suffix: 5
492: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
493: output_file: output/ex4.out
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
501: test:
502: suffix: 7
503: requires: !single
504: args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1
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
511: TEST*/