Actual source code: ex34.c
petsc-3.4.5 2014-06-29
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 <petscdmda.h>
26: #include <petscksp.h>
28: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,MatStructure*,void*);
29: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
31: typedef enum {DIRICHLET, NEUMANN} BCType;
33: typedef struct {
34: BCType bcType;
35: } UserContext;
39: int main(int argc,char **argv)
40: {
41: KSP ksp;
42: DM da;
43: UserContext user;
44: PetscReal norm;
45: const char *bcTypes[2] = {"dirichlet","neumann"};
47: PetscInt bc;
49: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
50: PetscScalar Hx,Hy,Hz;
51: PetscScalar ***array;
52: Vec x,b,r;
53: Mat J;
55: PetscInitialize(&argc,&argv,(char*)0,help);
57: KSPCreate(PETSC_COMM_WORLD,&ksp);
58: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,12,12,12,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
59: DMDASetInterpolationType(da, DMDA_Q0);
61: KSPSetDM(ksp,da);
63: PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
64: bc = (PetscInt)NEUMANN;
65: PetscOptionsEList("-bc_type","Type of boundary condition","ex34.c",bcTypes,2,bcTypes[0],&bc,NULL);
66: user.bcType = (BCType)bc;
67: PetscOptionsEnd();
69: KSPSetComputeRHS(ksp,ComputeRHS,&user);
70: KSPSetComputeOperators(ksp,ComputeMatrix,&user);
71: KSPSetFromOptions(ksp);
72: KSPSolve(ksp,NULL,NULL);
73: KSPGetSolution(ksp,&x);
74: KSPGetRhs(ksp,&b);
75: KSPGetOperators(ksp,NULL,&J,NULL);
76: VecDuplicate(b,&r);
78: MatMult(J,x,r);
79: VecAXPY(r,-1.0,b);
80: VecNorm(r,NORM_2,&norm);
81: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
83: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
84: Hx = 1.0 / (PetscReal)(mx);
85: Hy = 1.0 / (PetscReal)(my);
86: Hz = 1.0 / (PetscReal)(mz);
87: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
88: DMDAVecGetArray(da, x, &array);
90: for (k=zs; k<zs+zm; k++) {
91: for (j=ys; j<ys+ym; j++) {
92: for (i=xs; i<xs+xm; i++) {
93: array[k][j][i] -=
94: PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))*
95: PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))*
96: PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz));
97: }
98: }
99: }
100: DMDAVecRestoreArray(da, x, &array);
101: VecAssemblyBegin(x);
102: VecAssemblyEnd(x);
104: VecNorm(x,NORM_INFINITY,&norm);
105: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm);
106: VecNorm(x,NORM_1,&norm);
107: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
108: VecNorm(x,NORM_2,&norm);
109: PetscPrintf(PETSC_COMM_WORLD,"Error norm %g\n",norm/((PetscReal)(mx)*(PetscReal)(my)*(PetscReal)(mz)));
111: VecDestroy(&r);
112: KSPDestroy(&ksp);
113: DMDestroy(&da);
114: PetscFinalize();
115: return 0;
116: }
120: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
121: {
122: UserContext *user = (UserContext*)ctx;
124: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
125: PetscScalar Hx,Hy,Hz;
126: PetscScalar ***array;
127: DM da;
130: KSPGetDM(ksp,&da);
131: DMDAGetInfo(da, 0, &mx, &my, &mz, 0,0,0,0,0,0,0,0,0);
132: Hx = 1.0 / (PetscReal)(mx);
133: Hy = 1.0 / (PetscReal)(my);
134: Hz = 1.0 / (PetscReal)(mz);
135: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
136: DMDAVecGetArray(da, b, &array);
137: for (k=zs; k<zs+zm; k++) {
138: for (j=ys; j<ys+ym; j++) {
139: for (i=xs; i<xs+xm; i++) {
140: array[k][j][i] = 12 * PETSC_PI * PETSC_PI
141: * PetscCosScalar(2*PETSC_PI*(((PetscReal)i+0.5)*Hx))
142: * PetscCosScalar(2*PETSC_PI*(((PetscReal)j+0.5)*Hy))
143: * PetscCosScalar(2*PETSC_PI*(((PetscReal)k+0.5)*Hz))
144: * Hx * Hy * Hz;
145: }
146: }
147: }
148: DMDAVecRestoreArray(da, b, &array);
149: VecAssemblyBegin(b);
150: VecAssemblyEnd(b);
152: /* force right hand side to be consistent for singular matrix */
153: /* note this is really a hack, normally the model would provide you with a consistent right handside */
154: if (user->bcType == NEUMANN) {
155: MatNullSpace nullspace;
157: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
158: MatNullSpaceRemove(nullspace,b,NULL);
159: MatNullSpaceDestroy(&nullspace);
160: }
161: return(0);
162: }
167: PetscErrorCode ComputeMatrix(KSP ksp, Mat J,Mat jac,MatStructure *str, void *ctx)
168: {
169: UserContext *user = (UserContext*)ctx;
171: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,num, numi, numj, numk;
172: PetscScalar v[7],Hx,Hy,Hz,HyHzdHx,HxHzdHy,HxHydHz;
173: MatStencil row, col[7];
174: DM da;
177: KSPGetDM(ksp,&da);
178: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
179: Hx = 1.0 / (PetscReal)(mx);
180: Hy = 1.0 / (PetscReal)(my);
181: Hz = 1.0 / (PetscReal)(mz);
182: HyHzdHx = Hy*Hz/Hx;
183: HxHzdHy = Hx*Hz/Hy;
184: HxHydHz = Hx*Hy/Hz;
185: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
186: for (k=zs; k<zs+zm; k++) {
187: for (j=ys; j<ys+ym; j++) {
188: for (i=xs; i<xs+xm; i++) {
189: row.i = i; row.j = j; row.k = k;
190: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1) {
191: if (user->bcType == DIRICHLET) {
192: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
193: v[0] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz);
194: MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
195: } else if (user->bcType == NEUMANN) {
196: num = 0; numi=0; numj=0; numk=0;
197: if (k!=0) {
198: v[num] = -HxHydHz;
199: col[num].i = i;
200: col[num].j = j;
201: col[num].k = k-1;
202: num++; numk++;
203: }
204: if (j!=0) {
205: v[num] = -HxHzdHy;
206: col[num].i = i;
207: col[num].j = j-1;
208: col[num].k = k;
209: num++; numj++;
210: }
211: if (i!=0) {
212: v[num] = -HyHzdHx;
213: col[num].i = i-1;
214: col[num].j = j;
215: col[num].k = k;
216: num++; numi++;
217: }
218: if (i!=mx-1) {
219: v[num] = -HyHzdHx;
220: col[num].i = i+1;
221: col[num].j = j;
222: col[num].k = k;
223: num++; numi++;
224: }
225: if (j!=my-1) {
226: v[num] = -HxHzdHy;
227: col[num].i = i;
228: col[num].j = j+1;
229: col[num].k = k;
230: num++; numj++;
231: }
232: if (k!=mz-1) {
233: v[num] = -HxHydHz;
234: col[num].i = i;
235: col[num].j = j;
236: col[num].k = k+1;
237: num++; numk++;
238: }
239: v[num] = (PetscReal)(numk)*HxHydHz + (PetscReal)(numj)*HxHzdHy + (PetscReal)(numi)*HyHzdHx;
240: col[num].i = i; col[num].j = j; col[num].k = k;
241: num++;
242: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
243: }
244: } else {
245: v[0] = -HxHydHz; col[0].i = i; col[0].j = j; col[0].k = k-1;
246: v[1] = -HxHzdHy; col[1].i = i; col[1].j = j-1; col[1].k = k;
247: v[2] = -HyHzdHx; col[2].i = i-1; col[2].j = j; col[2].k = k;
248: v[3] = 2.0*(HyHzdHx + HxHzdHy + HxHydHz); col[3].i = i; col[3].j = j; col[3].k = k;
249: v[4] = -HyHzdHx; col[4].i = i+1; col[4].j = j; col[4].k = k;
250: v[5] = -HxHzdHy; col[5].i = i; col[5].j = j+1; col[5].k = k;
251: v[6] = -HxHydHz; col[6].i = i; col[6].j = j; col[6].k = k+1;
252: MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);
253: }
254: }
255: }
256: }
257: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
258: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
259: if (user->bcType == NEUMANN) {
260: MatNullSpace nullspace;
262: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
263: MatSetNullSpace(jac,nullspace);
264: MatNullSpaceDestroy(&nullspace);
265: }
266: return(0);
267: }