Actual source code: ex5.c

  1: static char help[] = "Tests the multigrid code.  The input parameters are:\n\
  2:   -x N              Use a mesh in the x direction of N.  \n\
  3:   -c N              Use N V-cycles.  \n\
  4:   -l N              Use N Levels.  \n\
  5:   -smooths N        Use N pre smooths and N post smooths.  \n\
  6:   -j                Use Jacobi smoother.  \n\
  7:   -a use additive multigrid \n\
  8:   -f use full multigrid (preconditioner variant) \n\
  9: This example also demonstrates matrix-free methods\n\n";

 11: /*
 12:   This is not a good example to understand the use of multigrid with PETSc.
 13: */

 15: #include <petscksp.h>

 17: PetscErrorCode residual(Mat, Vec, Vec, Vec);
 18: PetscErrorCode gauss_seidel(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *);
 19: PetscErrorCode jacobi_smoother(PC, Vec, Vec, Vec, PetscReal, PetscReal, PetscReal, PetscInt, PetscBool, PetscInt *, PCRichardsonConvergedReason *);
 20: PetscErrorCode interpolate(Mat, Vec, Vec, Vec);
 21: PetscErrorCode restrct(Mat, Vec, Vec);
 22: PetscErrorCode Create1dLaplacian(PetscInt, Mat *);
 23: PetscErrorCode CalculateRhs(Vec);
 24: PetscErrorCode CalculateError(Vec, Vec, Vec, PetscReal *);
 25: PetscErrorCode CalculateSolution(PetscInt, Vec *);
 26: PetscErrorCode amult(Mat, Vec, Vec);
 27: PetscErrorCode apply_pc(PC, Vec, Vec);

 29: int main(int Argc, char **Args)
 30: {
 31:   PetscInt    x_mesh = 15, levels = 3, cycles = 1, use_jacobi = 0;
 32:   PetscInt    i, smooths = 1, *N, its;
 33:   PCMGType    am = PC_MG_MULTIPLICATIVE;
 34:   Mat         cmat, mat[20], fmat;
 35:   KSP         cksp, ksp[20], kspmg;
 36:   PetscReal   e[3]; /* l_2 error,max error, residual */
 37:   const char *shellname;
 38:   Vec         x, solution, X[20], R[20], B[20];
 39:   PC          pcmg, pc;
 40:   PetscBool   flg;

 42:   PetscCall(PetscInitialize(&Argc, &Args, NULL, help));
 43:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-x", &x_mesh, NULL));
 44:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-l", &levels, NULL));
 45:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-c", &cycles, NULL));
 46:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-smooths", &smooths, NULL));
 47:   PetscCall(PetscOptionsHasName(NULL, NULL, "-a", &flg));

 49:   if (flg) am = PC_MG_ADDITIVE;
 50:   PetscCall(PetscOptionsHasName(NULL, NULL, "-f", &flg));
 51:   if (flg) am = PC_MG_FULL;
 52:   PetscCall(PetscOptionsHasName(NULL, NULL, "-j", &flg));
 53:   if (flg) use_jacobi = 1;

 55:   PetscCall(PetscMalloc1(levels, &N));
 56:   N[0] = x_mesh;
 57:   for (i = 1; i < levels; i++) {
 58:     N[i] = N[i - 1] / 2;
 59:     PetscCheck(N[i] >= 1, PETSC_COMM_WORLD, PETSC_ERR_USER, "Too many levels or N is not large enough");
 60:   }

 62:   PetscCall(Create1dLaplacian(N[levels - 1], &cmat));

 64:   PetscCall(KSPCreate(PETSC_COMM_WORLD, &kspmg));
 65:   PetscCall(KSPGetPC(kspmg, &pcmg));
 66:   PetscCall(KSPSetFromOptions(kspmg));
 67:   PetscCall(PCSetType(pcmg, PCMG));
 68:   PetscCall(PCMGSetLevels(pcmg, levels, NULL));
 69:   PetscCall(PCMGSetType(pcmg, am));

 71:   PetscCall(PCMGGetCoarseSolve(pcmg, &cksp));
 72:   PetscCall(KSPSetOperators(cksp, cmat, cmat));
 73:   PetscCall(KSPGetPC(cksp, &pc));
 74:   PetscCall(PCSetType(pc, PCLU));
 75:   PetscCall(KSPSetType(cksp, KSPPREONLY));

 77:   /* zero is finest level */
 78:   for (i = 0; i < levels - 1; i++) {
 79:     Mat dummy;

 81:     PetscCall(PCMGSetResidual(pcmg, levels - 1 - i, residual, NULL));
 82:     PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i + 1], N[i], N[i + 1], N[i], NULL, &mat[i]));
 83:     PetscCall(MatShellSetOperation(mat[i], MATOP_MULT, (void (*)(void))restrct));
 84:     PetscCall(MatShellSetOperation(mat[i], MATOP_MULT_TRANSPOSE_ADD, (void (*)(void))interpolate));
 85:     PetscCall(PCMGSetInterpolation(pcmg, levels - 1 - i, mat[i]));
 86:     PetscCall(PCMGSetRestriction(pcmg, levels - 1 - i, mat[i]));
 87:     PetscCall(PCMGSetCycleTypeOnLevel(pcmg, levels - 1 - i, (PCMGCycleType)cycles));

 89:     /* set smoother */
 90:     PetscCall(PCMGGetSmoother(pcmg, levels - 1 - i, &ksp[i]));
 91:     PetscCall(KSPGetPC(ksp[i], &pc));
 92:     PetscCall(PCSetType(pc, PCSHELL));
 93:     PetscCall(PCShellSetName(pc, "user_precond"));
 94:     PetscCall(PCShellGetName(pc, &shellname));
 95:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "level=%" PetscInt_FMT ", PCShell name is %s\n", i, shellname));

 97:     /* this is not used unless different options are passed to the solver */
 98:     PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[i], N[i], N[i], N[i], NULL, &dummy));
 99:     PetscCall(MatShellSetOperation(dummy, MATOP_MULT, (void (*)(void))amult));
