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: The same as ex22.c except it does not use DMMG, it uses its replacement.
14: See src/snes/examples/tutorials/ex50.c
16: Can also be run with -pc_type exotic -ksp_type fgmres
18: */
20: static char help[] = "Solves 3D Laplacian using multigrid.\n\n";
22: #include <petscksp.h>
23: #include <petscdmda.h>
31: int main(int argc,char **argv)
32: {
34: KSP ksp;
35: PetscReal norm;
36: DM da;
37: Vec x,b,r;
38: Mat A;
40: PetscInitialize(&argc,&argv,(char *)0,help);
42: KSPCreate(PETSC_COMM_WORLD,&ksp);
43: DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-7,-7,-7,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,0,&da);
44: DMSetInitialGuess(da,ComputeInitialGuess);
45: DMSetFunction(da,ComputeRHS);
46: DMSetJacobian(da,ComputeMatrix);
47: KSPSetDM(ksp,da);
48: /* KSPSetDMActive(ksp,PETSC_FALSE);*/
49: DMDestroy(&da);
51: KSPSetFromOptions(ksp);
52: KSPSetUp(ksp);
53: KSPSolve(ksp,PETSC_NULL,PETSC_NULL);
54: KSPGetSolution(ksp,&x);
55: KSPGetRhs(ksp,&b);
56: VecDuplicate(b,&r);
57: KSPGetOperators(ksp,&A,PETSC_NULL,PETSC_NULL);
59: MatMult(A,x,r);
60: VecAXPY(r,-1.0,b);
61: VecNorm(r,NORM_2,&norm);
62: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm);
64: VecDestroy(&r);
65: KSPDestroy(&ksp);
66: PetscFinalize();
68: return 0;
69: }
73: PetscErrorCode ComputeRHS(DM dm,Vec x,Vec b)
74: {
76: PetscInt mx,my,mz;
77: PetscScalar h;
78: PetscScalar ***xx,***bb;
79: PetscInt i,j,k,xm,ym,zm,xs,ys,zs;
80: PetscScalar Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
83: DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
84: h = 1.0/((mx-1)*(my-1)*(mz-1));
85: VecSet(b,h);
87: if (x) {
88: DMDAVecGetArray(dm,x,&xx);
89: DMDAVecGetArray(dm,b,&bb);
91: DMDAGetInfo(dm,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
92: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
93: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
94: DMDAGetCorners(dm,&xs,&ys,&zs,&xm,&ym,&zm);
95:
96: for (k=zs; k<zs+zm; k++){
97: for (j=ys; j<ys+ym; j++){
98: for(i=xs; i<xs+xm; i++){
99: if (!(i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1)){
100: bb[k][j][i] -= +HxHydHz*xx[k-1][j][i] + HxHzdHy*xx[k][j-1][i] + HyHzdHx*xx[k][j][i-1] + HyHzdHx*xx[k][j][i+1] + HxHzdHy*xx[k][j+1][i] + HxHydHz*xx[k+1][j][i];
101: }
102: bb[k][j][i] += 2.0*(HxHydHz + HxHzdHy + HyHzdHx)*xx[k][j][i];
103: }
104: }
105: }
106: DMDAVecRestoreArray(dm,x,&xx);
107: DMDAVecRestoreArray(dm,b,&bb);
108: }
109: return(0);
110: }
111:
114: PetscErrorCode ComputeInitialGuess(DM dm,Vec b)
115: {
119: VecSet(b,0);
120: return(0);
121: }
125: PetscErrorCode ComputeMatrix(DM dm,Vec x,Mat jac,Mat B,MatStructure *stflg)
126: {
127: DM da = dm;
129: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs;
130: PetscScalar v[7],Hx,Hy,Hz,HxHydHz,HyHzdHx,HxHzdHy;
131: MatStencil row,col[7];
133: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
134: Hx = 1.0 / (PetscReal)(mx-1); Hy = 1.0 / (PetscReal)(my-1); Hz = 1.0 / (PetscReal)(mz-1);
135: HxHydHz = Hx*Hy/Hz; HxHzdHy = Hx*Hz/Hy; HyHzdHx = Hy*Hz/Hx;
136: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
137:
138: for (k=zs; k<zs+zm; k++){
139: for (j=ys; j<ys+ym; j++){
140: for(i=xs; i<xs+xm; i++){
141: row.i = i; row.j = j; row.k = k;
142: if (i==0 || j==0 || k==0 || i==mx-1 || j==my-1 || k==mz-1){
143: v[0] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);
144: MatSetValuesStencil(B,1,&row,1,&row,v,INSERT_VALUES);
145: } else {
146: v[0] = -HxHydHz;col[0].i = i; col[0].j = j; col[0].k = k-1;
147: v[1] = -HxHzdHy;col[1].i = i; col[1].j = j-1; col[1].k = k;
148: v[2] = -HyHzdHx;col[2].i = i-1; col[2].j = j; col[2].k = k;
149: v[3] = 2.0*(HxHydHz + HxHzdHy + HyHzdHx);col[3].i = row.i; col[3].j = row.j; col[3].k = row.k;
150: v[4] = -HyHzdHx;col[4].i = i+1; col[4].j = j; col[4].k = k;
151: v[5] = -HxHzdHy;col[5].i = i; col[5].j = j+1; col[5].k = k;
152: v[6] = -HxHydHz;col[6].i = i; col[6].j = j; col[6].k = k+1;
153: MatSetValuesStencil(B,1,&row,7,col,v,INSERT_VALUES);
154: }
155: }
156: }
157: }
158: MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
159: MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
160: *stflg = SAME_NONZERO_PATTERN;
161: return 0;
162: }