Actual source code: ex34.c
petsc-3.8.4 2018-03-24
1: /*T
2: Concepts: KSP^solving a system of linear equations
3: Concepts: KSP^Laplacian, 3d
4: Processors: n
5: T*/
7: /*
8: Laplacian in 3D. Modeled by the partial differential equation
10: div grad u = f, 0 < x,y,z < 1,
12: with pure Neumann boundary conditions
14: u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1.
16: The functions are cell-centered
18: This uses multigrid to solve the linear system
20: Contributed by Jianming Yang <jianming-yang@uiowa.edu>
21: */
23: static char help[] = "Solves 3D Laplacian using multigrid.\n\n";
25: #include <petscdm.h>
26: #include <petscdmda.h>
27: #include <petscksp.h>
29: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
30: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
32: int main(int argc,char **argv)
33: {
34: KSP ksp;
35: DM da;
36: PetscReal norm;
38: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
39: PetscScalar Hx,Hy,Hz;
40: PetscScalar ***array;
41: Vec x,b,r;
42: Mat J;
44: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
45: KSPCreate(PETSC_COMM_WORLD,&ksp);
46: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,12,12,12,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
47: DMSetFromOptions(da);
48: DMSetUp(da);
49: DMDASetInterpolationType(da, DMDA_Q0);
51: KSPSetDM(ksp,da);
53: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
54: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
55: KSPSetFromOptions(ksp);
56: KSPSolve(ksp,NULL,NULL);
57: KSPGetSolution(ksp,&x);
58: KSPGetRhs(ksp,&b);
59: KSPGetOperators(ksp,NULL,&J);
60: VecDuplicate(b,&r);
62: MatMult(J,x,r);
63: VecAXPY(r,-1.0,b);
64: VecNorm(r,NORM_2,&norm);
65: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
67: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
68: Hx = 1.0 / (PetscReal)(mx);
69: Hy = 1.0 / (PetscReal)(my);
70: Hz = 1.0 / (PetscReal)(mz);
71: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
72: DMDAVecGetArray(da, x, &array);
74: for (k=zs; k<zs+zm; k++) {
75: for (j=ys; j<ys+ym; j++) {
76: for (i=xs; i<xs+xm; i++) {
77: array[k][j][i] -=
78: PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
79: PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
80: PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
81: }
82: }
83: }
84: DMDAVecRestoreArray(da, x, &array);
85: VecAssemblyBegin(x);
86: VecAssemblyEnd(x);
88: VecNorm(x,NORM_INFINITY,&norm);
89: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)norm);
90: VecNorm(x,NORM_1,&norm);
91: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
92: VecNorm(x,NORM_2,&norm);
93: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",(double)(norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz))));
95: VecDestroy(&r);
96: KSPDestroy(&ksp);
97: DMDestroy(&da);
98: PetscFinalize();
99: return ierr;
100: }
102: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
103: {
105: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
106: PetscScalar Hx,Hy,Hz;
107: PetscScalar ***array;
108: DM da;
109: MatNullSpace nullspace;
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);
115: Hy = 1.0 / (PetscReal)(my);
116: Hz = 1.0 / (PetscReal)(mz);
117: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
118: DMDAVecGetArray(da, b, &array);
119: for (k=zs; k<zs+zm; k++) {
120: for (j=ys; j<ys+ym; j++) {
121: for (i=xs; i<xs+xm; i++) {
122: array[k][j][i] = 12 * PETSC_PI * PETSC_PI
123: * PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
124: * PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
125: * PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
126: * Hx * Hy * Hz;
127: }
128: }
129: }
130: DMDAVecRestoreArray(da, b, &array);
131: VecAssemblyBegin(b);
132: VecAssemblyEnd(b);
134: /* force right hand side to be consistent for singular matrix */
135: /* note this is really a hack, normally the model would provide you with a consistent right handside */
137: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
138: MatNullSpaceRemove(nullspace,b);
139: MatNullSpaceDestroy(&nullspace);
140: return(0);
141: }
144: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac, void *ctx)
145: {
147: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
148: PetscScalar v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
149: MatStencil row, col[7];
150: DM da;
151: MatNullSpace nullspace;
154: KSPGetDM(ksp,&da);
155: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
156: Hx = 1.0 / (PetscReal)(mx);
157: Hy = 1.0 / (PetscReal)(my);
158: Hz = 1.0 / (PetscReal)(mz);
159: HyHzdHx = Hy*Hz/Hx;
160: HxHzdHy = Hx*Hz/Hy;
161: HxHydHz = Hx*Hy/Hz;
162: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
163: for (k=zs; k<zs+zm; k++) {
164: for (j=ys; j<ys+ym; j++) {
165: for (i=xs; i<xs+xm; i++) {
166: row.i = i; row.j = j; row.k = k;
167: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
168: num = 0; numi=0; numj=0; numk=0;
169: if (k!=0) {
170: v[num] = -HxHydHz;
171: col[num].i = i;
172: col[num].j = j;
173: col[num].k = k-1;
174: num++; numk++;
175: }
176: if (j!=0) {
177: v[num] = -HxHzdHy;
178: col[num].i = i;
179: col[num].j = j-1;
180: col[num].k = k;
181: num++; numj++;
182: }
183: if (i!=0) {
184: v[num] = -HyHzdHx;
185: col[num].i = i-1;
186: col[num].j = j;
187: col[num].k = k;
188: num++; numi++;
189: }
190: if (i!=mx-1) {
191: v[num] = -HyHzdHx;
192: col[num].i = i+1;
193: col[num].j = j;
194: col[num].k = k;
195: num++; numi++;
196: }
197: if (j!=my-1) {
198: v[num] = -HxHzdHy;
199: col[num].i = i;
200: col[num].j = j+1;
201: col[num].k = k;
202: num++; numj++;
203: }
204: if (k!=mz-1) {
205: v[num] = -HxHydHz;
206: col[num].i = i;
207: col[num].j = j;
208: col[num].k = k+1;
209: num++; numk++;
210: }
211: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
212: col[num].i = i; col[num].j = j; col[num].k = k;
213: num++;
214: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
215: } else {
216: v[0] = -HxHydHz; col[0].i = i; col[0].j = j; col[0].k = k-1;
217: v[1] = -HxHzdHy; col[1].i = i; col[1].j = j-1; col[1].k = k;
218: v[2] = -HyHzdHx; col[2].i = i-1; col[2].j = j; col[2].k = k;
219: v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i; col[3].j = j; col[3].k = k;
220: v[4] = -HyHzdHx; col[4].i = i+1; col[4].j = j; col[4].k = k;
221: v[5] = -HxHzdHy; col[5].i = i; col[5].j = j+1; col[5].k = k;
222: v[6] = -HxHydHz; col[6].i = i; col[6].j = j; col[6].k = k+1;
223: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
224: }
225: }
226: }
227: }
228: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
229: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
230: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
231: MatSetNullSpace(J,nullspace);
232: MatNullSpaceDestroy(&nullspace);
233: return(0);
234: }