Actual source code: ex32.c


  2: /*
  3: Laplacian in 2D. Modeled by the partial differential equation

  5:    div  grad u = f,  0 < x,y < 1,

  7: with forcing function

  9:    f = e^{-(1 - x)^2/\nu} e^{-(1 - y)^2/\nu}

 11: with pure Neumann boundary conditions

 13: The functions are cell-centered

 15: This uses multigrid to solve the linear system

 17:        Contributed by Andrei Draganescu <aidraga@sandia.gov>

 19: Note the nice multigrid convergence despite the fact it is only using
 20: piecewise constant interpolation/restriction. This is because cell-centered multigrid
 21: does not need the same rule:

 23:     polynomial degree(interpolation) + polynomial degree(restriction) + 2 > degree of PDE

 25: that vertex based multigrid needs.
 26: */

 28: static char help[] = "Solves 2D inhomogeneous Laplacian using multigrid.\n\n";

 30: #include <petscdm.h>
 31: #include <petscdmda.h>
 32: #include <petscksp.h>

 34: extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *);
 35: extern PetscErrorCode ComputeRHS(KSP, Vec, void *);

 37: typedef enum {
 38:   DIRICHLET,
 39:   NEUMANN
 40: } BCType;

 42: typedef struct {
 43:   PetscScalar nu;
 44:   BCType      bcType;
 45: } UserContext;

 47: int main(int argc, char **argv)
 48: {
 49:   KSP         ksp;
 50:   DM          da;
 51:   UserContext user;
 52:   const char *bcTypes[2] = {"dirichlet", "neumann"};
 53:   PetscInt    bc;

 56:   PetscInitialize(&argc, &argv, (char *)0, help);
 57:   KSPCreate(PETSC_COMM_WORLD, &ksp);
 58:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 12, 12, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, &da);
 59:   DMSetFromOptions(da);
 60:   DMSetUp(da);
 61:   DMDASetInterpolationType(da, DMDA_Q0);

 63:   KSPSetDM(ksp, da);

 65:   PetscOptionsBegin(PETSC_COMM_WORLD, "", "Options for the inhomogeneous Poisson equation", "DM");
 66:   user.nu = 0.1;
 67:   PetscOptionsScalar("-nu", "The width of the Gaussian source", "ex29.c", 0.1, &user.nu, NULL);
 68:   bc = (PetscInt)NEUMANN;
 69:   PetscOptionsEList("-bc_type", "Type of boundary condition", "ex29.c", bcTypes, 2, bcTypes[0], &bc, NULL);
 70:   user.bcType = (BCType)bc;
 71:   PetscOptionsEnd();

 73:   KSPSetComputeRHS(ksp, ComputeRHS, &user);
 74:   KSPSetComputeOperators(ksp, ComputeMatrix, &user);
 75:   KSPSetFromOptions(ksp);
 76:   KSPSolve(ksp, NULL, NULL);
 77:   KSPDestroy(&ksp);
 78:   DMDestroy(&da);
 79:   PetscFinalize();
 80:   return 0;
 81: }

 83: PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx)
 84: {
 85:   UserContext  *user = (UserContext *)ctx;
 86:   PetscInt      i, j, mx, my, xm, ym, xs, ys;
 87:   PetscScalar   Hx, Hy;
 88:   PetscScalar **array;
 89:   DM            da;

 92:   KSPGetDM(ksp, &da);
 93:   DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
 94:   Hx = 1.0 / (PetscReal)(mx);
 95:   Hy = 1.0 / (PetscReal)(my);
 96:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
 97:   DMDAVecGetArray(da, b, &array);
 98:   for (j = ys; j < ys + ym; j++) {
 99:     for (i = xs; i < xs + xm; i++) array[j][i] = PetscExpScalar(-(((PetscReal)i + 0.5) * Hx) * (((PetscReal)i + 0.5) * Hx) / user->nu) * PetscExpScalar(-(((PetscReal)j + 0.5) * Hy) * (((PetscReal)j + 0.5) * Hy) / user->nu) * Hx * Hy;
100:   }
101:   DMDAVecRestoreArray(da, b, &array);
102:   VecAssemblyBegin(b);
103:   VecAssemblyEnd(b);

105:   /* force right hand side to be consistent for singular matrix */
106:   /* note this is really a hack, normally the model would provide you with a consistent right handside */
107:   if (user->bcType == NEUMANN) {
108:     MatNullSpace nullspace;

110:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
111:     MatNullSpaceRemove(nullspace, b);
112:     MatNullSpaceDestroy(&nullspace);
113:   }
114:   return 0;
115: }

