2: /*$Id: ex22.c,v 1.19 2001/08/07 21:30:54 bsmith Exp $*/
  3: /*
  4:  Partial differential equation

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

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

 13:    This uses multigrid to solve the linear system

 15: */

 17: static char help[] = "Solves 1D variable coefficient Laplacian using multigrid.nn";

 19:  #include petscda.h
 20:  #include petscsles.h

 22: extern int ComputeJacobian(DMMG,Mat);
 23: extern int ComputeRHS(DMMG,Vec);

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

 30: #undef __FUNCT__
 32: int main(int argc,char **argv)
 33: {
 34:   int         ierr;
 35:   DMMG        *dmmg;
 36:   PetscScalar mone = -1.0;
 37:   PetscReal   norm;
 38:   DA          da;
 39:   AppCtx      user;

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

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

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

 53:   DMMGSetSLES(dmmg,ComputeRHS,ComputeJacobian);

 55:   DMMGSolve(dmmg);

 57:   MatMult(DMMGGetJ(dmmg),DMMGGetx(dmmg),DMMGGetr(dmmg));
 58:   VecAXPY(&mone,DMMGGetb(dmmg),DMMGGetr(dmmg));
 59:   VecNorm(DMMGGetr(dmmg),NORM_2,&norm);
 60:   /* PetscPrintf(PETSC_COMM_WORLD,"Residual norm %gn",norm); */

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

 65:   return 0;
 66: }

 68: #undef __FUNCT__
 70: int ComputeRHS(DMMG dmmg,Vec b)
 71: {
 72:   int         ierr,mx,idx[2];
 73:   PetscScalar h,v[2];

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

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

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