100:     PetscCall(KSPSetOperators(ksp[i], dummy, dummy));
101:     PetscCall(MatDestroy(&dummy));

103:     /*
104:         We override the matrix passed in by forcing it to use Richardson with
105:         a user provided application. This is non-standard and this practice
106:         should be avoided.
107:     */
108:     PetscCall(PCShellSetApply(pc, apply_pc));
109:     PetscCall(PCShellSetApplyRichardson(pc, gauss_seidel));
110:     if (use_jacobi) PetscCall(PCShellSetApplyRichardson(pc, jacobi_smoother));
111:     PetscCall(KSPSetType(ksp[i], KSPRICHARDSON));
112:     PetscCall(KSPSetInitialGuessNonzero(ksp[i], PETSC_TRUE));
113:     PetscCall(KSPSetTolerances(ksp[i], PETSC_CURRENT, PETSC_CURRENT, PETSC_CURRENT, smooths));

115:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));

117:     X[levels - 1 - i] = x;
118:     if (i > 0) PetscCall(PCMGSetX(pcmg, levels - 1 - i, x));
119:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));

121:     B[levels - 1 - i] = x;
122:     if (i > 0) PetscCall(PCMGSetRhs(pcmg, levels - 1 - i, x));
123:     PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[i], &x));

125:     R[levels - 1 - i] = x;

127:     PetscCall(PCMGSetR(pcmg, levels - 1 - i, x));
128:   }
129:   /* create coarse level vectors */
130:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x));
131:   PetscCall(PCMGSetX(pcmg, 0, x));
132:   X[0] = x;
133:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, N[levels - 1], &x));
134:   PetscCall(PCMGSetRhs(pcmg, 0, x));
135:   B[0] = x;

