Actual source code: ex4.c
petsc-3.8.4 2018-03-24
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 with SUNDIALS.
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;
58: PetscViewer viewer;
59: char pcinfo[120],tsinfo[120];
60: TSType tstype;
61: PetscBool sundials;
63: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
64: MPI_Comm_size(PETSC_COMM_WORLD,&size);
66: /* set data */
67: data.m = 9;
68: data.n = 9;
69: data.a = 1.0;
70: data.epsilon = 0.1;
71: data.dx = 1.0/(data.m+1.0);
72: data.dy = 1.0/(data.n+1.0);
73: mn = (data.m)*(data.n);
74: PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);
76: /* set initial conditions */
77: VecCreate(PETSC_COMM_WORLD,&global);
78: VecSetSizes(global,PETSC_DECIDE,mn);
79: VecSetFromOptions(global);
80: Initial(global,&data);
81: VecDuplicate(global,&x);
83: /* create timestep context */
84: TSCreate(PETSC_COMM_WORLD,&ts);
85: TSMonitorSet(ts,Monitor,&data,NULL);
86: #if defined(PETSC_HAVE_SUNDIALS)
87: TSSetType(ts,TSSUNDIALS);
88: #else
89: TSSetType(ts,TSEULER);
90: #endif
91: dt = 0.1;
92: ftime_original = data.tfinal = 1.0;
94: TSSetTimeStep(ts,dt);
95: TSSetMaxSteps(ts,time_steps);
96: TSSetMaxTime(ts,ftime_original);
97: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
98: TSSetSolution(ts,global);
100: /* set user provided RHSFunction and RHSJacobian */
101: TSSetRHSFunction(ts,NULL,RHSFunction,&data);
102: MatCreate(PETSC_COMM_WORLD,&J);
103: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
104: MatSetFromOptions(J);
105: MatSeqAIJSetPreallocation(J,5,NULL);
106: MatMPIAIJSetPreallocation(J,5,NULL,5,NULL);
108: PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);
109: if (!flg) {
110: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
111: } else {
112: TSGetSNES(ts,&snes);
113: PetscOptionsHasName(NULL,NULL,"-fd_color",&fd_jacobian_coloring);
114: if (fd_jacobian_coloring) { /* Use finite differences with coloring */
115: /* Get data structure of J */
116: PetscBool pc_diagonal;
117: PetscOptionsHasName(NULL,NULL,"-pc_diagonal",&pc_diagonal);
118: if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
119: PetscInt rstart,rend,i;
120: PetscScalar zero=0.0;
121: MatGetOwnershipRange(J,&rstart,&rend);
122: for (i=rstart; i<rend; i++) {
123: MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
124: }
125: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
127: } else {
128: /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
129: TSSetType(ts,TSBEULER);
130: TSSetUp(ts);
131: SNESComputeJacobianDefault(snes,x,J,J,ts);
132: }
134: /* create coloring context */
135: MatColoringCreate(J,&mc);
136: MatColoringSetType(mc,MATCOLORINGSL);
137: MatColoringSetFromOptions(mc);
138: MatColoringApply(mc,&iscoloring);
139: MatColoringDestroy(&mc);
140: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
141: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
142: MatFDColoringSetFromOptions(matfdcoloring);
143: MatFDColoringSetUp(J,iscoloring,matfdcoloring);
144: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
145: ISColoringDestroy(&iscoloring);
146: } else { /* Use finite differences (slow) */
147: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);
148: }
149: }
151: /* Pick up a Petsc preconditioner */
152: /* one can always set method or preconditioner during the run time */
153: TSGetSNES(ts,&snes);
154: SNESGetKSP(snes,&ksp);
155: KSPGetPC(ksp,&pc);
156: PCSetType(pc,PCJACOBI);
157: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
159: TSSetFromOptions(ts);
160: TSSetUp(ts);
162: /* Test TSSetPostStep() */
163: PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);
164: if (flg) {
165: TSSetPostStep(ts,PostStep);
166: }
168: PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);
169: for (iout=1; iout<=NOUT; iout++) {
170: TSSetMaxSteps(ts,time_steps);
171: TSSetMaxTime(ts,iout*ftime_original/NOUT);
172: TSSolve(ts,global);
173: TSGetSolveTime(ts,&ftime);
174: TSSetTime(ts,ftime);
175: TSSetTimeStep(ts,dt);
176: }
177: /* Interpolate solution at tfinal */
178: TSGetSolution(ts,&global);
179: TSInterpolate(ts,ftime_original,global);
181: PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);
182: if (flg) { /* print solution into a MATLAB file */
183: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
184: PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
185: VecView(global,viewfile);
186: PetscViewerPopFormat(viewfile);
187: PetscViewerDestroy(&viewfile);
188: }
190: /* display solver info for Sundials */
191: TSGetType(ts,&tstype);
192: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
193: if (sundials) {
194: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
195: TSView(ts,viewer);
196: PetscViewerDestroy(&viewer);
197: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
198: PCView(pc,viewer);
199: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",size,tsinfo,pcinfo);
200: PetscViewerDestroy(&viewer);
201: }
203: /* free the memories */
204: TSDestroy(&ts);
205: VecDestroy(&global);
206: VecDestroy(&x);
207: MatDestroy(&J);
208: if (fd_jacobian_coloring) {MatFDColoringDestroy(&matfdcoloring);}
209: PetscFinalize();
210: return ierr;
211: }
213: /* -------------------------------------------------------------------*/
214: /* the initial function */
215: PetscReal f_ini(PetscReal x,PetscReal y)
216: {
217: PetscReal f;
219: f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
220: return f;
221: }
223: PetscErrorCode Initial(Vec global,void *ctx)
224: {
225: Data *data = (Data*)ctx;
226: PetscInt m,row,col;
227: PetscReal x,y,dx,dy;
228: PetscScalar *localptr;
229: PetscInt i,mybase,myend,locsize;
233: /* make the local copies of parameters */
234: m = data->m;
235: dx = data->dx;
236: dy = data->dy;
238: /* determine starting point of each processor */
239: VecGetOwnershipRange(global,&mybase,&myend);
240: VecGetLocalSize(global,&locsize);
242: /* Initialize the array */
243: VecGetArray(global,&localptr);
245: for (i=0; i<locsize; i++) {
246: row = 1+(mybase+i)-((mybase+i)/m)*m;
247: col = (mybase+i)/m+1;
248: x = dx*row;
249: y = dy*col;
250: localptr[i] = f_ini(x,y);
251: }
253: VecRestoreArray(global,&localptr);
254: return(0);
255: }
257: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
258: {
259: VecScatter scatter;
260: IS from,to;
261: PetscInt i,n,*idx,nsteps,maxsteps;
262: Vec tmp_vec;
263: PetscErrorCode ierr;
264: const PetscScalar *tmp;
267: TSGetStepNumber(ts,&nsteps);
268: /* display output at selected time steps */
269: TSGetMaxSteps(ts, &maxsteps);
270: if (nsteps % 10 != 0) return(0);
272: /* Get the size of the vector */
273: VecGetSize(global,&n);
275: /* Set the index sets */
276: PetscMalloc1(n,&idx);
277: for (i=0; i<n; i++) idx[i]=i;
279: /* Create local sequential vectors */
280: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
282: /* Create scatter context */
283: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
284: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
285: VecScatterCreate(global,from,tmp_vec,to,&scatter);
286: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
287: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
289: VecGetArrayRead(tmp_vec,&tmp);
290: PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));
291: VecRestoreArrayRead(tmp_vec,&tmp);
293: PetscFree(idx);
294: ISDestroy(&from);
295: ISDestroy(&to);
296: VecScatterDestroy(&scatter);
297: VecDestroy(&tmp_vec);
298: return(0);
299: }
301: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr)
302: {
303: Data *data = (Data*)ptr;
304: PetscScalar v[5];
305: PetscInt idx[5],i,j,row;
307: PetscInt m,n,mn;
308: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
311: m = data->m;
312: n = data->n;
313: mn = m*n;
314: dx = data->dx;
315: dy = data->dy;
316: a = data->a;
317: epsilon = data->epsilon;
319: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
320: xl = 0.5*a/dx+epsilon/(dx*dx);
321: xr = -0.5*a/dx+epsilon/(dx*dx);
322: yl = 0.5*a/dy+epsilon/(dy*dy);
323: yr = -0.5*a/dy+epsilon/(dy*dy);
325: row = 0;
326: v[0] = xc; v[1] = xr; v[2] = yr;
327: idx[0] = 0; idx[1] = 2; idx[2] = m;
328: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
330: row = m-1;
331: v[0] = 2.0*xl; v[1] = xc; v[2] = yr;
332: idx[0] = m-2; idx[1] = m-1; idx[2] = m-1+m;
333: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
335: for (i=1; i<m-1; i++) {
336: row = i;
337: v[0] = xl; v[1] = xc; v[2] = xr; v[3] = yr;
338: idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m;
339: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
340: }
342: for (j=1; j<n-1; j++) {
343: row = j*m;
344: v[0] = xc; v[1] = xr; v[2] = yl; v[3] = yr;
345: idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m;
346: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
348: row = j*m+m-1;
349: v[0] = xc; v[1] = 2.0*xl; v[2] = yl; v[3] = yr;
350: 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;
351: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
353: for (i=1; i<m-1; i++) {
354: row = j*m+i;
355: v[0] = xc; v[1] = xl; v[2] = xr; v[3] = yl; v[4]=yr;
356: idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m;
357: idx[4] = j*m+i+m;
358: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
359: }
360: }
362: row = mn-m;
363: v[0] = xc; v[1] = xr; v[2] = 2.0*yl;
364: idx[0] = mn-m; idx[1] = mn-m+1; idx[2] = mn-m-m;
365: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
367: row = mn-1;
368: v[0] = xc; v[1] = 2.0*xl; v[2] = 2.0*yl;
369: idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m;
370: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
372: for (i=1; i<m-1; i++) {
373: row = mn-m+i;
374: v[0] = xl; v[1] = xc; v[2] = xr; v[3] = 2.0*yl;
375: idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m;
376: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
377: }
379: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
380: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
382: return(0);
383: }
385: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
386: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
387: {
388: Data *data = (Data*)ctx;
389: PetscInt m,n,mn;
390: PetscReal dx,dy;
391: PetscReal xc,xl,xr,yl,yr;
392: PetscReal a,epsilon;
393: PetscScalar *outptr;
394: const PetscScalar *inptr;
395: PetscInt i,j,len;
396: PetscErrorCode ierr;
397: IS from,to;
398: PetscInt *idx;
399: VecScatter scatter;
400: Vec tmp_in,tmp_out;
403: m = data->m;
404: n = data->n;
405: mn = m*n;
406: dx = data->dx;
407: dy = data->dy;
408: a = data->a;
409: epsilon = data->epsilon;
411: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
412: xl = 0.5*a/dx+epsilon/(dx*dx);
413: xr = -0.5*a/dx+epsilon/(dx*dx);
414: yl = 0.5*a/dy+epsilon/(dy*dy);
415: yr = -0.5*a/dy+epsilon/(dy*dy);
417: /* Get the length of parallel vector */
418: VecGetSize(globalin,&len);
420: /* Set the index sets */
421: PetscMalloc1(len,&idx);
422: for (i=0; i<len; i++) idx[i]=i;
424: /* Create local sequential vectors */
425: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
426: VecDuplicate(tmp_in,&tmp_out);
428: /* Create scatter context */
429: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
430: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
431: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
432: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
433: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
434: VecScatterDestroy(&scatter);
436: /*Extract income array - include ghost points */
437: VecGetArrayRead(tmp_in,&inptr);
439: /* Extract outcome array*/
440: VecGetArray(tmp_out,&outptr);
442: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
443: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
444: for (i=1; i<m-1; i++) {
445: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m];
446: }
448: for (j=1; j<n-1; j++) {
449: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m];
450: 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];
451: for (i=1; i<m-1; i++) {
452: 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];
453: }
454: }
456: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
457: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
458: for (i=1; i<m-1; i++) {
459: 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];
460: }
462: VecRestoreArrayRead(tmp_in,&inptr);
463: VecRestoreArray(tmp_out,&outptr);
465: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
466: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
467: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
469: /* Destroy idx aand scatter */
470: VecDestroy(&tmp_in);
471: VecDestroy(&tmp_out);
472: ISDestroy(&from);
473: ISDestroy(&to);
474: VecScatterDestroy(&scatter);
476: PetscFree(idx);
477: return(0);
478: }
480: PetscErrorCode PostStep(TS ts)
481: {
483: PetscReal t;
486: TSGetTime(ts,&t);
487: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t);
488: return(0);
489: }