Actual source code: ex45.c
2: /*
3: Laplacian in 3D. Modeled by the partial differential equation
5: - Laplacian u = 1,0 < x,y,z < 1,
7: with boundary conditions
9: u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
11: This uses multigrid to solve the linear system
13: See src/snes/tutorials/ex50.c
15: Can also be run with -pc_type exotic -ksp_type fgmres
17: */
19: static char help[] = "Solves 3D Laplacian using multigrid.\n\n";
21: #include <petscksp.h>
22: #include <petscdm.h>
23: #include <petscdmda.h>
25: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
26: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
27: extern PetscErrorCode ComputeInitialGuess(KSP,Vec,void*);
29: int main(int argc,char **argv)
30: {
31: KSP ksp;
32: PetscReal norm;
33: DM da;
34: Vec x,b,r;
35: Mat A;
37: PetscInitialize(&argc,&argv,(char*)0,help);
39: KSPCreate(PETSC_COMM_WORLD,&ksp);
40: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,7,7,7,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
41: DMSetFromOptions(da);
42: DMSetUp(da);
43: KSPSetDM(ksp,da);
44: KSPSetComputeInitialGuess(ksp,ComputeInitialGuess,NULL);
45: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
46: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
47: DMDestroy(&da);
49: KSPSetFromOptions(ksp);
50: KSPSolve(ksp,NULL,NULL);
51: KSPGetSolution(ksp,&x);
52: KSPGetRhs(ksp,&b);
53: VecDuplicate(b,&r);
54: KSPGetOperators(ksp,&A,NULL);
56: MatMult(A,x,r);
57: VecAXPY(r,-1.0,b);
58: VecNorm(r,NORM_2,&norm);
59: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
61: VecDestroy(&r);
62: KSPDestroy(&ksp);
63: PetscFinalize();
64: return 0;
65: }
67: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
68: {
69: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
70: DM dm;
71: PetscScalar Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
72: PetscScalar ***barray;
75: KSPGetDM(ksp,&dm);
76: DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
77: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
78: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
79: DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);
80: DMDAVecGetArray(dm,b,&barray);
82: for (k=zs; k<zs+zm; k++) {
83: for (j=ys; j<ys+ym; j++) {
84: for (i=xs; i<xs+xm; i++) {
85: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
86: barray[k][j][i] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
87: } else {
88: barray[k][j][i] = Hx*Hy*Hz;
89: }
90: }
91: }
92: }
93: DMDAVecRestoreArray(dm,b,&barray);
94: return 0;
95: }
97: PetscErrorCode ComputeInitialGuess(KSP ksp,Vec b,void *ctx)
98: {
100: VecSet(b,0);
101: return 0;
102: }
104: PetscErrorCode ComputeMatrix(KSP ksp,Mat jac,Mat B,void *ctx)
105: {
106: DM da;
107: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
108: PetscScalar v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
109: MatStencil row,col[7];
112: KSPGetDM(ksp,&da);
113: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
114: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
115: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
116: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
118: for (k=zs; k<zs+zm; k++) {
119: for (j=ys; j<ys+ym; j++) {
120: for (i=xs; i<xs+xm; i++) {
121: row.i = i; row.j = j; row.k = k;
122: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
123: v[0] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
124: MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
125: } else {
126: v[0] = -HxHydHz;col[0].i = i; col[0].j = j; col[0].k = k-1;
127: v[1] = -HxHzdHy;col[1].i = i; col[1].j = j-1; col[1].k = k;
128: v[2] = -HyHzdHx;col[2].i = i-1; col[2].j = j; col[2].k = k;
129: v[3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);col[3].i = row.i; col[3].j = row.j; col[3].k = row.k;
130: v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
131: v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
132: v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
133: MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);
134: }
135: }
136: }
137: }
138: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
139: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
140: return 0;
141: }
143: /*TEST
145: test:
146: nsize: 4
147: args: -pc_type exotic -ksp_monitor_short -ksp_type fgmres -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 1 -mg_levels_pc_type bjacobi
148: output_file: output/ex45_1.out
150: test:
151: suffix: 2
152: nsize: 4
153: args: -ksp_monitor_short -da_grid_x 21 -da_grid_y 21 -da_grid_z 21 -pc_type mg -pc_mg_levels 3 -mg_levels_ksp_type richardson -mg_levels_ksp_max_it 1 -mg_levels_pc_type bjacobi
155: test:
156: suffix: telescope
157: nsize: 4
158: args: -ksp_type fgmres -ksp_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -pc_mg_levels 2 -da_grid_x 65 -da_grid_y 65 -da_grid_z 65 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_ignore_kspcomputeoperators -mg_coarse_pc_telescope_reduction_factor 4 -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_galerkin pmat -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_levels_ksp_type richardson -mg_coarse_telescope_mg_levels_pc_type jacobi -mg_levels_ksp_type richardson -mg_coarse_telescope_mg_levels_ksp_type richardson -ksp_rtol 1.0e-4
160: test:
161: suffix: telescope_2
162: nsize: 4
163: args: -ksp_type fgmres -ksp_monitor_short -pc_type mg -mg_levels_ksp_type richardson -mg_levels_pc_type jacobi -pc_mg_levels 2 -da_grid_x 65 -da_grid_y 65 -da_grid_z 65 -mg_coarse_pc_type telescope -mg_coarse_pc_telescope_reduction_factor 2 -mg_coarse_telescope_pc_type mg -mg_coarse_telescope_pc_mg_galerkin pmat -mg_coarse_telescope_pc_mg_levels 3 -mg_coarse_telescope_mg_levels_ksp_type richardson -mg_coarse_telescope_mg_levels_pc_type jacobi -mg_levels_ksp_type richardson -mg_coarse_telescope_mg_levels_ksp_type richardson -ksp_rtol 1.0e-4
165: TEST*/