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: /*T

     20: T*/

     22: #include <petscdm.h>
     23: #include <petscdmda.h>
     24: #include <petscksp.h>

     26: static PetscErrorCode ComputeMatrix(KSP,Mat,Mat,void*);
     27: static PetscErrorCode ComputeRHS(KSP,Vec,void*);

     29: typedef struct {
     30:   PetscInt    k;
     31:   PetscScalar e;
     32: } AppCtx;

     34: int main(int argc,char **argv)
     35: {
     37:   KSP            ksp;
     38:   DM             da;
     39:   AppCtx         user;
     40:   Mat            A;
     41:   Vec            b,b2;
     42:   Vec            x;
     43:   PetscReal      nrm;

     45:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
     46:   user.k = 1;
     47:   user.e = .99;
     48:   PetscOptionsGetInt(NULL,0,"-k",&user.k,0);
     49:   PetscOptionsGetScalar(NULL,0,"-e",&user.e,0);

     51:   KSPCreate(PETSC_COMM_WORLD,&ksp);
     52:   DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,128,1,1,0,&da);
     53:   DMSetFromOptions(da);
     54:   DMSetUp(da);
     55:   KSPSetDM(ksp,da);
     56:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
     57:   KSPSetComputeOperators(ksp,ComputeMatrix,&user);
     58:   KSPSetFromOptions(ksp);
     59:   KSPSolve(ksp,NULL,NULL);

     61:   KSPGetOperators(ksp,&A,NULL);
     62:   KSPGetSolution(ksp,&x);
     63:   KSPGetRhs(ksp,&b);
     64:   VecDuplicate(b,&b2);
     65:   MatMult(A,x,b2);
     66:   VecAXPY(b2,-1.0,b);
     67:   VecNorm(b2,NORM_MAX,&nrm);
     68:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)nrm);

     70:   VecDestroy(&b2);
     71:   KSPDestroy(&ksp);
     72:   DMDestroy(&da);
     73:   PetscFinalize();
     74:   return ierr;
     75: }

     77: static PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
     78: {
     80:   PetscInt       mx,idx[2];
     81:   PetscScalar    h,v[2];
     82:   DM             da;

     85:   KSPGetDM(ksp,&da);
     86:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
     87:   h      = 1.0/((mx-1));
     88:   VecSet(b,h);
     89:   idx[0] = 0; idx[1] = mx -1;
     90:   v[0]   = v[1] = 0.0;
     91:   VecSetValues(b,2,idx,v,INSERT_VALUES);
     92:   VecAssemblyBegin(b);
     93:   VecAssemblyEnd(b);
     94:   return(0);
     95: }

     97: static PetscErrorCode ComputeMatrix(KSP ksp,Mat J,Mat jac,void *ctx)
     98: {
     99:   AppCtx         *user = (AppCtx*)ctx;
    101:   PetscInt       i,mx,xm,xs;
    102:   PetscScalar    v[3],h,xlow,xhigh;
    103:   MatStencil     row,col[3];
    104:   DM             da;

    107:   KSPGetDM(ksp,&da);
    108:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
    109:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
    110:   h    = 1.0/(mx-1);

    112:   for (i=xs; i<xs+xm; i++) {
    113:     row.i = i;
    114:     if (i==0 || i==mx-1) {
    115:       v[0] = 2.0/h;
    116:       MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);
    117:     } else {
    118:       xlow  = h*(PetscReal)i - .5*h;
    119:       xhigh = xlow + h;
    120:       v[0]  = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xlow))/h;col[0].i = i-1;
    121:       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;
    122:       v[2]  = (-1.0 - user->e*PetscSinScalar(2.0*PETSC_PI*user->k*xhigh))/h;col[2].i = i+1;
    123:       MatSetValuesStencil(jac,1,&row,3,col,v,INSERT_VALUES);
    124:     }
    125:   }
    126:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
    127:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
    128:   return(0);
    129: }

    132: /*TEST

    134:    test:
    135:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full
    136:       requires: !single

    138:    test:
    139:       suffix: 2
    140:       nsize: 2
    141:       args: -pc_type mg -ksp_type fgmres -da_refine 2 -ksp_monitor_short -mg_levels_ksp_monitor_short -mg_levels_ksp_norm_type unpreconditioned -ksp_view -pc_mg_type full
    142:       requires: !single

    144: TEST*/