137:   /* create matrix multiply for finest level */
138:   PetscCall(MatCreateShell(PETSC_COMM_WORLD, N[0], N[0], N[0], N[0], NULL, &fmat));
139:   PetscCall(MatShellSetOperation(fmat, MATOP_MULT, (void (*)(void))amult));
140:   PetscCall(KSPSetOperators(kspmg, fmat, fmat));

142:   PetscCall(CalculateSolution(N[0], &solution));
143:   PetscCall(CalculateRhs(B[levels - 1]));
144:   PetscCall(VecSet(X[levels - 1], 0.0));

146:   PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1]));
147:   PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e));
148:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "l_2 error %g max error %g resi %g\n", (double)e[0], (double)e[1], (double)e[2]));

150:   PetscCall(KSPSolve(kspmg, B[levels - 1], X[levels - 1]));
151:   PetscCall(KSPGetIterationNumber(kspmg, &its));
152:   PetscCall(residual((Mat)0, B[levels - 1], X[levels - 1], R[levels - 1]));
153:   PetscCall(CalculateError(solution, X[levels - 1], R[levels - 1], e));
154:   PetscCall(PetscPrintf(PETSC_COMM_SELF, "its %" PetscInt_FMT " l_2 error %g max error %g resi %g\n", its, (double)e[0], (double)e[1], (double)e[2]));

156:   PetscCall(PetscFree(N));
157:   PetscCall(VecDestroy(&solution));

159:   /* note we have to keep a list of all vectors allocated, this is
160:      not ideal, but putting it in MGDestroy is not so good either*/
161:   for (i = 0; i < levels; i++) {
162:     PetscCall(VecDestroy(&X[i]));
163:     PetscCall(VecDestroy(&B[i]));
164:     if (i) PetscCall(VecDestroy(&R[i]));
165:   }
166:   for (i = 0; i < levels - 1; i++) PetscCall(MatDestroy(&mat[i]));
167:   PetscCall(MatDestroy(&cmat));
168:   PetscCall(MatDestroy(&fmat));
169:   PetscCall(KSPDestroy(&kspmg));
170:   PetscCall(PetscFinalize());
171:   return 0;
172: }

174: PetscErrorCode residual(Mat mat, Vec bb, Vec xx, Vec rr)
175: {
176:   PetscInt           i, n1;
177:   PetscScalar       *x, *r;
178:   const PetscScalar *b;

180:   PetscFunctionBegin;
181:   PetscCall(VecGetSize(bb, &n1));
182:   PetscCall(VecGetArrayRead(bb, &b));
183:   PetscCall(VecGetArray(xx, &x));
184:   PetscCall(VecGetArray(rr, &r));
185:   n1--;
186:   r[0]  = b[0] + x[1] - 2.0 * x[0];
187:   r[n1] = b[n1] + x[n1 - 1] - 2.0 * x[n1];
188:   for (i = 1; i < n1; i++) r[i] = b[i] + x[i + 1] + x[i - 1] - 2.0 * x[i];
189:   PetscCall(VecRestoreArrayRead(bb, &b));
190:   PetscCall(VecRestoreArray(xx, &x));
191:   PetscCall(VecRestoreArray(rr, &r));
192:   PetscFunctionReturn(PETSC_SUCCESS);
193: }

195: PetscErrorCode amult(Mat mat, Vec xx, Vec yy)
196: {
197:   PetscInt           i, n1;
198:   PetscScalar       *y;
199:   const PetscScalar *x;

201:   PetscFunctionBegin;
202:   PetscCall(VecGetSize(xx, &n1));
203:   PetscCall(VecGetArrayRead(xx, &x));
204:   PetscCall(VecGetArray(yy, &y));
205:   n1--;
206:   y[0]  = -x[1] + 2.0 * x[0];
207:   y[n1] = -x[n1 - 1] + 2.0 * x[n1];
208:   for (i = 1; i < n1; i++) y[i] = -x[i + 1] - x[i - 1] + 2.0 * x[i];
209:   PetscCall(VecRestoreArrayRead(xx, &x));
210:   PetscCall(VecRestoreArray(yy, &y));
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }

