Actual source code: ex50.c

petsc-master 2016-09-23
Report Typos and Errors
1: /*   DMDA/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 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 3 -ksp_monitor -ksp_view -dm_view draw -draw_pause -1
16:           ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg  -pc_mg_levels 1 -mg_levels_0_pc_type ilu -mg_levels_0_pc_factor_levels 1 -ksp_monitor -ksp_view
17:           ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor
18:           mpiexec -n 4 ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 10 -ksp_monitor -ksp_view -log_summary
19: */

21: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";

23:  #include <petscdm.h>
24:  #include <petscdmda.h>
25:  #include <petscksp.h>
26:  #include <petscsys.h>
27:  #include <petscvec.h>

29: extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,void*);
30: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
31: extern PetscErrorCode ComputeTrueSolution(DM, Vec);
32: extern PetscErrorCode VecView_VTK(Vec, const char [], const char []);

34: typedef enum {DIRICHLET, NEUMANN} BCType;

36: typedef struct {
37:   PetscScalar uu, tt;
38:   BCType      bcType;
39: } UserContext;

43: int main(int argc,char **argv)
44: {
45:   KSP            ksp;
46:   DM             da;
47:   UserContext    user;
48:   PetscInt       bc;

51:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
52:   KSPCreate(PETSC_COMM_WORLD,&ksp);
53:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-11,-11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
54:   KSPSetDM(ksp,(DM)da);
55:   DMSetApplicationContext(da,&user);

57:   user.uu     = 1.0;
58:   user.tt     = 1.0;
59:   bc          = (PetscInt)NEUMANN; /* Use Neumann Boundary Conditions */
60:   user.bcType = (BCType)bc;

63:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
64:   KSPSetComputeOperators(ksp,ComputeJacobian,&user);
65:   KSPSetFromOptions(ksp);
66:   KSPSolve(ksp,NULL,NULL);

68:   DMDestroy(&da);
69:   KSPDestroy(&ksp);
70:   PetscFinalize();
71:   return ierr;
72: }

76: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
77: {
78:   UserContext    *user = (UserContext*)ctx;
80:   PetscInt       i,j,M,N,xm,ym,xs,ys;
81:   PetscScalar    Hx,Hy,pi,uu,tt;
82:   PetscScalar    **array;
83:   DM             da;

86:   KSPGetDM(ksp,&da);
87:   DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
88:   uu   = user->uu; tt = user->tt;
89:   pi   = 4*atan(1.0);
90:   Hx   = 1.0/(PetscReal)(M);
91:   Hy   = 1.0/(PetscReal)(N);

93:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
94:   /* printf(" M N: %d %d; xm ym: %d %d; xs ys: %d %d\n",M,N,xm,ym,xs,ys); */
95:   DMDAVecGetArray(da, b, &array);
96:   for (j=ys; j<ys+ym; j++) {
97:     for (i=xs; i<xs+xm; i++) {
98:       array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*PetscCosScalar(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
99:     }
100:   }
101:   DMDAVecRestoreArray(da, b, &array);
102:   VecAssemblyBegin(b);
103:   VecAssemblyEnd(b);

105:   /* force right hand side to be consistent for singular matrix */
106:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
107:   if (user->bcType == NEUMANN) {
108:     MatNullSpace nullspace;

110:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
111:     MatNullSpaceRemove(nullspace,b);
112:     MatNullSpaceDestroy(&nullspace);
113:   }
114:   return(0);
115: }

119: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
120: {
121:   UserContext    *user = (UserContext*)ctx;
123:   PetscInt       i, j, M, N, xm, ym, xs, ys, num, numi, numj;
124:   PetscScalar    v[5], Hx, Hy, HydHx, HxdHy;
125:   MatStencil     row, col[5];
126:   DM             da;

129:   KSPGetDM(ksp,&da);
130:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
131:   Hx    = 1.0 / (PetscReal)(M);
132:   Hy    = 1.0 / (PetscReal)(N);
133:   HxdHy = Hx/Hy;
134:   HydHx = Hy/Hx;
135:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
136:   for (j=ys; j<ys+ym; j++) {
137:     for (i=xs; i<xs+xm; i++) {
138:       row.i = i; row.j = j;

140:       if (i==0 || j==0 || i==M-1 || j==N-1) {
141:         if (user->bcType == DIRICHLET) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Dirichlet boundary conditions not supported !\n");
142:         else if (user->bcType == NEUMANN) {
143:           num=0; numi=0; numj=0;
144:           if (j!=0) {
145:             v[num] = -HxdHy;              col[num].i = i;   col[num].j = j-1;
146:             num++; numj++;
147:           }
148:           if (i!=0) {
149:             v[num] = -HydHx;              col[num].i = i-1; col[num].j = j;
150:             num++; numi++;
151:           }
152:           if (i!=M-1) {
153:             v[num] = -HydHx;              col[num].i = i+1; col[num].j = j;
154:             num++; numi++;
155:           }
156:           if (j!=N-1) {
157:             v[num] = -HxdHy;              col[num].i = i;   col[num].j = j+1;
158:             num++; numj++;
159:           }
160:           v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i;   col[num].j = j;
161:           num++;
162:           MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
163:         }
164:       } else {
165:         v[0] = -HxdHy;              col[0].i = i;   col[0].j = j-1;
166:         v[1] = -HydHx;              col[1].i = i-1; col[1].j = j;
167:         v[2] = 2.0*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
168:         v[3] = -HydHx;              col[3].i = i+1; col[3].j = j;
169:         v[4] = -HxdHy;              col[4].i = i;   col[4].j = j+1;
170:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
171:       }
172:     }
173:   }
174:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
175:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
176:   if (user->bcType == NEUMANN) {
177:     MatNullSpace nullspace;

179:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
180:     MatSetNullSpace(J,nullspace);
181:     MatNullSpaceDestroy(&nullspace);
182:   }
183:   return(0);
184: }