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
 7:  #include petscdmmg.h


 15: int main(int argc,char **argv)
 16: {
 18:   PetscInt       i;
 19:   DMMG           *dmmg;
 20:   PetscReal      norm;
 21:   DA             da;

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

 25:   DMMGCreate(PETSC_COMM_WORLD,3,PETSC_NULL,&dmmg);
 26:   DACreate1d(PETSC_COMM_WORLD,DA_XPERIODIC,-3,2,1,0,&da);
 27:   DMMGSetDM(dmmg,(DM)da);
 28:   DADestroy(da);

 30:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);

 32:   ComputeInitialSolution(dmmg);

 34:   VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
 35:   for (i=0; i<1000; i++) {
 36:     DMMGSolve(dmmg);
 37:     VecView(DMMGGetx(dmmg),PETSC_VIEWER_DRAW_WORLD);
 38:   }
 39:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
 40:   VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
 41:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
 42:   /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */

 44:   DMMGDestroy(dmmg);
 45:   PetscFinalize();

 47:   return 0;
 48: }

 52: PetscErrorCode ComputeInitialSolution(DMMG *dmmg)
 53: {
 55:   PetscInt       mx,col[2],xs,xm,i;
 56:   PetscScalar    Hx,val[2];
 57:   Vec            x = DMMGGetx(dmmg);

 60:   DAGetInfo(DMMGGetDA(dmmg),0,&mx,0,0,0,0,0,0,0,0,0);
 61:   Hx = 2.0*PETSC_PI / (PetscReal)(mx);
 62:   DAGetCorners(DMMGGetDA(dmmg),&xs,0,0,&xm,0,0);
 63: 
 64:   for(i=xs; i<xs+xm; i++){
 65:     col[0] = 2*i; col[1] = 2*i + 1;
 66:     val[0] = val[1] = PetscSinScalar(((PetscScalar)i)*Hx);
 67:     VecSetValues(x,2,col,val,INSERT_VALUES);
 68:   }
 69:   VecAssemblyBegin(x);
 70:   VecAssemblyEnd(x);
 71:   return(0);
 72: }
 73: 
 76: PetscErrorCode ComputeRHS(DMMG dmmg,Vec b)
 77: {
 79:   PetscInt       mx;
 80:   PetscScalar    h;

 83:   DAGetInfo((DA)dmmg->dm,0,&mx,0,0,0,0,0,0,0,0,0);
 84:   h    = 2.0*PETSC_PI/((mx));
 85:   VecCopy(dmmg->x,b);
 86:   VecScale(b,h);
 87:   return(0);
 88: }

 92: PetscErrorCode ComputeJacobian(DMMG dmmg,Mat J,Mat jac)
 93: {
 94:   DA             da = (DA)dmmg->dm;
 96:   PetscInt       i,mx,xm,xs;
 97:   PetscScalar    v[7],Hx;
 98:   MatStencil     row,col[7];
 99:   PetscScalar    lambda;

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

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