Actual source code: ex66.c
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: Periodic boundary conditions
8: p(x=0) = p(x=1)
9: Neuman boundary conditions
10: dp/dy = 0 for y = 0, y = 1.
12: This example uses the DM_BOUNDARY_MIRROR to implement the Neumann boundary conditions, see the manual pages for DMBoundaryType
14: Compare to ex50.c
15: */
17: static char help[] = "Solves 2D Poisson equation,\n\
18: using mirrored boundaries to implement Neumann boundary conditions,\n\
19: using multigrid.\n\n";
21: #include <petscdm.h>
22: #include <petscdmda.h>
23: #include <petscksp.h>
24: #include <petscsys.h>
25: #include <petscvec.h>
27: extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,void*);
28: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
30: typedef struct {
31: PetscScalar uu, tt;
32: } UserContext;
34: int main(int argc,char **argv)
35: {
36: KSP ksp;
37: DM da;
38: UserContext user;
40: PetscInitialize(&argc,&argv,(char*)0,help);
41: KSPCreate(PETSC_COMM_WORLD,&ksp);
42: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_MIRROR,DMDA_STENCIL_STAR,11,11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
43: DMSetFromOptions(da);
44: DMSetUp(da);
45: KSPSetDM(ksp,(DM)da);
46: DMSetApplicationContext(da,&user);
48: user.uu = 1.0;
49: user.tt = 1.0;
51: KSPSetComputeRHS(ksp,ComputeRHS,&user);
52: KSPSetComputeOperators(ksp,ComputeJacobian,&user);
53: KSPSetFromOptions(ksp);
54: KSPSolve(ksp,NULL,NULL);
56: DMDestroy(&da);
57: KSPDestroy(&ksp);
58: PetscFinalize();
59: return 0;
60: }
62: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
63: {
64: UserContext *user = (UserContext*)ctx;
65: PetscInt i,j,M,N,xm,ym,xs,ys;
66: PetscScalar Hx,Hy,pi,uu,tt;
67: PetscScalar **array;
68: DM da;
69: MatNullSpace nullspace;
72: KSPGetDM(ksp,&da);
73: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
74: uu = user->uu; tt = user->tt;
75: pi = 4*PetscAtanReal(1.0);
76: Hx = 1.0/(PetscReal)(M);
77: Hy = 1.0/(PetscReal)(N);
79: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
80: DMDAVecGetArray(da, b, &array);
81: for (j=ys; j<ys+ym; j++) {
82: for (i=xs; i<xs+xm; i++) {
83: array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*PetscCosScalar(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
84: }
85: }
86: DMDAVecRestoreArray(da, b, &array);
87: VecAssemblyBegin(b);
88: VecAssemblyEnd(b);
90: /* force right hand side to be consistent for singular matrix */
91: /* note this is really a hack, normally the model would provide you with a consistent right handside */
92: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
93: MatNullSpaceRemove(nullspace,b);
94: MatNullSpaceDestroy(&nullspace);
95: return 0;
96: }
98: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
99: {
100: PetscInt i, j, M, N, xm, ym, xs, ys;
101: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
102: MatStencil row, col[5];
103: DM da;
104: MatNullSpace nullspace;
107: KSPGetDM(ksp,&da);
108: DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
109: Hx = 1.0 / (PetscReal)(M);
110: Hy = 1.0 / (PetscReal)(N);
111: HxdHy = Hx/Hy;
112: HydHx = Hy/Hx;
113: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
114: for (j=ys; j<ys+ym; j++) {
115: for (i=xs; i<xs+xm; i++) {
116: row.i = i; row.j = j;
117: v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
118: v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
119: v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
120: v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
121: v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
122: MatSetValuesStencil(jac,1,&row,5,col,v,ADD_VALUES);
123: }
124: }
125: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
126: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
128: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
129: MatSetNullSpace(J,nullspace);
130: MatNullSpaceDestroy(&nullspace);
131: return 0;
132: }
134: /*TEST
136: build:
137: requires: !complex !single
139: test:
140: args: -pc_type mg -pc_mg_type full -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type svd -ksp_view
142: test:
143: suffix: 2
144: nsize: 4
145: args: -pc_type mg -pc_mg_type full -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type redundant -mg_coarse_redundant_pc_type svd -ksp_view
147: TEST*/