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

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


 10: extern PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
 11: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);
 12: extern PetscErrorCode ComputeInitialSolution(DM,Vec);

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

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

 26:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 27:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,-3,2,1,0,&da);
 28:   KSPSetDM(ksp,da);
 29:   KSPSetComputeRHS(ksp,ComputeRHS,NULL);
 30:   KSPSetComputeOperators(ksp,ComputeMatrix,NULL);

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

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

 58:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
 59:   Hx   = 2.0*PETSC_PI / (PetscReal)(mx);
 60:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);

 62:   for (i=xs; i<xs+xm; i++) {
 63:     col[0] = 2*i; col[1] = 2*i + 1;
 64:     val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
 65:     VecSetValues(x,2,col,val,INSERT_VALUES);
 66:   }
 67:   VecAssemblyBegin(x);
 68:   VecAssemblyEnd(x);
 69:   return(0);
 70: }

 74: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
 75: {
 77:   PetscInt       mx;
 78:   PetscScalar    h;
 79:   Vec            x;
 80:   DM             da;

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

 94: PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
 95: {
 97:   PetscInt       i,mx,xm,xs;
 98:   PetscScalar    v[7],Hx;
 99:   MatStencil     row,col[7];
100:   PetscScalar    lambda;
101:   DM             da;

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

117:     row.i = i; row.j = 0; row.k = 0; row.c = 1;
118:     v[0]  = lambda; col[0].i = i-1;   col[0].c = 0;
119:     v[1]  = Hx;     col[1].i = i;   col[1].c = 1;
120:     v[2]  = -lambda;col[2].i = i+1; col[2].c = 0;
121:     MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
122:   }
123:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
124:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
125:   MatView(jac,PETSC_VIEWER_BINARY_(PETSC_COMM_SELF));
126:   return(0);
127: }