Actual source code: ex50.c
petsc-3.7.5 2017-01-01
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);
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 0;
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: }