Actual source code: ex50.c
1: /* DMDA/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/tutorials/ex29.c and ex32.c
14: Compare to ex66.c
16: Example of Usage:
17: ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 3 -ksp_monitor -ksp_view -dm_view draw -draw_pause -1
18: ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type ilu -mg_levels_0_pc_factor_levels 1 -ksp_monitor -ksp_view
19: ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor
20: mpiexec -n 4 ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 10 -ksp_monitor -ksp_view -log_view
21: */
23: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";
25: #include <petscdm.h>
26: #include <petscdmda.h>
27: #include <petscksp.h>
28: #include <petscsys.h>
29: #include <petscvec.h>
31: extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,void*);
32: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
34: typedef struct {
35: PetscScalar uu, tt;
36: } UserContext;
38: int main(int argc,char **argv)
39: {
40: KSP ksp;
41: DM da;
42: UserContext user;
44: PetscInitialize(&argc,&argv,(char*)0,help);
45: KSPCreate(PETSC_COMM_WORLD,&ksp);
46: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,11,11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
47: DMSetFromOptions(da);
48: DMSetUp(da);
49: KSPSetDM(ksp,(DM)da);
50: DMSetApplicationContext(da,&user);
52: user.uu = 1.0;
53: user.tt = 1.0;
55: KSPSetComputeRHS(ksp,ComputeRHS,&user);
56: KSPSetComputeOperators(ksp,ComputeJacobian,&user);
57: KSPSetFromOptions(ksp);
58: KSPSolve(ksp,NULL,NULL);
60: DMDestroy(&da);
61: KSPDestroy(&ksp);
62: PetscFinalize();
63: return 0;
64: }
66: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
67: {
68: UserContext *user = (UserContext*)ctx;
69: PetscInt i,j,M,N,xm,ym,xs,ys;
70: PetscScalar Hx,Hy,pi,uu,tt;
71: PetscScalar **array;
72: DM da;
73: MatNullSpace nullspace;
76: KSPGetDM(ksp,&da);
77: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
78: uu = user->uu; tt = user->tt;
79: pi = 4*atan(1.0);
80: Hx = 1.0/(PetscReal)(M);
81: Hy = 1.0/(PetscReal)(N);
83: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
84: DMDAVecGetArray(da, b, &array);
85: for (j=ys; j<ys+ym; j++) {
86: for (i=xs; i<xs+xm; i++) {
87: array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*PetscCosScalar(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
88: }
89: }
90: DMDAVecRestoreArray(da, b, &array);
91: VecAssemblyBegin(b);
92: VecAssemblyEnd(b);
94: /* force right hand side to be consistent for singular matrix */
95: /* note this is really a hack, normally the model would provide you with a consistent right handside */
96: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
97: MatNullSpaceRemove(nullspace,b);
98: MatNullSpaceDestroy(&nullspace);
99: return 0;
100: }
102: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
103: {
104: PetscInt i, j, M, N, xm, ym, xs, ys, num, numi, numj;
105: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
106: MatStencil row, col[5];
107: DM da;
108: MatNullSpace nullspace;
111: KSPGetDM(ksp,&da);
112: DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
113: Hx = 1.0 / (PetscReal)(M);
114: Hy = 1.0 / (PetscReal)(N);
115: HxdHy = Hx/Hy;
116: HydHx = Hy/Hx;
117: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
118: for (j=ys; j<ys+ym; j++) {
119: for (i=xs; i<xs+xm; i++) {
120: row.i = i; row.j = j;
122: if (i==0 || j==0 || i==M-1 || j==N-1) {
123: num=0; numi=0; numj=0;
124: if (j!=0) {
125: v[num] = -HxdHy; col[num].i = i; col[num].j = j-1;
126: num++; numj++;
127: }
128: if (i!=0) {
129: v[num] = -HydHx; col[num].i = i-1; col[num].j = j;
130: num++; numi++;
131: }
132: if (i!=M-1) {
133: v[num] = -HydHx; col[num].i = i+1; col[num].j = j;
134: num++; numi++;
135: }
136: if (j!=N-1) {
137: v[num] = -HxdHy; col[num].i = i; col[num].j = j+1;
138: num++; numj++;
139: }
140: v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i; col[num].j = j;
141: num++;
142: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
143: } else {
144: v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
145: v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
146: v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
147: v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
148: v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
149: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
150: }
151: }
152: }
153: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
154: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
156: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
157: MatSetNullSpace(J,nullspace);
158: MatNullSpaceDestroy(&nullspace);
159: return 0;
160: }
162: /*TEST
164: build:
165: requires: !complex !single
167: test:
168: args: -pc_type mg -pc_mg_type full -ksp_type cg -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type svd -ksp_view
170: test:
171: suffix: 2
172: nsize: 4
173: args: -pc_type mg -pc_mg_type full -ksp_type cg -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd -ksp_view
175: test:
176: suffix: 3
177: nsize: 2
178: args: -pc_type mg -pc_mg_type full -ksp_monitor_short -da_refine 5 -mg_coarse_ksp_type cg -mg_coarse_ksp_converged_reason -mg_coarse_ksp_rtol 1e-2 -mg_coarse_ksp_max_it 5 -mg_coarse_pc_type none -pc_mg_levels 2 -ksp_type pipefgmres -ksp_pipefgmres_shift 1.5
180: test:
181: suffix: tut_1
182: nsize: 1
183: args: -da_grid_x 4 -da_grid_y 4 -mat_view
185: test:
186: suffix: tut_2
187: requires: superlu_dist parmetis
188: nsize: 4
189: args: -da_grid_x 120 -da_grid_y 120 -pc_type lu -pc_factor_mat_solver_type superlu_dist -ksp_monitor -ksp_view
191: test:
192: suffix: tut_3
193: nsize: 4
194: args: -da_grid_x 1025 -da_grid_y 1025 -pc_type mg -pc_mg_levels 9 -ksp_monitor
196: TEST*/