214: PetscErrorCode apply_pc(PC pc, Vec bb, Vec xx)
215: {
216:   PetscFunctionBegin;
217:   SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Not implemented");
218: }

220: PetscErrorCode gauss_seidel(PC pc, Vec bb, Vec xx, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt m, PetscBool guesszero, PetscInt *its, PCRichardsonConvergedReason *reason)
221: {
222:   PetscInt           i, n1;
223:   PetscScalar       *x;
224:   const PetscScalar *b;

226:   PetscFunctionBegin;
227:   *its    = m;
228:   *reason = PCRICHARDSON_CONVERGED_ITS;
229:   PetscCall(VecGetSize(bb, &n1));
230:   n1--;
231:   PetscCall(VecGetArrayRead(bb, &b));
232:   PetscCall(VecGetArray(xx, &x));
233:   while (m--) {
234:     x[0] = .5 * (x[1] + b[0]);
235:     for (i = 1; i < n1; i++) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
236:     x[n1] = .5 * (x[n1 - 1] + b[n1]);
237:     for (i = n1 - 1; i > 0; i--) x[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
238:     x[0] = .5 * (x[1] + b[0]);
239:   }
240:   PetscCall(VecRestoreArrayRead(bb, &b));
241:   PetscCall(VecRestoreArray(xx, &x));
242:   PetscFunctionReturn(PETSC_SUCCESS);
243: }

245: PetscErrorCode jacobi_smoother(PC pc, Vec bb, Vec xx, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt m, PetscBool guesszero, PetscInt *its, PCRichardsonConvergedReason *reason)
246: {
247:   PetscInt           i, n, n1;
248:   PetscScalar       *r, *x;
249:   const PetscScalar *b;

251:   PetscFunctionBegin;
252:   *its    = m;
253:   *reason = PCRICHARDSON_CONVERGED_ITS;
254:   PetscCall(VecGetSize(bb, &n));
255:   n1 = n - 1;
256:   PetscCall(VecGetArrayRead(bb, &b));
257:   PetscCall(VecGetArray(xx, &x));
258:   PetscCall(VecGetArray(w, &r));

260:   while (m--) {
261:     r[0] = .5 * (x[1] + b[0]);
262:     for (i = 1; i < n1; i++) r[i] = .5 * (x[i + 1] + x[i - 1] + b[i]);
263:     r[n1] = .5 * (x[n1 - 1] + b[n1]);
264:     for (i = 0; i < n; i++) x[i] = (2.0 * r[i] + x[i]) / 3.0;
265:   }
266:   PetscCall(VecRestoreArrayRead(bb, &b));
267:   PetscCall(VecRestoreArray(xx, &x));
268:   PetscCall(VecRestoreArray(w, &r));
269:   PetscFunctionReturn(PETSC_SUCCESS);
270: }
271: /*
272:    We know for this application that yy  and zz are the same
273: */

275: PetscErrorCode interpolate(Mat mat, Vec xx, Vec yy, Vec zz)
276: {
277:   PetscInt           i, n, N, i2;
278:   PetscScalar       *y;
279:   const PetscScalar *x;

281:   PetscFunctionBegin;
282:   PetscCall(VecGetSize(yy, &N));
283:   PetscCall(VecGetArrayRead(xx, &x));
284:   PetscCall(VecGetArray(yy, &y));
285:   n = N / 2;
286:   for (i = 0; i < n; i++) {
287:     i2 = 2 * i;
288:     y[i2] += .5 * x[i];
289:     y[i2 + 1] += x[i];
290:     y[i2 + 2] += .5 * x[i];
291:   }
292:   PetscCall(VecRestoreArrayRead(xx, &x));
293:   PetscCall(VecRestoreArray(yy, &y));
294:   PetscFunctionReturn(PETSC_SUCCESS);
295: }

297: PetscErrorCode restrct(Mat mat, Vec rr, Vec bb)
298: {
299:   PetscInt           i, n, N, i2;
300:   PetscScalar       *b;
301:   const PetscScalar *r;

303:   PetscFunctionBegin;
304:   PetscCall(VecGetSize(rr, &N));
305:   PetscCall(VecGetArrayRead(rr, &r));
306:   PetscCall(VecGetArray(bb, &b));
307:   n = N / 2;

309:   for (i = 0; i < n; i++) {
310:     i2   = 2 * i;
311:     b[i] = (r[i2] + 2.0 * r[i2 + 1] + r[i2 + 2]);
312:   }
313:   PetscCall(VecRestoreArrayRead(rr, &r));
314:   PetscCall(VecRestoreArray(bb, &b));
315:   PetscFunctionReturn(PETSC_SUCCESS);
316: }

318: PetscErrorCode Create1dLaplacian(PetscInt n, Mat *mat)
319: {
320:   PetscScalar mone = -1.0, two = 2.0;
321:   PetscInt    i, idx;

323:   PetscFunctionBegin;
324:   PetscCall(MatCreateSeqAIJ(PETSC_COMM_SELF, n, n, 3, NULL, mat));

326:   idx = n - 1;
327:   PetscCall(MatSetValues(*mat, 1, &idx, 1, &idx, &two, INSERT_VALUES));
328:   for (i = 0; i < n - 1; i++) {
329:     PetscCall(MatSetValues(*mat, 1, &i, 1, &i, &two, INSERT_VALUES));
330:     idx = i + 1;
331:     PetscCall(MatSetValues(*mat, 1, &idx, 1, &i, &mone, INSERT_VALUES));
332:     PetscCall(MatSetValues(*mat, 1, &i, 1, &idx, &mone, INSERT_VALUES));
333:   }
334:   PetscCall(MatAssemblyBegin(*mat, MAT_FINAL_ASSEMBLY));
335:   PetscCall(MatAssemblyEnd(*mat, MAT_FINAL_ASSEMBLY));
336:   PetscFunctionReturn(PETSC_SUCCESS);
337: }

339: PetscErrorCode CalculateRhs(Vec u)
340: {
341:   PetscInt    i, n;
342:   PetscReal   h;
343:   PetscScalar uu;

345:   PetscFunctionBegin;
346:   PetscCall(VecGetSize(u, &n));
347:   h = 1.0 / ((PetscReal)(n + 1));
348:   for (i = 0; i < n; i++) {
349:     uu = 2.0 * h * h;
350:     PetscCall(VecSetValues(u, 1, &i, &uu, INSERT_VALUES));
351:   }
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: PetscErrorCode CalculateSolution(PetscInt n, Vec *solution)
356: {
357:   PetscInt    i;
358:   PetscReal   h, x = 0.0;
359:   PetscScalar uu;

361:   PetscFunctionBegin;
362:   PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, solution));
363:   h = 1.0 / ((PetscReal)(n + 1));
364:   for (i = 0; i < n; i++) {
365:     x += h;
366:     uu = x * (1. - x);
367:     PetscCall(VecSetValues(*solution, 1, &i, &uu, INSERT_VALUES));
368:   }
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: PetscErrorCode CalculateError(Vec solution, Vec u, Vec r, PetscReal *e)
373: {
374:   PetscFunctionBegin;
375:   PetscCall(VecNorm(r, NORM_2, e + 2));
376:   PetscCall(VecWAXPY(r, -1.0, u, solution));
377:   PetscCall(VecNorm(r, NORM_2, e));
378:   PetscCall(VecNorm(r, NORM_1, e + 1));
379:   PetscFunctionReturn(PETSC_SUCCESS);
380: }

382: /*TEST

384:    test:

386: TEST*/