Actual source code: ex50.c

petsc-master 2017-07-26
Report Typos and Errors
  1: /*   DMDA/KSP solving a system of linear equations.
  2:      Poisson equation in 2D:

  4:      div(grad p) = f,  0 < x,y < 1
  5:      with
  6:        forcing function f = -cos(m*pi*x)*cos(n*pi*y),
  7:        Neuman boundary conditions
  8:         dp/dx = 0 for x = 0, x = 1.
  9:         dp/dy = 0 for y = 0, y = 1.

 11:      Contributed by Michael Boghosian <boghmic@iit.edu>, 2008,
 12:          based on petsc/src/ksp/ksp/examples/tutorials/ex29.c and ex32.c

 14:      Compare to ex66.c

 16:      Example of Usage:
 17:           ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 3 -ksp_monitor -ksp_view -dm_view draw -draw_pause -1
 18:           ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg  -pc_mg_levels 1 -mg_levels_0_pc_type ilu -mg_levels_0_pc_factor_levels 1 -ksp_monitor -ksp_view
 19:           ./ex50 -da_grid_x 100 -da_grid_y 100 -pc_type mg -pc_mg_levels 1 -mg_levels_0_pc_type lu -mg_levels_0_pc_factor_shift_type NONZERO -ksp_monitor
 20:           mpiexec -n 4 ./ex50 -da_grid_x 3 -da_grid_y 3 -pc_type mg -da_refine 10 -ksp_monitor -ksp_view -log_view
 21: */

 23: static char help[] = "Solves 2D Poisson equation using multigrid.\n\n";

 25:  #include <petscdm.h>
 26:  #include <petscdmda.h>
 27:  #include <petscksp.h>
 28:  #include <petscsys.h>
 29:  #include <petscvec.h>

 31: extern PetscErrorCode ComputeJacobian(KSP,Mat,Mat,void*);
 32: extern PetscErrorCode ComputeRHS(KSP,Vec,void*);

 34: typedef struct {
 35:   PetscScalar uu, tt;
 36: } UserContext;

 38: int main(int argc,char **argv)
 39: {
 40:   KSP            ksp;
 41:   DM             da;
 42:   UserContext    user;

 45:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 46:   KSPCreate(PETSC_COMM_WORLD,&ksp);
 47:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,11,11,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da);
 48:   DMSetFromOptions(da);
 49:   DMSetUp(da);
 50:   KSPSetDM(ksp,(DM)da);
 51:   DMSetApplicationContext(da,&user);

 53:   user.uu     = 1.0;
 54:   user.tt     = 1.0;

 56:   KSPSetComputeRHS(ksp,ComputeRHS,&user);
 57:   KSPSetComputeOperators(ksp,ComputeJacobian,&user);
 58:   KSPSetFromOptions(ksp);
 59:   KSPSolve(ksp,NULL,NULL);

 61:   DMDestroy(&da);
 62:   KSPDestroy(&ksp);
 63:   PetscFinalize();
 64:   return ierr;
 65: }

 67: PetscErrorCode ComputeRHS(KSP ksp,Vec b,void *ctx)
 68: {
 69:   UserContext    *user = (UserContext*)ctx;
 71:   PetscInt       i,j,M,N,xm,ym,xs,ys;
 72:   PetscScalar    Hx,Hy,pi,uu,tt;
 73:   PetscScalar    **array;
 74:   DM             da;
 75:   MatNullSpace   nullspace;

 78:   KSPGetDM(ksp,&da);
 79:   DMDAGetInfo(da, 0, &M, &N, 0,0,0,0,0,0,0,0,0,0);
 80:   uu   = user->uu; tt = user->tt;
 81:   pi   = 4*atan(1.0);
 82:   Hx   = 1.0/(PetscReal)(M);
 83:   Hy   = 1.0/(PetscReal)(N);

 85:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0); /* Fine grid */
 86:   DMDAVecGetArray(da, b, &array);
 87:   for (j=ys; j<ys+ym; j++) {
 88:     for (i=xs; i<xs+xm; i++) {
 89:       array[j][i] = -PetscCosScalar(uu*pi*((PetscReal)i+0.5)*Hx)*PetscCosScalar(tt*pi*((PetscReal)j+0.5)*Hy)*Hx*Hy;
 90:     }
 91:   }
 92:   DMDAVecRestoreArray(da, b, &array);
 93:   VecAssemblyBegin(b);
 94:   VecAssemblyEnd(b);

 96:   /* force right hand side to be consistent for singular matrix */
 97:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
 98:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
 99:   MatNullSpaceRemove(nullspace,b);
100:   MatNullSpaceDestroy(&nullspace);
101:   return(0);
102: }

104: PetscErrorCode ComputeJacobian(KSP ksp,Mat J, Mat jac,void *ctx)
105: {
107:   PetscInt       i, j, M, N, xm, ym, xs, ys, num, numi, numj;
108:   PetscScalar    v[5], Hx, Hy, HydHx, HxdHy;
109:   MatStencil     row, col[5];
110:   DM             da;
111:   MatNullSpace   nullspace;

114:   KSPGetDM(ksp,&da);
115:   DMDAGetInfo(da,0,&M,&N,0,0,0,0,0,0,0,0,0,0);
116:   Hx    = 1.0 / (PetscReal)(M);
117:   Hy    = 1.0 / (PetscReal)(N);
118:   HxdHy = Hx/Hy;
119:   HydHx = Hy/Hx;
120:   DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0);
121:   for (j=ys; j<ys+ym; j++) {
122:     for (i=xs; i<xs+xm; i++) {
123:       row.i = i; row.j = j;

125:       if (i==0 || j==0 || i==M-1 || j==N-1) {
126:         num=0; numi=0; numj=0;
127:         if (j!=0) {
128:           v[num] = -HxdHy;              col[num].i = i;   col[num].j = j-1;
129:           num++; numj++;
130:         }
131:         if (i!=0) {
132:           v[num] = -HydHx;              col[num].i = i-1; col[num].j = j;
133:           num++; numi++;
134:         }
135:         if (i!=M-1) {
136:           v[num] = -HydHx;              col[num].i = i+1; col[num].j = j;
137:           num++; numi++;
138:         }
139:         if (j!=N-1) {
140:           v[num] = -HxdHy;              col[num].i = i;   col[num].j = j+1;
141:           num++; numj++;
142:         }
143:         v[num] = ((PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx); col[num].i = i;   col[num].j = j;
144:         num++;
145:         MatSetValuesStencil(jac,1,&row,num,col,v,INSERT_VALUES);
146:       } else {
147:         v[0] = -HxdHy;              col[0].i = i;   col[0].j = j-1;
148:         v[1] = -HydHx;              col[1].i = i-1; col[1].j = j;
149:         v[2] = 2.0*(HxdHy + HydHx); col[2].i = i;   col[2].j = j;
150:         v[3] = -HydHx;              col[3].i = i+1; col[3].j = j;
151:         v[4] = -HxdHy;              col[4].i = i;   col[4].j = j+1;
152:         MatSetValuesStencil(jac,1,&row,5,col,v,INSERT_VALUES);
153:       }
154:     }
155:   }
156:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
157:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

159:   MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_TRUE,0,0,&nullspace);
160:   MatSetNullSpace(J,nullspace);
161:   MatNullSpaceDestroy(&nullspace);
162:   return(0);
163: }