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;
45: PetscInitialize(&argc, &argv, (char *)0, help);
46: KSPCreate(PETSC_COMM_WORLD, &ksp);
47: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 11, 11, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
48: DMSetFromOptions(da);
49: DMSetUp(da);
50: KSPSetDM(ksp, (DM)da);
51: DMSetApplicationContext(da, &user);
53: user.uu = 1.0;
54: user.tt = 1.0;
56: KSPSetComputeRHS(ksp, ComputeRHS, &user);
57: KSPSetComputeOperators(ksp, ComputeJacobian, &user);
58: KSPSetFromOptions(ksp);
59: KSPSolve(ksp, NULL, NULL);
61: DMDestroy(&da);
62: KSPDestroy(&ksp);
63: PetscFinalize();
64: return 0;
65: }
67: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
68: {
69: UserContext *user = (UserContext *)ctx;
70: PetscInt i, j, M, N, xm, ym, xs, ys;
71: PetscScalar Hx, Hy, pi, uu, tt;
72: PetscScalar **array;
73: DM da;
74: MatNullSpace nullspace;
77: KSPGetDM(ksp, &da);
78: DMDAGetInfo(da, 0, &M, &N, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
79: uu = user->uu;
80: tt = user->tt;
81: pi = 4 * atan(1.0);
82: Hx = 1.0 / (PetscReal)(M);
83: Hy = 1.0 / (PetscReal)(N);
85: DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0); /* Fine grid */
86: DMDAVecGetArray(da, b, &array);
87: for (j = ys; j < ys + ym; j++) {
88: for (i = xs; i < xs + xm; i++) array[j][i] = -PetscCosScalar(uu * pi * ((PetscReal)i + 0.5) * Hx) * PetscCosScalar(tt * pi * ((PetscReal)j + 0.5) * Hy) * Hx * Hy;
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;
121: row.j = j;
123: if (i == 0 || j == 0 || i == M - 1 || j == N - 1) {
124: num = 0;
125: numi = 0;
126: numj = 0;
127: if (j != 0) {
128: v[num] = -HxdHy;
129: col[num].i = i;
130: col[num].j = j - 1;
131: num++;
132: numj++;
133: }
134: if (i != 0) {
135: v[num] = -HydHx;
136: col[num].i = i - 1;
137: col[num].j = j;
138: num++;
139: numi++;
140: }
141: if (i != M - 1) {
142: v[num] = -HydHx;
143: col[num].i = i + 1;
144: col[num].j = j;
145: num++;
146: numi++;
147: }
148: if (j != N - 1) {
149: v[num] = -HxdHy;
150: col[num].i = i;
151: col[num].j = j + 1;
152: num++;
153: numj++;
154: }
155: v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx);
156: col[num].i = i;
157: col[num].j = j;
158: num++;
159: MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES);
160: } else {
161: v[0] = -HxdHy;
162: col[0].i = i;
163: col[0].j = j - 1;
164: v[1] = -HydHx;
165: col[1].i = i - 1;
166: col[1].j = j;
167: v[2] = 2.0 * (HxdHy + HydHx);
168: col[2].i = i;
169: col[2].j = j;
170: v[3] = -HydHx;
171: col[3].i = i + 1;
172: col[3].j = j;
173: v[4] = -HxdHy;
174: col[4].i = i;
175: col[4].j = j + 1;
176: MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES);
177: }
178: }
179: }
180: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
181: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
183: MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
184: MatSetNullSpace(J, nullspace);
185: MatNullSpaceDestroy(&nullspace);
186: return 0;
187: }
189: /*TEST
191: build:
192: requires: !complex !single
194: test:
195: args: -pc_type mg -pc_mg_type full -ksp_type cg -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type svd -ksp_view
197: test:
198: suffix: 2
199: nsize: 4
200: 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
202: test:
203: suffix: 3
204: nsize: 2
205: 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
207: test:
208: suffix: tut_1
209: nsize: 1
210: args: -da_grid_x 4 -da_grid_y 4 -mat_view
212: test:
213: suffix: tut_2
214: requires: superlu_dist parmetis
215: nsize: 4
216: args: -da_grid_x 120 -da_grid_y 120 -pc_type lu -pc_factor_mat_solver_type superlu_dist -ksp_monitor -ksp_view
218: test:
219: suffix: tut_3
220: nsize: 4
221: args: -da_grid_x 1025 -da_grid_y 1025 -pc_type mg -pc_mg_levels 9 -ksp_monitor
223: TEST*/