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:   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:   DMMGSetKSP(dmmg,ComputeRHS,ComputeJacobian);

     55:   DMMGSolve(dmmg);

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

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

     65:   return 0;
     66: }

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

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

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

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