Actual source code: ex28.c
3: static char help[] = "Solves 1D wave equation using multigrid.\n\n";
5: #include petscda.h
6: #include petscsles.h
8: extern int ComputeJacobian(DMMG,Mat);
9: extern int ComputeRHS(DMMG,Vec);
10: extern int ComputeInitialSolution(DMMG*);
14: int main(int argc,char **argv)
15: {
16: int ierr,i;
17: DMMG *dmmg;
18: PetscScalar mone = -1.0;
19: PetscReal norm;
20: DA da;
22: PetscInitialize(&argc,&argv,(char *)0,help);
24: DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
25: DACreate1d(PETSC_COMM_WORLD,DA_XPERIODIC,-3,2,1,0,&da);
26: DMMGSetDM(dmmg,(DM)da);
27: DADestroy(da);
29: DMMGSetSLES(dmmg,ComputeRHS,ComputeJacobian);
31: ComputeInitialSolution(dmmg);
33: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
34: for (i=0; i<1000; i++) {
35: DMMGSolve(dmmg);
36: VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
37: }
44: MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
45: VecAXPY(&mone,DMMGGetb(dmmg),DMMGGetr(dmmg));
46: VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
47: /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",norm); */
49: DMMGDestroy(dmmg);
50: PetscFinalize();
52: return 0;
53: }
57: int ComputeInitialSolution(DMMG *dmmg)
58: {
59: int ierr,mx,col[2],xs,xm,i;
60: PetscScalar Hx,val[2];
61: Vec x = DMMGGetx(dmmg);
64: DAGetInfo(DMMGGetDA(dmmg),0,&mx,0,0,0,0,0,0,0,0,0);
65: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
66: DAGetCorners(DMMGGetDA(dmmg),&xs,0,0,&xm,0,0);
67:
68: for(i=xs; i<xs+xm; i++){
69: col[0] = 2*i; col[1] = 2*i + 1;
70: val[0] = val[1] = PetscSinScalar(i*Hx);
71: VecSetValues(x,2,col,val,INSERT_VALUES);
72: }
73: VecAssemblyBegin(x);
74: VecAssemblyEnd(x);
75: return(0);
76: }
77:
80: int ComputeRHS(DMMG dmmg,Vec b)
81: {
82: int ierr,mx;
83: PetscScalar h;
86: DAGetInfo((DA)dmmg->dm,0,&mx,0,0,0,0,0,0,0,0,0);
87: h = 2.0*PETSC_PI/((mx));
88: VecCopy(dmmg->x,b);
89: VecScale(&h,b);
90: return(0);
91: }
95: int ComputeJacobian(DMMG dmmg,Mat jac)
96: {
97: DA da = (DA)dmmg->dm;
98: int ierr,i,j,k,mx,xm,xs;
99: PetscScalar v[7],Hx;
100: MatStencil row,col[7];
101: PetscScalar lambda;
103: PetscMemzero(col,7*sizeof(MatStencil));
104: DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
105: Hx = 2.0*PETSC_PI / (PetscReal)(mx);
106: DAGetCorners(da,&xs,0,0,&xm,0,0);
107: lambda = 2*Hx;
108: for(i=xs; i<xs+xm; i++){
109: row.i = i; row.j = 0; row.k = 0; row.c = 0;
110: v[0] = Hx; col[0].i = i; col[0].c = 0;
111: v[1] = lambda; col[1].i = i-1; col[1].c = 1;
112: v[2] = -lambda;col[2].i = i+1; col[2].c = 1;
113: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
115: row.i = i; row.j = 0; row.k = 0; row.c = 1;
116: v[0] = lambda; col[0].i = i-1; col[0].c = 0;
117: v[1] = Hx; col[1].i = i; col[1].c = 1;
118: v[2] = -lambda;col[2].i = i+1; col[2].c = 0;
119: MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
120: }
121: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
122: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
123: MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
124: return 0;
125: }