Actual source code: ex28.c
2: static char help[] = "Solves 1D wave equation using multigrid.\n\n";
4: #include <petscdm.h>
5: #include <petscdmda.h>
6: #include <petscksp.h>
8: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
9: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
10: extern PetscErrorCode ComputeInitialSolution(DM,Vec);
12: int main(int argc,char **argv)
13: {
14: PetscInt i;
15: KSP ksp;
16: DM da;
17: Vec x;
19: PetscInitialize(&argc,&argv,(char*)0,help);
20: KSPCreate(PETSC_COMM_WORLD,&ksp);
21: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,3,2,1,0,&da);
22: DMSetFromOptions(da);
23: DMSetUp(da);
24: KSPSetDM(ksp,da);
25: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
26: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
28: KSPSetFromOptions(ksp);
29: DMCreateGlobalVector(da,&x);
30: ComputeInitialSolution(da,x);
31: DMSetApplicationContext(da,x);
32: KSPSetUp(ksp);
33: VecView(x,PETSC_VIEWER_DRAW_WORLD);
34: for (i=0; i<10; i++) {
35: KSPSolve(ksp,NULL,x);
36: VecView(x,PETSC_VIEWER_DRAW_WORLD);
37: }
38: VecDestroy(&x);
39: KSPDestroy(&ksp);
40: DMDestroy(&da);
41: PetscFinalize();
42: return 0;
43: }
45: PetscErrorCode ComputeInitialSolution(DM da,Vec x)
46: {
47: PetscInt mx,col[2],xs,xm,i;
48: PetscScalar Hx,val[2];
51: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
52: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
53: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
55: for (i=xs; i<xs+xm; i++) {
56: col[0] = 2*i; col[1] = 2*i + 1;
57: val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
58: VecSetValues(x,2,col,val,INSERT_VALUES);
59: }
60: VecAssemblyBegin(x);
61: VecAssemblyEnd(x);
62: return 0;
63: }
65: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
66: {
67: PetscInt mx;
68: PetscScalar h;
69: Vec x;
70: DM da;
73: KSPGetDM(ksp,&da);
74: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
75: DMGetApplicationContext(da,&x);
76: h = 2.0*PETSC_PI/((mx));
77: VecCopy(x,b);
78: VecScale(b,h);
79: return 0;
80: }
82: PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
83: {
84: PetscInt i,mx,xm,xs;
85: PetscScalar v[7],Hx;
86: MatStencil row,col[7];
87: PetscScalar lambda;
88: DM da;
91: KSPGetDM(ksp,&da);
92: PetscArrayzero(col,7);
93: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
94: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
95: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
96: lambda = 2.0*Hx;
97: for (i=xs; i<xs+xm; i++) {
98: row.i = i; row.j = 0; row.k = 0; row.c = 0;
99: v[0] = Hx; col[0].i = i; col[0].c = 0;
100: v[1] = lambda; col[1].i = i-1; col[1].c = 1;
101: v[2] = -lambda;col[2].i = i+1; col[2].c = 1;
102: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
104: row.i = i; row.j = 0; row.k = 0; row.c = 1;
105: v[0] = lambda; col[0].i = i-1; col[0].c = 0;
106: v[1] = Hx; col[1].i = i; col[1].c = 1;
107: v[2] = -lambda;col[2].i = i+1; col[2].c = 0;
108: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
109: }
110: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
111: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
112: MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
113: return 0;
114: }
116: /*TEST
118: test:
119: args: -ksp_monitor_short -pc_type mg -pc_mg_type full -ksp_type fgmres -da_refine 2 -mg_levels_ksp_type gmres -mg_levels_ksp_max_it 1 -mg_levels_pc_type ilu
121: TEST*/