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: {
 36:   KSP            ksp;
 37:   DM             da;
 38:   AppCtx         user;
 39:   Mat            A;
 40:   Vec            b,b2;
 41:   Vec            x;
 42:   PetscReal      nrm;

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

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

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

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

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

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

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

104:   KSPGetDM(ksp,&da);
105:   DMDAGetInfo(da,0,&mx,0,0,0,0,0,0,0,0,0,0,0);
106:   DMDAGetCorners(da,&xs,0,0,&xm,0,0);
107:   h    = 1.0/(mx-1);

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

128: /*TEST

130:    test:
131:       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
132:       requires: !single

134:    test:
135:       suffix: 2
136:       nsize: 2
137:       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
138:       requires: !single

140: TEST*/