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: {
15: PetscInt i;
16: KSP ksp;
17: DM da;
18: Vec x;
20: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
21: KSPCreate(PETSC_COMM_WORLD,&ksp);
22: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,3,2,1,0,&da);
23: DMSetFromOptions(da);
24: DMSetUp(da);
25: KSPSetDM(ksp,da);
26: KSPSetComputeRHS(ksp,ComputeRHS,NULL);
27: KSPSetComputeOperators(ksp,ComputeMatrix,NULL);
29: KSPSetFromOptions(ksp);
30: DMCreateGlobalVector(da,&x);
31: ComputeInitialSolution(da,x);
32: DMSetApplicationContext(da,x);
33: KSPSetUp(ksp);
34: VecView(x,PETSC_VIEWER_DRAW_WORLD);
35: for (i=0; i<10; i++) {
36: KSPSolve(ksp,NULL,x);
37: VecView(x,PETSC_VIEWER_DRAW_WORLD);
38: }
39: VecDestroy(&x);
40: KSPDestroy(&ksp);
41: DMDestroy(&da);
42: PetscFinalize();
43: return ierr;
44: }
46: PetscErrorCode ComputeInitialSolution(DM da,Vec x)
47: {
49: PetscInt mx,col[2],xs,xm,i;
50: PetscScalar Hx,val[2];
53: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
54: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
55: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
57: for (i=xs; i<xs+xm; i++) {
58: col[0] = 2*i; col[1] = 2*i + 1;
59: val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
60: VecSetValues(x,2,col,val,INSERT_VALUES);
61: }
62: VecAssemblyBegin(x);
63: VecAssemblyEnd(x);
64: return(0);
65: }
67: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
68: {
70: PetscInt mx;
71: PetscScalar h;
72: Vec x;
73: DM da;
76: KSPGetDM(ksp,&da);
77: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
78: DMGetApplicationContext(da,&x);
79: h = 2.0*PETSC_PI/((mx));
80: VecCopy(x,b);
81: VecScale(b,h);
82: return(0);
83: }
85: PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
86: {
88: PetscInt i,mx,xm,xs;
89: PetscScalar v[7],Hx;
90: MatStencil row,col[7];
91: PetscScalar lambda;
92: DM da;
95: KSPGetDM(ksp,&da);
96: PetscArrayzero(col,7);
97: DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
98: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
99: DMDAGetCorners(da,&xs,0,0,&xm,0,0);
100: lambda = 2.0*Hx;
101: for (i=xs; i<xs+xm; i++) {
102: row.i = i; row.j = 0; row.k = 0; row.c = 0;
103: v[0] = Hx; col[0].i = i; col[0].c = 0;
104: v[1] = lambda; col[1].i = i-1; col[1].c = 1;
105: v[2] = -lambda;col[2].i = i+1; col[2].c = 1;
106: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
108: row.i = i; row.j = 0; row.k = 0; row.c = 1;
109: v[0] = lambda; col[0].i = i-1; col[0].c = 0;
110: v[1] = Hx; col[1].i = i; col[1].c = 1;
111: v[2] = -lambda;col[2].i = i+1; col[2].c = 0;
112: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
113: }
114: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
116: MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
117: return(0);
118: }
120: /*TEST
122: test:
123: 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
125: TEST*/