117: PetscErrorCode ComputeMatrix(KSP ksp, Mat J, Mat jac, void *ctx)
118: {
119:   UserContext *user = (UserContext *)ctx;
120:   PetscInt     i, j, mx, my, xm, ym, xs, ys, num, numi, numj;
121:   PetscScalar  v[5], Hx, Hy, HydHx, HxdHy;
122:   MatStencil   row, col[5];
123:   DM           da;

126:   KSPGetDM(ksp, &da);
127:   DMDAGetInfo(da, 0, &mx, &my, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
128:   Hx    = 1.0 / (PetscReal)(mx);
129:   Hy    = 1.0 / (PetscReal)(my);
130:   HxdHy = Hx / Hy;
131:   HydHx = Hy / Hx;
132:   DMDAGetCorners(da, &xs, &ys, 0, &xm, &ym, 0);
133:   for (j = ys; j < ys + ym; j++) {
134:     for (i = xs; i < xs + xm; i++) {
135:       row.i = i;
136:       row.j = j;
137:       if (i == 0 || j == 0 || i == mx - 1 || j == my - 1) {
138:         if (user->bcType == DIRICHLET) {
139:           v[0] = 2.0 * (HxdHy + HydHx);
140:           MatSetValuesStencil(jac, 1, &row, 1, &row, v, INSERT_VALUES);
141:           SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Dirichlet boundary conditions not supported !");
142:         } else if (user->bcType == NEUMANN) {
143:           num  = 0;
144:           numi = 0;
145:           numj = 0;
146:           if (j != 0) {
147:             v[num]     = -HxdHy;
148:             col[num].i = i;
149:             col[num].j = j - 1;
150:             num++;
151:             numj++;
152:           }
153:           if (i != 0) {
154:             v[num]     = -HydHx;
155:             col[num].i = i - 1;
156:             col[num].j = j;
157:             num++;
158:             numi++;
159:           }
160:           if (i != mx - 1) {
161:             v[num]     = -HydHx;
162:             col[num].i = i + 1;
163:             col[num].j = j;
164:             num++;
165:             numi++;
166:           }
167:           if (j != my - 1) {
168:             v[num]     = -HxdHy;
169:             col[num].i = i;
170:             col[num].j = j + 1;
171:             num++;
172:             numj++;
173:           }
174:           v[num]     = (PetscReal)(numj)*HxdHy + (PetscReal)(numi)*HydHx;
175:           col[num].i = i;
176:           col[num].j = j;
177:           num++;
178:           MatSetValuesStencil(jac, 1, &row, num, col, v, INSERT_VALUES);
179:         }
180:       } else {
181:         v[0]     = -HxdHy;
182:         col[0].i = i;
183:         col[0].j = j - 1;
184:         v[1]     = -HydHx;
185:         col[1].i = i - 1;
186:         col[1].j = j;
187:         v[2]     = 2.0 * (HxdHy + HydHx);
188:         col[2].i = i;
189:         col[2].j = j;
190:         v[3]     = -HydHx;
191:         col[3].i = i + 1;
192:         col[3].j = j;
193:         v[4]     = -HxdHy;
194:         col[4].i = i;
195:         col[4].j = j + 1;
196:         MatSetValuesStencil(jac, 1, &row, 5, col, v, INSERT_VALUES);
197:       }
198:     }
199:   }
200:   MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
201:   MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
202:   if (user->bcType == NEUMANN) {
203:     MatNullSpace nullspace;

205:     MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, 0, &nullspace);
206:     MatSetNullSpace(J, nullspace);
207:     MatNullSpaceDestroy(&nullspace);
208:   }
209:   return 0;
210: }

212: /*TEST

214:    test:
215:       args: -pc_type mg -pc_mg_type full -ksp_type fgmres -ksp_monitor_short -pc_mg_levels 3 -mg_coarse_pc_factor_shift_type nonzero

217: TEST*/