Actual source code: ex5.c
petsc-3.14.5 2021-03-03
2: static char help[] = "Tests the multigrid code. The input parameters are:\n\
3: -x N Use a mesh in the x direction of N. \n\
4: -c N Use N V-cycles. \n\
5: -l N Use N Levels. \n\
6: -smooths N Use N pre smooths and N post smooths. \n\
7: -j Use Jacobi smoother. \n\
8: -a use additive multigrid \n\
9: -f use full multigrid (preconditioner variant) \n\
10: This example also demonstrates matrix-free methods\n\n";
12: /*
13: This is not a good example to understand the use of multigrid with PETSc.
14: */
16: #include <petscksp.h>
18: PetscErrorCode residual(Mat,Vec,Vec,Vec);
19: PetscErrorCode gauss_seidel(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscBool,PetscInt*,PCRichardsonConvergedReason*);
20: PetscErrorCode jacobi_smoother(PC,Vec,Vec,Vec,PetscReal,PetscReal,PetscReal,PetscInt,PetscBool,PetscInt*,PCRichardsonConvergedReason*);
21: PetscErrorCode interpolate(Mat,Vec,Vec,Vec);
22: PetscErrorCode restrct(Mat,Vec,Vec);
23: PetscErrorCode Create1dLaplacian(PetscInt,Mat*);
24: PetscErrorCode CalculateRhs(Vec);
25: PetscErrorCode CalculateError(Vec,Vec,Vec,PetscReal*);
26: PetscErrorCode CalculateSolution(PetscInt,Vec*);
27: PetscErrorCode amult(Mat,Vec,Vec);
28: PetscErrorCode apply_pc(PC,Vec,Vec);
30: int main(int Argc,char **Args)
31: {
32: PetscInt x_mesh = 15,levels = 3,cycles = 1,use_jacobi = 0;
33: PetscInt i,smooths = 1,*N,its;
35: PCMGType am = PC_MG_MULTIPLICATIVE;
36: Mat cmat,mat[20],fmat;
37: KSP cksp,ksp[20],kspmg;
38: PetscReal e[3]; /* l_2 error,max error, residual */
39: const char *shellname;
40: Vec x,solution,X[20],R[20],B[20];
41: PC pcmg,pc;
42: PetscBool flg;
44: PetscInitialize(&Argc,&Args,(char*)0,help);if (ierr) return ierr;
45: PetscOptionsGetInt(NULL,NULL,"-x",&x_mesh,NULL);
46: PetscOptionsGetInt(NULL,NULL,"-l",&levels,NULL);
47: PetscOptionsGetInt(NULL,NULL,"-c",&cycles,NULL);
48: PetscOptionsGetInt(NULL,NULL,"-smooths",&smooths,NULL);
49: PetscOptionsHasName(NULL,NULL,"-a",&flg);
51: if (flg) am = PC_MG_ADDITIVE;
52: PetscOptionsHasName(NULL,NULL,"-f",&flg);
53: if (flg) am = PC_MG_FULL;
54: PetscOptionsHasName(NULL,NULL,"-j",&flg);
55: if (flg) use_jacobi = 1;
57: PetscMalloc1(levels,&N);
58: N[0] = x_mesh;
59: for (i=1; i<levels; i++) {
60: N[i] = N[i-1]/2;
61: if (N[i] < 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Too many levels or N is not large enough");
62: }
64: Create1dLaplacian(N[levels-1],&cmat);
66: KSPCreate(PETSC_COMM_WORLD,&kspmg);
67: KSPGetPC(kspmg,&pcmg);
68: KSPSetFromOptions(kspmg);
69: PCSetType(pcmg,PCMG);
70: PCMGSetLevels(pcmg,levels,NULL);
71: PCMGSetType(pcmg,am);
73: PCMGGetCoarseSolve(pcmg,&cksp);
74: KSPSetOperators(cksp,cmat,cmat);
75: KSPGetPC(cksp,&pc);
76: PCSetType(pc,PCLU);
77: KSPSetType(cksp,KSPPREONLY);
79: /* zero is finest level */
80: for (i=0; i<levels-1; i++) {
81: Mat dummy;
83: PCMGSetResidual(pcmg,levels - 1 - i,residual,NULL);
84: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],NULL,&mat[i]);
85: MatShellSetOperation(mat[i],MATOP_MULT,(void (*)(void))restrct);
86: MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void (*)(void))interpolate);
87: PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);
88: PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);
89: PCMGSetCycleTypeOnLevel(pcmg,levels - 1 - i,(PCMGCycleType)cycles);
91: /* set smoother */
92: PCMGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);
93: KSPGetPC(ksp[i],&pc);
94: PCSetType(pc,PCSHELL);
95: PCShellSetName(pc,"user_precond");
96: PCShellGetName(pc,&shellname);
97: PetscPrintf(PETSC_COMM_WORLD,"level=%D, PCShell name is %s\n",i,shellname);
99: /* this is not used unless different options are passed to the solver */
100: MatCreateShell(PETSC_COMM_WORLD,N[i],N[i],N[i],N[i],NULL,&dummy);
101: MatShellSetOperation(dummy,MATOP_MULT,(void (*)(void))amult);
102: KSPSetOperators(ksp[i],dummy,dummy);
103: MatDestroy(&dummy);
105: /*
106: We override the matrix passed in by forcing it to use Richardson with
107: a user provided application. This is non-standard and this practice
108: should be avoided.
109: */
110: PCShellSetApply(pc,apply_pc);
111: PCShellSetApplyRichardson(pc,gauss_seidel);
112: if (use_jacobi) {
113: PCShellSetApplyRichardson(pc,jacobi_smoother);
114: }
115: KSPSetType(ksp[i],KSPRICHARDSON);
116: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
117: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
119: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
121: X[levels - 1 - i] = x;
122: if (i > 0) {
123: PCMGSetX(pcmg,levels - 1 - i,x);
124: }
125: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
127: B[levels -1 - i] = x;
128: if (i > 0) {
129: PCMGSetRhs(pcmg,levels - 1 - i,x);
130: }
131: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
133: R[levels - 1 - i] = x;
135: PCMGSetR(pcmg,levels - 1 - i,x);
136: }
137: /* create coarse level vectors */
138: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
139: PCMGSetX(pcmg,0,x); X[0] = x;
140: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
141: PCMGSetRhs(pcmg,0,x); B[0] = x;
143: /* create matrix multiply for finest level */
144: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],NULL,&fmat);
145: MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);
146: KSPSetOperators(kspmg,fmat,fmat);
148: CalculateSolution(N[0],&solution);
149: CalculateRhs(B[levels-1]);
150: VecSet(X[levels-1],0.0);
152: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
153: CalculateError(solution,X[levels-1],R[levels-1],e);
154: PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);
156: KSPSolve(kspmg,B[levels-1],X[levels-1]);
157: KSPGetIterationNumber(kspmg,&its);
158: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
159: CalculateError(solution,X[levels-1],R[levels-1],e);
160: PetscPrintf(PETSC_COMM_SELF,"its %D l_2 error %g max error %g resi %g\n",its,(double)e[0],(double)e[1],(double)e[2]);
162: PetscFree(N);
163: VecDestroy(&solution);
165: /* note we have to keep a list of all vectors allocated, this is
166: not ideal, but putting it in MGDestroy is not so good either*/
167: for (i=0; i<levels; i++) {
168: VecDestroy(&X[i]);
169: VecDestroy(&B[i]);
170: if (i) {VecDestroy(&R[i]);}
171: }
172: for (i=0; i<levels-1; i++) {
173: MatDestroy(&mat[i]);
174: }
175: MatDestroy(&cmat);
176: MatDestroy(&fmat);
177: KSPDestroy(&kspmg);
178: PetscFinalize();
179: return ierr;
180: }
182: /* --------------------------------------------------------------------- */
183: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
184: {
185: PetscInt i,n1;
186: PetscErrorCode ierr;
187: PetscScalar *x,*r;
188: const PetscScalar *b;
191: VecGetSize(bb,&n1);
192: VecGetArrayRead(bb,&b);
193: VecGetArray(xx,&x);
194: VecGetArray(rr,&r);
195: n1--;
196: r[0] = b[0] + x[1] - 2.0*x[0];
197: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
198: for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
199: VecRestoreArrayRead(bb,&b);
200: VecRestoreArray(xx,&x);
201: VecRestoreArray(rr,&r);
202: return(0);
203: }
205: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
206: {
207: PetscInt i,n1;
208: PetscErrorCode ierr;
209: PetscScalar *y;
210: const PetscScalar *x;
213: VecGetSize(xx,&n1);
214: VecGetArrayRead(xx,&x);
215: VecGetArray(yy,&y);
216: n1--;
217: y[0] = -x[1] + 2.0*x[0];
218: y[n1] = -x[n1-1] + 2.0*x[n1];
219: for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
220: VecRestoreArrayRead(xx,&x);
221: VecRestoreArray(yy,&y);
222: return(0);
223: }
225: /* --------------------------------------------------------------------- */
226: PetscErrorCode apply_pc(PC pc,Vec bb,Vec xx)
227: {
229: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Not implemented");
230: return(0);
231: }
233: 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)
234: {
235: PetscInt i,n1;
236: PetscErrorCode ierr;
237: PetscScalar *x;
238: const PetscScalar *b;
241: *its = m;
242: *reason = PCRICHARDSON_CONVERGED_ITS;
243: VecGetSize(bb,&n1); n1--;
244: VecGetArrayRead(bb,&b);
245: VecGetArray(xx,&x);
246: while (m--) {
247: x[0] = .5*(x[1] + b[0]);
248: for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
249: x[n1] = .5*(x[n1-1] + b[n1]);
250: for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
251: x[0] = .5*(x[1] + b[0]);
252: }
253: VecRestoreArrayRead(bb,&b);
254: VecRestoreArray(xx,&x);
255: return(0);
256: }
257: /* --------------------------------------------------------------------- */
258: 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)
259: {
260: PetscInt i,n,n1;
261: PetscErrorCode ierr;
262: PetscScalar *r,*x;
263: const PetscScalar *b;
266: *its = m;
267: *reason = PCRICHARDSON_CONVERGED_ITS;
268: VecGetSize(bb,&n); n1 = n - 1;
269: VecGetArrayRead(bb,&b);
270: VecGetArray(xx,&x);
271: VecGetArray(w,&r);
273: while (m--) {
274: r[0] = .5*(x[1] + b[0]);
275: for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]);
276: r[n1] = .5*(x[n1-1] + b[n1]);
277: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
278: }
279: VecRestoreArrayRead(bb,&b);
280: VecRestoreArray(xx,&x);
281: VecRestoreArray(w,&r);
282: return(0);
283: }
284: /*
285: We know for this application that yy and zz are the same
286: */
287: /* --------------------------------------------------------------------- */
288: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
289: {
290: PetscInt i,n,N,i2;
291: PetscErrorCode ierr;
292: PetscScalar *y;
293: const PetscScalar *x;
296: VecGetSize(yy,&N);
297: VecGetArrayRead(xx,&x);
298: VecGetArray(yy,&y);
299: n = N/2;
300: for (i=0; i<n; i++) {
301: i2 = 2*i;
302: y[i2] += .5*x[i];
303: y[i2+1] += x[i];
304: y[i2+2] += .5*x[i];
305: }
306: VecRestoreArrayRead(xx,&x);
307: VecRestoreArray(yy,&y);
308: return(0);
309: }
310: /* --------------------------------------------------------------------- */
311: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
312: {
313: PetscInt i,n,N,i2;
314: PetscErrorCode ierr;
315: PetscScalar *b;
316: const PetscScalar *r;
319: VecGetSize(rr,&N);
320: VecGetArrayRead(rr,&r);
321: VecGetArray(bb,&b);
322: n = N/2;
324: for (i=0; i<n; i++) {
325: i2 = 2*i;
326: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
327: }
328: VecRestoreArrayRead(rr,&r);
329: VecRestoreArray(bb,&b);
330: return(0);
331: }
332: /* --------------------------------------------------------------------- */
333: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
334: {
335: PetscScalar mone = -1.0,two = 2.0;
336: PetscInt i,idx;
340: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);
342: idx = n-1;
343: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
344: for (i=0; i<n-1; i++) {
345: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
346: idx = i+1;
347: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
348: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
349: }
350: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
351: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
352: return(0);
353: }
354: /* --------------------------------------------------------------------- */
355: PetscErrorCode CalculateRhs(Vec u)
356: {
358: PetscInt i,n;
359: PetscReal h,x = 0.0;
360: PetscScalar uu;
363: VecGetSize(u,&n);
364: h = 1.0/((PetscReal)(n+1));
365: for (i=0; i<n; i++) {
366: x += h; uu = 2.0*h*h;
367: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
368: }
369: return(0);
370: }
371: /* --------------------------------------------------------------------- */
372: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
373: {
375: PetscInt i;
376: PetscReal h,x = 0.0;
377: PetscScalar uu;
380: VecCreateSeq(PETSC_COMM_SELF,n,solution);
381: h = 1.0/((PetscReal)(n+1));
382: for (i=0; i<n; i++) {
383: x += h; uu = x*(1.-x);
384: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
385: }
386: return(0);
387: }
388: /* --------------------------------------------------------------------- */
389: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
390: {
394: VecNorm(r,NORM_2,e+2);
395: VecWAXPY(r,-1.0,u,solution);
396: VecNorm(r,NORM_2,e);
397: VecNorm(r,NORM_1,e+1);
398: return(0);
399: }
401: /*TEST
403: test:
405: TEST*/