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:   PetscScalar    mone = -1.0;
     21:   PetscReal      norm;
     22:   DA             da;

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

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

     31:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);

     33:   ComputeInitialSolution(dmmg);

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

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

     48:   return 0;
     49: }

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

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

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

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

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

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