Actual source code: ex28.c


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

 5:  #include petscda.h
 6:  #include petscksp.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:   DMMGSetKSP(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: }