Actual source code: ex66.c
petsc-3.15.0 2021-04-05
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;
41: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
42: KSPCreate(PETSC_COMM_WORLD,&ksp);
43: DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_MIRROR,DMDA_STENCIL_STAR,11,11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
44: DMSetFromOptions(da);
45: DMSetUp(da);
46: KSPSetDM(ksp,(DM)da);
47: DMSetApplicationContext(da,&user);
49: user.uu = 1.0;
50: user.tt = 1.0;
52: KSPSetComputeRHS(ksp,ComputeRHS,&user);
53: KSPSetComputeOperators(ksp,ComputeJacobian,&user);
54: KSPSetFromOptions(ksp);
55: KSPSolve(ksp,NULL,NULL);
57: DMDestroy(&da);
58: KSPDestroy(&ksp);
59: PetscFinalize();
60: return ierr;
61: }
63: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
64: {
65: UserContext *user = (UserContext*)ctx;
67: PetscInt i,j,M,N,xm,ym,xs,ys;
68: PetscScalar Hx,Hy,pi,uu,tt;
69: PetscScalar **array;
70: DM da;
71: MatNullSpace nullspace;
74: KSPGetDM(ksp,&da);
75: DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
76: uu = user->uu; tt = user->tt;
77: pi = 4*PetscAtanReal(1.0);
78: Hx = 1.0/(PetscReal)(M);
79: Hy = 1.0/(PetscReal)(N);
81: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
82: DMDAVecGetArray(da, b, &array);
83: for (j=ys; j<ys+ym; j++) {
84: for (i=xs; i<xs+xm; i++) {
85: array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*PetscCosScalar(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
86: }
87: }
88: DMDAVecRestoreArray(da, b, &array);
89: VecAssemblyBegin(b);
90: VecAssemblyEnd(b);
92: /* force right hand side to be consistent for singular matrix */
93: /* note this is really a hack, normally the model would provide you with a consistent right handside */
94: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
95: MatNullSpaceRemove(nullspace,b);
96: MatNullSpaceDestroy(&nullspace);
97: return(0);
98: }
100: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
101: {
103: PetscInt i, j, M, N, xm, ym, xs, ys;
104: PetscScalar v[5], Hx, Hy, HydHx, HxdHy;
105: MatStencil row, col[5];
106: DM da;
107: MatNullSpace nullspace;
110: KSPGetDM(ksp,&da);
111: DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
112: Hx = 1.0 / (PetscReal)(M);
113: Hy = 1.0 / (PetscReal)(N);
114: HxdHy = Hx/Hy;
115: HydHx = Hy/Hx;
116: DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
117: for (j=ys; j<ys+ym; j++) {
118: for (i=xs; i<xs+xm; i++) {
119: row.i = i; row.j = j;
120: v[0] = -HxdHy; col[0].i = i; col[0].j = j-1;
121: v[1] = -HydHx; col[1].i = i-1; col[1].j = j;
122: v[2] = 2.0*(HxdHy + HydHx); col[2].i = i; col[2].j = j;
123: v[3] = -HydHx; col[3].i = i+1; col[3].j = j;
124: v[4] = -HxdHy; col[4].i = i; col[4].j = j+1;
125: MatSetValuesStencil(jac,1,&row,5,col,v,ADD_VALUES);
126: }
127: }
128: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
129: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
131: MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
132: MatSetNullSpace(J,nullspace);
133: MatNullSpaceDestroy(&nullspace);
134: return(0);
135: }
139: /*TEST
141: build:
142: requires: !complex !single
144: test:
145: args: -pc_type mg -pc_mg_type full -ksp_monitor_short -da_refine 3 -mg_coarse_pc_type svd -ksp_view
147: test:
148: suffix: 2
149: nsize: 4
150: 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
152: TEST*/