3: static char help[] = "Solves 1D wave equation using multigrid.\n\n";

  5: #include <petscdmda.h>
  6: #include <petscksp.h>


  9: extern PetscErrorCode ComputeMatrix(DM,Vec,Mat,Mat,MatStructure*);
 10: extern PetscErrorCode ComputeRHS(DM,Vec,Vec);
 11: extern PetscErrorCode ComputeInitialSolution(DM,Vec);

 15: int main(int argc,char **argv)
 16: {
 18:   PetscInt       i;
 19:   KSP            ksp;
 20:   DM             da;
 21:   Vec            x;

 23:   PetscInitialize(&argc,&argv,(char *)0,help);

 25:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 26:   DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,-3,2,1,0,&da);
 27:   KSPSetDM(ksp,(DM)da);
 28:   DMSetFunction(da,ComputeRHS);
 29:   DMSetJacobian(da,ComputeMatrix);

 31:   KSPSetFromOptions(ksp);
 32:   DMCreateGlobalVector(da,&x);
 33:   ComputeInitialSolution(da,x);
 34:   DMSetApplicationContext(da,x);
 35:   KSPSetUp(ksp);
 36:   VecView(x,PETSC_VIEWER_DRAW_WORLD);
 37:   for (i=0; i<10; i++) {
 38:     KSPSolve(ksp,PETSC_NULL,x);
 39:     VecView(x,PETSC_VIEWER_DRAW_WORLD);
 40:   }
 41:   VecDestroy(&x);
 42:   KSPDestroy(&ksp);
 43:   DMDestroy(&da);
 44:   PetscFinalize();
 45:   return 0;
 46: }

 50: PetscErrorCode ComputeInitialSolution(DM da,Vec x)
 51: {
 53:   PetscInt       mx,col[2],xs,xm,i;
 54:   PetscScalar    Hx,val[2];

 57:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
 58:   Hx = 2.0*PETSC_PI / (PetscReal)(mx);
 59:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
 60: 
 61:   for(i=xs; i<xs+xm; i++){
 62:     col[0] = 2*i; col[1] = 2*i + 1;
 63:     val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
 64:     VecSetValues(x,2,col,val,INSERT_VALUES);
 65:   }
 66:   VecAssemblyBegin(x);
 67:   VecAssemblyEnd(x);
 68:   return(0);
 69: }
 70: 
 73: PetscErrorCode ComputeRHS(DM da,Vec dummy,Vec b)
 74: {
 76:   PetscInt       mx;
 77:   PetscScalar    h;
 78:   Vec            x;

 81:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
 82:   DMGetApplicationContext(da,&x);
 83:   h    = 2.0*PETSC_PI/((mx));
 84:   VecCopy(x,b);
 85:   VecScale(b,h);
 86:   return(0);
 87: }

 91: PetscErrorCode ComputeMatrix(DM da,Vec x,Mat J,Mat jac,MatStructure *str)
 92: {
 94:   PetscInt       i,mx,xm,xs;
 95:   PetscScalar    v[7],Hx;
 96:   MatStencil     row,col[7];
 97:   PetscScalar    lambda;

 99:   PetscMemzero(col,7*sizeof(MatStencil));
100:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
101:   Hx = 2.0*PETSC_PI / (PetscReal)(mx);
102:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
103:   lambda = 2.0*Hx;
104:   for(i=xs; i<xs+xm; i++){
105:     row.i = i; row.j = 0; row.k = 0; row.c = 0;
106:     v[0] = Hx;     col[0].i = i;   col[0].c = 0;
107:     v[1] = lambda; col[1].i = i-1;   col[1].c = 1;
108:     v[2] = -lambda;col[2].i = i+1; col[2].c = 1;
109:     MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);

111:     row.i = i; row.j = 0; row.k = 0; row.c = 1;
112:     v[0] = lambda; col[0].i = i-1;   col[0].c = 0;
113:     v[1] = Hx;     col[1].i = i;   col[1].c = 1;
114:     v[2] = -lambda;col[2].i = i+1; col[2].c = 0;
115:     MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
116:   }
117:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
118:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
119:   MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
120:   return 0;
121: }