Actual source code: ex25.c

  2: /*
  3:  Partial differential equation

  5:    d  (1 + e*sine(2*pi*k*x)) d u = 1, 0 < x < 1,
  6:    --                        ---
  7:    dx                        dx
  8: with boundary conditions

 10:    u = 0 for x = 0, x = 1

 12:    This uses multigrid to solve the linear system

 14: */

 16: static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.\n\n";

 18:  #include petscda.h
 19:  #include petscksp.h
 20:  #include petscdmmg.h


 25: typedef struct {
 26:   PetscInt    k;
 27:   PetscScalar e;
 28: } AppCtx;

 32: int main(int argc,char **argv)
 33: {
 35:   DMMG           *dmmg;
 36:   PetscReal      norm;
 37:   DA             da;
 38:   AppCtx         user;

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

 42:   user.k = 1;
 43:   user.e = .99;
 44:   PetscOptionsGetInt(0,"-k",&user.k,0);
 45:   PetscOptionsGetScalar(0,"-e",&user.e,0);

 47:   DMMGCreate(PETSC_COMM_WORLD,3,&user,&dmmg);
 48:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,-3,1,1,0,&da);
 49:   DMMGSetDM(dmmg,(DM)da);
 50:   DADestroy(da);

 52:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);

 54:   DMMGSolve(dmmg);

 56:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
 57:   VecAXPY(DMMGGetr(dmmg),-1.0,DMMGGetRHS(dmmg));
 58:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
 59:   /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %G\n",norm); */

 61:   DMMGDestroy(dmmg);
 62:   PetscFinalize();

 64:   return 0;
 65: }

 69: PetscErrorCode ComputeRHS(DMMG dmmg,Vec b)
 70: {
 72:   PetscInt       mx,idx[2];
 73:   PetscScalar    h,v[2];

 76:   DAGetInfo((DA)dmmg->dm,0,&mx,0,0,0,0,0,0,0,0,0);
 77:   h      = 1.0/((mx-1));
 78:   VecSet(b,h);
 79:   idx[0] = 0; idx[1] = mx -1;
 80:   v[0]   = v[1] = 0.0;
 81:   VecSetValues(b,2,idx,v,INSERT_VALUES);
 82:   VecAssemblyBegin(b);
 83:   VecAssemblyEnd(b);
 84:   return(0);
 85: }
 86: 
 89: PetscErrorCode ComputeJacobian(DMMG dmmg,Mat J,Mat jac)
 90: {
 91:   DA             da = (DA)dmmg->dm;
 93:   PetscInt       i,mx,xm,xs;
 94:   PetscScalar    v[3],h,xlow,xhigh;
 95:   MatStencil     row,col[3];
 96:   AppCtx         *user = (AppCtx*)dmmg->user;

 98:   DAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0);
 99:   DAGetCorners(da,&xs,0,0,&xm,0,0);
100:   h    = 1.0/(mx-1);

102:   for(i=xs; i<xs+xm; i++){
103:     row.i = i;
104:     if (i==0 || i==mx-1){
105:       v[0] = 2.0;
106:       MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
107:     } else {
108:        xlow  = h*(PetscReal)i - .5*h;
109:        xhigh = xlow + h;
110:        v[0] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow))/h;col[0].i = i-1;
111:        v[1] = (2.0 + user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow) + user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[1].i = row.i;
112:        v[2] = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1;
113:       MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
114:     }
115:   }
116:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
117:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
118:   return 0;
119: }