Actual source code: ex50.c
petsc-3.4.5 2014-06-29
1: /* DM/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/examples/tutorials/ex29.c and ex32.c
14: Example of Usage:
15: ./ex50 -mglevels 3 -ksp_monitor -M 3 -N 3 -ksp_view -dm_view draw -draw_pause -1
16: ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_factor_levels <ilu_levels> -ksp_monitor -cmp_solu
17: ./ex50 -M 100 -N 100 -mglevels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor -cmp_solu
18: mpiexec -n 4 ./ex50 -M 3 -N 3 -ksp_monitor -ksp_view -mglevels 10 -log_summary
19: */
21: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";
23: #include <petscdmda.h>
24: #include <petscksp.h>
25: #include <petscsys.h>
26: #include <petscvec.h>
28: extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,MatStructure*,void*);
29: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
30: extern PetscErrorCode ComputeTrueSolution(DM, Vec);
31: extern PetscErrorCode VecView_VTK(Vec, const char [], const char []);
33: typedef enum {DIRICHLET, NEUMANN} BCType;
35: typedef struct {
36: PetscScalar uu, tt;
37: BCType bcType;
38: } UserContext;
42: int main(int argc,char **argv)
43: {
44: KSP ksp;
45: DM da;
46: UserContext user;
47: PetscInt bc;
50: PetscInitialize(&argc,&argv,(char*)0,help);
51: KSPCreate(PETSC_COMM_WORLD,&ksp);
52: DMDACreate2d(PETSC_COMM_WORLD, DMDA_BOUNDARY_NONE, DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
53: KSPSetDM(ksp,(DM)da);
54: DMSetApplicationContext(da,&user);
56: user.uu = 1.0;
57: user.tt = 1.0;
58: bc = (PetscInt)NEUMANN; /* Use Neumann Boundary Conditions */
59: user.bcType = (BCType)bc;
62: KSPSetComputeRHS(ksp,ComputeRHS,&user);
63: KSPSetComputeOperators(ksp,ComputeJacobian,&user);
64: KSPSetFromOptions(ksp);
65: KSPSolve(ksp,NULL,NULL);
67: DMDestroy(&da);
68: KSPDestroy(&ksp);
69: PetscFinalize();
70: return 0;
71: }
75: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
76: {
77: UserContext *user = (UserContext*)ctx;
79: PetscInt i,j,M,N,xm,ym,xs,ys;
80: PetscScalar Hx,Hy,pi,uu,tt;
81: PetscScalar **array;
82: DM da;
85: KSPGetDM(ksp,&da);
86: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
87: uu = user->uu; tt = user->tt;
88: pi = 4*atan(1.0);
89: Hx = 1.0/(PetscReal)(M);
90: Hy = 1.0/(PetscReal)(N);
92: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
93: /* printf(" M N: %d %d; xm ym: %d %d; xs ys: %d %d\n",M,N,xm,ym,xs,ys); */
94: DMDAVecGetArray(da, b, &array);
95: for (j=ys; j<ys+ym; j++) {
96: for (i=xs; i<xs+xm; i++) {
97: array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*cos(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
98: }
99: }
100: DMDAVecRestoreArray(da, b, &array);
101: VecAssemblyBegin(b);
102: VecAssemblyEnd(b);
104: /* force right hand side to be consistent for singular matrix */
105: /* note this is really a hack, normally the model would provide you with a consistent right handside */
106: if (user->bcType == NEUMANN) {
107: MatNullSpace nullspace;
109: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
110: MatNullSpaceRemove(nullspace,b,NULL);
111: MatNullSpaceDestroy(&nullspace);
112: }
113: return(0);
114: }
118: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,MatStructure *str,void *ctx)
119: {
120: UserContext *user = (UserContext*)ctx;
122: PetscInt i, j, M, N, xm, ym, xs, ys, num, numi, numj;
123: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
124: MatStencil row, col[5];
125: DM da;
128: KSPGetDM(ksp,&da);
129: DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
130: Hx = 1.0 / (PetscReal)(M);
131: Hy = 1.0 / (PetscReal)(N);
132: HxdHy = Hx/Hy;
133: HydHx = Hy/Hx;
134: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
135: for (j=ys; j<ys+ym; j++) {
136: for (i=xs; i<xs+xm; i++) {
137: row.i = i; row.j = j;
139: if (i==0 || j==0 || i==M-1 || j==N-1) {
140: if (user->bcType == DIRICHLET) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
141: else if (user->bcType == NEUMANN) {
142: num=0; numi=0; numj=0;
143: if (j!=0) {
144: v[num] = -HxdHy; col[num].i = i; col[num].j = j-1;
145: num++; numj++;
146: }
147: if (i!=0) {
148: v[num] = -HydHx; col[num].i = i-1; col[num].j = j;
149: num++; numi++;
150: }
151: if (i!=M-1) {
152: v[num] = -HydHx; col[num].i = i+1; col[num].j = j;
153: num++; numi++;
154: }
155: if (j!=N-1) {
156: v[num] = -HxdHy; col[num].i = i; col[num].j = j+1;
157: num++; numj++;
158: }
159: v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i; col[num].j = j;
160: num++;
161: MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
162: }
163: } else {
164: v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
165: v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
166: v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
167: v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
168: v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
169: MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
170: }
171: }
172: }
173: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
174: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
175: if (user->bcType == NEUMANN) {
176: MatNullSpace nullspace;
178: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
179: MatSetNullSpace(jac,nullspace);
180: MatNullSpaceDestroy(&nullspace);
181: }
182: return(0);
183: }