Actual source code: ex19.c
1: /* "$Id: ex19.c,v 1.12 2001/08/07 21:30:50 bsmith Exp $" */
3: static char help[] ="Solvers Laplacian with multigrid, bad way.\n\
4: -mx <xg>, where <xg> = number of grid points in the x-direction\n\
5: -my <yg>, where <yg> = number of grid points in the y-direction\n\
6: -Nx <npx>, where <npx> = number of processors in the x-direction\n\
7: -Ny <npy>, where <npy> = number of processors in the y-direction\n\n";
9: /*
10: This problem is modeled by
11: the partial differential equation
12:
13: -Laplacian u = g, 0 < x,y < 1,
14:
15: with boundary conditions
16:
17: u = 0 for x = 0, x = 1, y = 0, y = 1.
18:
19: A finite difference approximation with the usual 5-point stencil
20: is used to discretize the boundary value problem to obtain a nonlinear
21: system of equations.
22: */
24: #include petscsles.h
25: #include petscda.h
26: #include petscmg.h
28: /* User-defined application contexts */
30: typedef struct {
31: int mx,my; /* number grid points in x and y direction */
32: Vec localX,localF; /* local vectors with ghost region */
33: DA da;
34: Vec x,b,r; /* global vectors */
35: Mat J; /* Jacobian on grid */
36: } GridCtx;
38: typedef struct {
39: GridCtx fine;
40: GridCtx coarse;
41: SLES sles_coarse;
42: int ratio;
43: Mat I; /* interpolation from coarse to fine */
44: } AppCtx;
46: #define COARSE_LEVEL 0
47: #define FINE_LEVEL 1
49: extern int FormJacobian_Grid(AppCtx *,GridCtx *,Mat *);
51: /*
52: Mm_ratio - ration of grid lines between fine and coarse grids.
53: */
56: int main(int argc,char **argv)
57: {
58: AppCtx user;
59: int ierr,its,N,n,Nx = PETSC_DECIDE,Ny = PETSC_DECIDE;
60: int size,nlocal,Nlocal;
61: SLES sles,sles_fine;
62: PC pc;
63: PetscScalar one = 1.0;
65: PetscInitialize(&argc,&argv,PETSC_NULL,help);
67: user.ratio = 2;
68: user.coarse.mx = 5; user.coarse.my = 5;
69: PetscOptionsGetInt(PETSC_NULL,"-Mx",&user.coarse.mx,PETSC_NULL);
70: PetscOptionsGetInt(PETSC_NULL,"-My",&user.coarse.my,PETSC_NULL);
71: PetscOptionsGetInt(PETSC_NULL,"-ratio",&user.ratio,PETSC_NULL);
72: user.fine.mx = user.ratio*(user.coarse.mx-1)+1; user.fine.my = user.ratio*(user.coarse.my-1)+1;
74: PetscPrintf(PETSC_COMM_WORLD,"Coarse grid size %d by %d\n",user.coarse.mx,user.coarse.my);
75: PetscPrintf(PETSC_COMM_WORLD,"Fine grid size %d by %d\n",user.fine.mx,user.fine.my);
77: n = user.fine.mx*user.fine.my; N = user.coarse.mx*user.coarse.my;
79: MPI_Comm_size(PETSC_COMM_WORLD,&size);
80: PetscOptionsGetInt(PETSC_NULL,"-Nx",&Nx,PETSC_NULL);
81: PetscOptionsGetInt(PETSC_NULL,"-Ny",&Ny,PETSC_NULL);
83: /* Set up distributed array for fine grid */
84: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.fine.mx,
85: user.fine.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.fine.da);
86: DACreateGlobalVector(user.fine.da,&user.fine.x);
87: VecDuplicate(user.fine.x,&user.fine.r);
88: VecDuplicate(user.fine.x,&user.fine.b);
89: VecGetLocalSize(user.fine.x,&nlocal);
90: DACreateLocalVector(user.fine.da,&user.fine.localX);
91: VecDuplicate(user.fine.localX,&user.fine.localF);
92: MatCreateMPIAIJ(PETSC_COMM_WORLD,nlocal,nlocal,n,n,5,PETSC_NULL,3,PETSC_NULL,&user.fine.J);
94: /* Set up distributed array for coarse grid */
95: DACreate2d(PETSC_COMM_WORLD,DA_NONPERIODIC,DA_STENCIL_STAR,user.coarse.mx,
96: user.coarse.my,Nx,Ny,1,1,PETSC_NULL,PETSC_NULL,&user.coarse.da);
97: DACreateGlobalVector(user.coarse.da,&user.coarse.x);
98: VecDuplicate(user.coarse.x,&user.coarse.b);
99: VecGetLocalSize(user.coarse.x,&Nlocal);
100: DACreateLocalVector(user.coarse.da,&user.coarse.localX);
101: VecDuplicate(user.coarse.localX,&user.coarse.localF);
102: MatCreateMPIAIJ(PETSC_COMM_WORLD,Nlocal,Nlocal,N,N,5,PETSC_NULL,3,PETSC_NULL,&user.coarse.J);
104: /* Create linear solver */
105: SLESCreate(PETSC_COMM_WORLD,&sles);
107: /* set two level additive Schwarz preconditioner */
108: SLESGetPC(sles,&pc);
109: PCSetType(pc,PCMG);
110: MGSetLevels(pc,2,PETSC_NULL);
111: MGSetType(pc,MGADDITIVE);
113: FormJacobian_Grid(&user,&user.coarse,&user.coarse.J);
114: FormJacobian_Grid(&user,&user.fine,&user.fine.J);
116: /* Create coarse level */
117: MGGetCoarseSolve(pc,&user.sles_coarse);
118: SLESSetOptionsPrefix(user.sles_coarse,"coarse_");
119: SLESSetFromOptions(user.sles_coarse);
120: SLESSetOperators(user.sles_coarse,user.coarse.J,user.coarse.J,DIFFERENT_NONZERO_PATTERN);
121: MGSetX(pc,COARSE_LEVEL,user.coarse.x);
122: MGSetRhs(pc,COARSE_LEVEL,user.coarse.b);
124: /* Create fine level */
125: MGGetSmoother(pc,FINE_LEVEL,&sles_fine);
126: SLESSetOptionsPrefix(sles_fine,"fine_");
127: SLESSetFromOptions(sles_fine);
128: SLESSetOperators(sles_fine,user.fine.J,user.fine.J,DIFFERENT_NONZERO_PATTERN);
129: MGSetR(pc,FINE_LEVEL,user.fine.r);
130: MGSetResidual(pc,FINE_LEVEL,MGDefaultResidual,user.fine.J);
132: /* Create interpolation between the levels */
133: DAGetInterpolation(user.coarse.da,user.fine.da,&user.I,PETSC_NULL);
134: MGSetInterpolate(pc,FINE_LEVEL,user.I);
135: MGSetRestriction(pc,FINE_LEVEL,user.I);
137: SLESSetOperators(sles,user.fine.J,user.fine.J,DIFFERENT_NONZERO_PATTERN);
139: VecSet(&one,user.fine.b);
140: {
141: PetscRandom rand;
142: PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rand);
143: VecSetRandom(rand,user.fine.b);
144: PetscRandomDestroy(rand);
145: }
147: /* Set options, then solve nonlinear system */
148: SLESSetFromOptions(sles);
150: SLESSolve(sles,user.fine.b,user.fine.x,&its);
151: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %d\n",its);
153: /* Free data structures */
154: MatDestroy(user.fine.J);
155: VecDestroy(user.fine.x);
156: VecDestroy(user.fine.r);
157: VecDestroy(user.fine.b);
158: DADestroy(user.fine.da);
159: VecDestroy(user.fine.localX);
160: VecDestroy(user.fine.localF);
162: MatDestroy(user.coarse.J);
163: VecDestroy(user.coarse.x);
164: VecDestroy(user.coarse.b);
165: DADestroy(user.coarse.da);
166: VecDestroy(user.coarse.localX);
167: VecDestroy(user.coarse.localF);
169: SLESDestroy(sles);
170: MatDestroy(user.I);
171: PetscFinalize();
173: return 0;
174: }
178: int FormJacobian_Grid(AppCtx *user,GridCtx *grid,Mat *J)
179: {
180: Mat jac = *J;
181: int ierr,i,j,row,mx,my,xs,ys,xm,ym,Xs,Ys,Xm,Ym,col[5];
182: int nloc,*ltog,grow;
183: PetscScalar two = 2.0,one = 1.0,v[5],hx,hy,hxdhy,hydhx,value;
185: mx = grid->mx; my = grid->my;
186: hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1);
187: hxdhy = hx/hy; hydhx = hy/hx;
189: /* Get ghost points */
190: DAGetCorners(grid->da,&xs,&ys,0,&xm,&ym,0);
191: DAGetGhostCorners(grid->da,&Xs,&Ys,0,&Xm,&Ym,0);
192: DAGetGlobalIndices(grid->da,&nloc,<og);
194: /* Evaluate Jacobian of function */
195: for (j=ys; j<ys+ym; j++) {
196: row = (j - Ys)*Xm + xs - Xs - 1;
197: for (i=xs; i<xs+xm; i++) {
198: row++;
199: grow = ltog[row];
200: if (i > 0 && i < mx-1 && j > 0 && j < my-1) {
201: v[0] = -hxdhy; col[0] = ltog[row - Xm];
202: v[1] = -hydhx; col[1] = ltog[row - 1];
203: v[2] = two*(hydhx + hxdhy); col[2] = grow;
204: v[3] = -hydhx; col[3] = ltog[row + 1];
205: v[4] = -hxdhy; col[4] = ltog[row + Xm];
206: MatSetValues(jac,1,&grow,5,col,v,INSERT_VALUES);
207: } else if ((i > 0 && i < mx-1) || (j > 0 && j < my-1)){
208: value = .5*two*(hydhx + hxdhy);
209: MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);
210: } else {
211: value = .25*two*(hydhx + hxdhy);
212: MatSetValues(jac,1,&grow,1,&grow,&value,INSERT_VALUES);
213: }
214: }
215: }
216: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
217: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
219: return 0;
220: }