Actual source code: ex50.c
1: /* DMMG/KSP solving a system of linear equations.
2: Poisson equation in 2D:
4: div(grad p) = f, 0 < x,y < 1
5: with
6: forcing function f = -cos(m*pi*x)*cos(n*pi*y),
7: Neuman boundary conditions
8: dp/dx = 0 for x = 0, x = 1.
9: dp/dy = 0 for y = 0, y = 1.
11: Contributed by Michael Boghosian <boghmic@iit.edu>, 2008,
12: based on petsc/src/ksp/ksp/examples/tutorials/ex29.c and ex32.c
14: Example of Usage:
15: ./ex50 -mglevels 3 -ksp_monitor -M 3 -N 3 -ksp_view -da_view_draw -draw_pause -1
16: ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_factor_levels <ilu_levels> -ksp_monitor -cmp_solu
17: ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor -cmp_solu
18: mpiexec -n 4 ./ex50 -M 3 -N 3 -ksp_monitor -ksp_view -mglevels 10 -log_summary
19: */
21: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";
23: #include <petscdmda.h>
24: #include <petscksp.h>
25: #include <petscpcmg.h>
26: #include <petscdmmg.h>
27: #include <petscsys.h>
28: #include <petscvec.h>
35: typedef enum {DIRICHLET, NEUMANN} BCType;
37: typedef struct {
38: PetscScalar uu, tt;
39: BCType bcType;
40: } UserContext;
44: int main(int argc,char **argv)
45: {
46: DMMG *dmmg;
47: DM da;
48: UserContext user;
49: PetscInt l, bc, mglevels, M, N, stages[3];
50: PetscReal norm;
52: PetscMPIInt rank,nproc;
53: PetscBool flg;
54:
55: PetscInitialize(&argc,&argv,(char *)0,help);
56: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
57: MPI_Comm_rank(PETSC_COMM_WORLD,&nproc);
58: PetscLogStageRegister("DMMG Setup",&stages[0]);
59: PetscLogStageRegister("DMMG Solve",&stages[1]);
60:
61: PetscLogStagePush(stages[0]); /* Start DMMG Setup */
62:
63: /* SET VARIABLES: */
64: mglevels = 1;
65: PetscOptionsGetInt(PETSC_NULL,"-mglevels",&mglevels,PETSC_NULL);
66: M = 11; /* number of grid points in x dir. on coarse grid */
67: N = 11; /* number of grid points in y dir. on coarse grid */
68: PetscOptionsGetInt(PETSC_NULL,"-M",&M,PETSC_NULL);
69: PetscOptionsGetInt(PETSC_NULL,"-N",&N,PETSC_NULL);
70:
71: DMMGCreate(PETSC_COMM_WORLD,mglevels,PETSC_NULL,&dmmg);
72: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,M,N,PETSC_DECIDE,PETSC_DECIDE,1,1,PETSC_NULL,PETSC_NULL,&da);
73: DMMGSetDM(dmmg,(DM)da);
74: DMDestroy(&da);
75:
76: /* Set user contex */
77: user.uu = 1.0;
78: user.tt = 1.0;
79: bc = (PetscInt)NEUMANN; // Use Neumann Boundary Conditions
80: user.bcType = (BCType)bc;
81: for (l = 0; l < DMMGGetLevels(dmmg); l++) {
82: DMMGSetUser(dmmg,l,&user);
83: }
84:
85: DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);
86: if (user.bcType == NEUMANN){
87: DMMGSetNullSpace(dmmg,PETSC_TRUE,0,PETSC_NULL);
88: }
89: PetscLogStagePop(); /* Finish DMMG Setup */
91: /* DMMG SOLVE: */
92: PetscLogStagePush(stages[1]); /* Start DMMG Solve */
93: DMMGSolve(dmmg);
94: PetscLogStagePop(); /* Finish DMMG Solve */
96: /* Compare solution with the true solution p */
97: PetscOptionsHasName(PETSC_NULL, "-cmp_solu", &flg);
98: if (flg){
99: Vec x,p;
100: if (mglevels != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"mglevels must equls 1 for comparison");
101: if (nproc > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"num of proc must equls 1 for comparison");
102: x = DMMGGetx(dmmg);
103: VecDuplicate(x,&p);
104: ComputeTrueSolution(dmmg,p);
106: VecAXPY(p,-1.0,x); /* p <- (p-x) */
107: VecNorm(p, NORM_2, &norm);
108: if (!rank){printf("| solu_compt - solu_true | = %g\n",norm);}
109: VecDestroy(&p);
110: }
111:
112: DMMGDestroy(dmmg);
113: PetscFinalize();
114: return 0;
115: }
117: // COMPUTE RHS:--------------------------------------------------------------
120: PetscErrorCode ComputeRHS(DMMG dmmg, Vec b)
121: {
122: DM da = dmmg->dm;
123: UserContext *user = (UserContext *) dmmg->user;
125: PetscInt i, j, M, N, xm ,ym ,xs, ys;
126: PetscScalar Hx, Hy, pi, uu, tt;
127: PetscScalar **array;
130: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
131: uu = user->uu; tt = user->tt;
132: pi = 4*atan(1.0);
133: Hx = 1.0/(PetscReal)(M);
134: Hy = 1.0/(PetscReal)(N);
136: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); // Fine grid
137: //printf(" M N: %d %d; xm ym: %d %d; xs ys: %d %d\n",M,N,xm,ym,xs,ys);
138: DMDAVecGetArray(da, b, &array);
139: for (j=ys; j<ys+ym; j++){
140: for(i=xs; i<xs+xm; i++){
141: array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*cos(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
142: }
143: }
144: DMDAVecRestoreArray(da, b, &array);
145: VecAssemblyBegin(b);
146: VecAssemblyEnd(b);
148: /* force right hand side to be consistent for singular matrix */
149: /* note this is really a hack, normally the model would provide you with a consistent right handside */
150: if (user->bcType == NEUMANN) {
151: MatNullSpace nullspace;
152: KSPGetNullSpace(dmmg->ksp,&nullspace);
153: MatNullSpaceRemove(nullspace,b,PETSC_NULL);
154: }
155: return(0);
156: }
158: // COMPUTE JACOBIAN:--------------------------------------------------------------
161: PetscErrorCode ComputeJacobian(DMMG dmmg, Mat J, Mat jac)
162: {
163: DM da = dmmg->dm;
164: UserContext *user = (UserContext *) dmmg->user;
166: PetscInt i, j, M, N, xm, ym, xs, ys, num, numi, numj;
167: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
168: MatStencil row, col[5];
171: DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
172: Hx = 1.0 / (PetscReal)(M);
173: Hy = 1.0 / (PetscReal)(N);
174: HxdHy = Hx/Hy;
175: HydHx = Hy/Hx;
176: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
177: for (j=ys; j<ys+ym; j++){
178: for(i=xs; i<xs+xm; i++){
179: row.i = i; row.j = j;
180:
181: if (i==0 || j==0 || i==M-1 || j==N-1) {
182: if (user->bcType == DIRICHLET){
183: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
184: } else if (user->bcType == NEUMANN){
185: num=0; numi=0; numj=0;
186: if (j!=0) {
187: v[num] = -HxdHy; col[num].i = i; col[num].j = j-1;
188: num++; numj++;
189: }
190: if (i!=0) {
191: v[num] = -HydHx; col[num].i = i-1; col[num].j = j;
192: num++; numi++;
193: }
194: if (i!=M-1) {
195: v[num] = -HydHx; col[num].i = i+1; col[num].j = j;
196: num++; numi++;
197: }
198: if (j!=N-1) {
199: v[num] = -HxdHy; col[num].i = i; col[num].j = j+1;
200: num++; numj++;
201: }
202: v[num] = ( (PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx ); col[num].i = i; col[num].j = j;
203: num++;
204: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
205: }
206: } else {
207: v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
208: v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
209: v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
210: v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
211: v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
212: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
213: }
214: }
215: }
216: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
217: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
218: return(0);
219: }
221: // COMPUTE TrueSolution:--------------------------------------------------------------
224: PetscErrorCode ComputeTrueSolution(DMMG *dmmg, Vec b)
225: {
226: DM da = (*dmmg)->dm;
227: UserContext *user = (UserContext *) (*dmmg)->user;
229: PetscInt i, j, M, N, xm ,ym ,xs, ys;
230: PetscScalar Hx, Hy, pi, uu, tt, cc;
231: PetscScalar *array;
234: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0); /* level_0 ! */
235: //printf("ComputeTrueSolution - M N: %d %d;\n",M,N);
237: uu = user->uu; tt = user->tt;
238: pi = 4*atan(1.0);
239: cc = -1.0/( (uu*pi)*(uu*pi) + (tt*pi)*(tt*pi) );
240: Hx = 1.0/(PetscReal)(M);
241: Hy = 1.0/(PetscReal)(N);
243: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
244: VecGetArray(b, &array);
245: for (j=ys; j<ys+ym; j++){
246: for(i=xs; i<xs+xm; i++){
247: array[j*xm+i] = cos(uu*pi*((PetscReal)i+0.5)*Hx)*cos(tt*pi*((PetscReal)j+0.5)*Hy)*cc;
248: }
249: }
250: VecRestoreArray(b, &array);
251: return(0);
252: }
254: // VECVIEW_VTK:--------------------------------------------------------------
257: PetscErrorCode VecView_VTK(Vec x, const char filename[], const char bcName[])
258: {
259: MPI_Comm comm;
260: DM da;
261: Vec coords;
262: PetscViewer viewer;
263: PetscScalar *array, *values;
264: PetscInt nn, NN, maxn, M, N;
265: PetscInt i, p, dof = 1;
266: MPI_Status status;
267: PetscMPIInt rank, size, tag;
268: PetscErrorCode ierr;
269: PetscBool flg;
272: PetscObjectGetComm((PetscObject) x, &comm);
273: PetscViewerASCIIOpen(comm, filename, &viewer);
275: VecGetSize(x, &NN);
276: VecGetLocalSize(x, &nn);
277: PetscObjectQuery((PetscObject) x, "DM", (PetscObject *) &da);
278: if (!da) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
279: PetscTypeCompare((PetscObject)da,DMDA,&flg);
280: if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_WRONG,"Vector not generated from a DMDA");
282: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,&dof,0,0,0,0,0);
284: PetscViewerASCIIPrintf(viewer, "# vtk DataFile Version 2.0\n");
285: PetscViewerASCIIPrintf(viewer, "Inhomogeneous Poisson Equation with %s boundary conditions\n", bcName);
286: PetscViewerASCIIPrintf(viewer, "ASCII\n");
287: // get coordinates of nodes
288: DMDAGetCoordinates(da, &coords);
289: if (!coords) {
290: DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0);
291: DMDAGetCoordinates(da, &coords);
292: }
293: PetscViewerASCIIPrintf(viewer, "DATASET RECTILINEAR_GRID\n");
294: PetscViewerASCIIPrintf(viewer, "DIMENSIONS %d %d %d\n", M, N, 1);
295: VecGetArray(coords, &array);
296: PetscViewerASCIIPrintf(viewer, "X_COORDINATES %d double\n", M);
297: for(i = 0; i < M; i++) {
298: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*2]));
299: }
300: PetscViewerASCIIPrintf(viewer, "\n");
301: PetscViewerASCIIPrintf(viewer, "Y_COORDINATES %d double\n", N);
302: for(i = 0; i < N; i++) {
303: PetscViewerASCIIPrintf(viewer, "%G ", PetscRealPart(array[i*M*2+1]));
304: }
305: PetscViewerASCIIPrintf(viewer, "\n");
306: PetscViewerASCIIPrintf(viewer, "Z_COORDINATES %d double\n", 1);
307: PetscViewerASCIIPrintf(viewer, "%G\n", 0.0);
308: VecRestoreArray(coords, &array);
309: PetscViewerASCIIPrintf(viewer, "POINT_DATA %d\n", N);
310: PetscViewerASCIIPrintf(viewer, "SCALARS scalars double %d\n", 1);
311: PetscViewerASCIIPrintf(viewer, "LOOKUP_TABLE default\n");
312: VecGetArray(x, &array);
314: // Determine maximum message to arrive:
315: MPI_Comm_rank(comm, &rank);
316: MPI_Comm_size(comm, &size) ;
317: MPI_Reduce(&nn, &maxn, 1, MPIU_INT, MPI_MAX, 0, comm);
318: tag = ((PetscObject) viewer)->tag;
319: if (!rank) {
320: PetscMalloc((maxn+1) * sizeof(PetscScalar), &values);
321: for(i = 0; i < nn; i++) {
322: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
323: }
324: for(p = 1; p < size; p++) {
325: MPI_Recv(values, (PetscMPIInt) nn, MPIU_SCALAR, p, tag, comm, &status);
326: MPI_Get_count(&status, MPIU_SCALAR, &nn);
327: for(i = 0; i < nn; i++) {
328: PetscViewerASCIIPrintf(viewer, "%G\n", PetscRealPart(array[i]));
329: }
330: }
331: PetscFree(values);
332: } else {
333: MPI_Send(array, nn, MPIU_SCALAR, 0, tag, comm);
334: }
335: VecRestoreArray(x, &array);
336: PetscViewerFlush(viewer);
337: PetscViewerDestroy(&viewer);
338:
339: return(0);
340: }