学习poisson.c

  1 static char help[] = "A structured-grid Poisson problem with DMDA+KSP.\n\n";
  2 
  3 #include <petsc.h>
  4 
  5 extern PetscErrorCode formMatrix(DM, Mat);
  6 extern PetscErrorCode formExact(DM, Vec);
  7 extern PetscErrorCode formRHS(DM, Vec);
  8 
  9 //STARTMAIN
 10 int main(int argc,char **args) {
 11   PetscErrorCode ierr;
 12   DM            da;
 13   Mat           A;
 14   Vec           b,u,uexact;
 15   KSP           ksp;
 16   double        errnorm;
 17   DMDALocalInfo info;
 18 
 19   PetscInitialize(&argc,&args,(char*)0,help);
 20 
 21   // default size (9 x 9) can be changed using -da_refine X or
 22   //     -da_grid_x M -da_grid_y N
 23   ierr = DMDACreate2d(PETSC_COMM_WORLD,
 24       DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR,
 25       9,9,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da); CHKERRQ(ierr);
 26 
 27   // create linear system matrix A
 28   ierr = DMSetFromOptions(da); CHKERRQ(ierr);
 29   ierr = DMSetUp(da); CHKERRQ(ierr);
 30   ierr = DMCreateMatrix(da,&A); CHKERRQ(ierr);
 31   ierr = MatSetFromOptions(A); CHKERRQ(ierr);
 32 
 33   // create right-hand-side (RHS) b, approx solution u, exact solution uexact
 34   ierr = DMCreateGlobalVector(da,&b); CHKERRQ(ierr);
 35   ierr = VecDuplicate(b,&u); CHKERRQ(ierr);
 36   ierr = VecDuplicate(b,&uexact); CHKERRQ(ierr);
 37 
 38   // fill vectors and assemble linear system
 39   ierr = formExact(da,uexact); CHKERRQ(ierr);
 40   ierr = formRHS(da,b); CHKERRQ(ierr);
 41   ierr = formMatrix(da,A); CHKERRQ(ierr);
 42 
 43   // create and solve the linear system
 44   ierr = KSPCreate(PETSC_COMM_WORLD,&ksp); CHKERRQ(ierr);
 45   ierr = KSPSetOperators(ksp,A,A); CHKERRQ(ierr);
 46   ierr = KSPSetFromOptions(ksp); CHKERRQ(ierr);
 47   ierr = KSPSolve(ksp,b,u); CHKERRQ(ierr);
 48 
 49   // report on grid and numerical error
 50   ierr = VecAXPY(u,-1.0,uexact); CHKERRQ(ierr);    // u <- u + (-1.0) uxact
 51   ierr = VecNorm(u,NORM_INFINITY,&errnorm); CHKERRQ(ierr);
 52   ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr);
 53   ierr = PetscPrintf(PETSC_COMM_WORLD,
 54              "on %d x %d grid:  error |u-uexact|_inf = %g\n",
 55              info.mx,info.my,errnorm); CHKERRQ(ierr);
 56 
 57   VecDestroy(&u);  VecDestroy(&uexact);  VecDestroy(&b);
 58   MatDestroy(&A);  KSPDestroy(&ksp);  DMDestroy(&da);
 59   return PetscFinalize();
 60 }
 61 //ENDMAIN
 62 
 63 //STARTMATRIX
 64 PetscErrorCode formMatrix(DM da, Mat A) {
 65   PetscErrorCode ierr;
 66   DMDALocalInfo info;
 67   MatStencil    row, col[5];
 68   double        hx, hy, v[5];
 69   int           i, j, ncols;
 70 
 71   ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
 72   hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
 73   for (j = info.ys; j < info.ys+info.ym; j++) {
 74     for (i = info.xs; i < info.xs+info.xm; i++) {
 75       row.j = j;           // row of A corresponding to (x_i,y_j)
 76       row.i = i;
 77       col[0].j = j;        // in this diagonal entry
 78       col[0].i = i;
 79       ncols = 1;
 80       if (i==0 || i==info.mx-1 || j==0 || j==info.my-1) {
 81         v[0] = 1.0;  // if on boundary, just insert diagonal entry
 82       } else {
 83         v[0] = 2*(hy/hx + hx/hy); // ... everywhere else we build a row
 84         // if neighbor is NOT a known boundary value then we put an entry
 85         if (i-1 > 0) {
 86           col[ncols].j = j;    col[ncols].i = i-1;  v[ncols++] = -hy/hx;  }
 87         if (i+1 < info.mx-1) {
 88           col[ncols].j = j;    col[ncols].i = i+1;  v[ncols++] = -hy/hx;  }
 89         if (j-1 > 0) {
 90           col[ncols].j = j-1;  col[ncols].i = i;    v[ncols++] = -hx/hy;  }
 91         if (j+1 < info.my-1) {
 92           col[ncols].j = j+1;  col[ncols].i = i;    v[ncols++] = -hx/hy;  }
 93       }
 94       ierr = MatSetValuesStencil(A,1,&row,ncols,col,v,INSERT_VALUES); CHKERRQ(ierr);
 95     }
 96   }
 97   ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
 98   ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
 99   return 0;
100 }
101 //ENDMATRIX
102 
103 //STARTEXACT
104 PetscErrorCode formExact(DM da, Vec uexact) {
105   PetscErrorCode ierr;
106   DMDALocalInfo info;
107   int           i, j;
108   double        hx, hy, x, y, **auexact;
109 
110   ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
111   hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
112   ierr = DMDAVecGetArray(da, uexact, &auexact);CHKERRQ(ierr);
113   for (j = info.ys; j < info.ys+info.ym; j++) {
114     y = j * hy;
115     for (i = info.xs; i < info.xs+info.xm; i++) {
116       x = i * hx;
117       auexact[j][i] = x*x * (1.0 - x*x) * y*y * (y*y - 1.0);
118     }
119   }
120   ierr = DMDAVecRestoreArray(da, uexact, &auexact);CHKERRQ(ierr);
121   return 0;
122 }
123 
124 PetscErrorCode formRHS(DM da, Vec b) {
125   PetscErrorCode ierr;
126   int           i, j;
127   double        hx, hy, x, y, f, **ab;
128   DMDALocalInfo info;
129 
130   ierr = DMDAGetLocalInfo(da,&info); CHKERRQ(ierr);
131   hx = 1.0/(info.mx-1);  hy = 1.0/(info.my-1);
132   ierr = DMDAVecGetArray(da, b, &ab);CHKERRQ(ierr);
133   for (j=info.ys; j<info.ys+info.ym; j++) {
134     y = j * hy;
135     for (i=info.xs; i<info.xs+info.xm; i++) {
136       x = i * hx;
137       if (i==0 || i==info.mx-1 || j==0 || j==info.my-1) {
138         ab[j][i] = 0.0;                    // on bdry the eqn is 1*u = 0
139       } else {  // if not bdry; note  f = - (u_xx + u_yy)  where u is exact
140         f = 2.0 * ( (1.0 - 6.0*x*x) * y*y * (1.0 - y*y)
141                     + (1.0 - 6.0*y*y) * x*x * (1.0 - x*x) );
142         ab[j][i] = hx * hy * f;
143       }
144     }
145   }
146   ierr = DMDAVecRestoreArray(da, b, &ab); CHKERRQ(ierr);
147   return 0;
148 }
149 //ENDEXACT