Actual source code: ex5.c
petsc-3.11.4 2019-09-28
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);
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;
34: PCMGType am = PC_MG_MULTIPLICATIVE;
35: Mat cmat,mat[20],fmat;
36: KSP cksp,ksp[20],kspmg;
37: PetscReal e[3]; /* l_2 error,max error, residual */
38: const char *shellname;
39: Vec x,solution,X[20],R[20],B[20];
40: PC pcmg,pc;
41: PetscBool flg;
43: PetscInitialize(&Argc,&Args,(char*)0,help);if (ierr) return ierr;
44: PetscOptionsGetInt(NULL,NULL,"-x",&x_mesh,NULL);
45: PetscOptionsGetInt(NULL,NULL,"-l",&levels,NULL);
46: PetscOptionsGetInt(NULL,NULL,"-c",&cycles,NULL);
47: PetscOptionsGetInt(NULL,NULL,"-smooths",&smooths,NULL);
48: PetscOptionsHasName(NULL,NULL,"-a",&flg);
50: if (flg) am = PC_MG_ADDITIVE;
51: PetscOptionsHasName(NULL,NULL,"-f",&flg);
52: if (flg) am = PC_MG_FULL;
53: PetscOptionsHasName(NULL,NULL,"-j",&flg);
54: if (flg) use_jacobi = 1;
56: PetscMalloc1(levels,&N);
57: N[0] = x_mesh;
58: for (i=1; i<levels; i++) {
59: N[i] = N[i-1]/2;
60: if (N[i] < 1) SETERRQ(PETSC_COMM_WORLD,1,"Too many levels");
61: }
63: Create1dLaplacian(N[levels-1],&cmat);
65: KSPCreate(PETSC_COMM_WORLD,&kspmg);
66: KSPGetPC(kspmg,&pcmg);
67: KSPSetFromOptions(kspmg);
68: PCSetType(pcmg,PCMG);
69: PCMGSetLevels(pcmg,levels,NULL);
70: PCMGSetType(pcmg,am);
72: PCMGGetCoarseSolve(pcmg,&cksp);
73: KSPSetOperators(cksp,cmat,cmat);
74: KSPGetPC(cksp,&pc);
75: PCSetType(pc,PCLU);
76: KSPSetType(cksp,KSPPREONLY);
78: /* zero is finest level */
79: for (i=0; i<levels-1; i++) {
80: PCMGSetResidual(pcmg,levels - 1 - i,residual,(Mat)0);
81: MatCreateShell(PETSC_COMM_WORLD,N[i+1],N[i],N[i+1],N[i],(void*)0,&mat[i]);
82: MatShellSetOperation(mat[i],MATOP_MULT,(void (*)(void))restrct);
83: MatShellSetOperation(mat[i],MATOP_MULT_TRANSPOSE_ADD,(void (*)(void))interpolate);
84: PCMGSetInterpolation(pcmg,levels - 1 - i,mat[i]);
85: PCMGSetRestriction(pcmg,levels - 1 - i,mat[i]);
86: PCMGSetCycleTypeOnLevel(pcmg,levels - 1 - i,(PCMGCycleType)cycles);
88: /* set smoother */
89: PCMGGetSmoother(pcmg,levels - 1 - i,&ksp[i]);
90: KSPGetPC(ksp[i],&pc);
91: PCSetType(pc,PCSHELL);
92: PCShellSetName(pc,"user_precond");
93: PCShellGetName(pc,&shellname);
94: PetscPrintf(PETSC_COMM_WORLD,"level=%D, PCShell name is %s\n",i,shellname);
96: /* this is a dummy! since KSP requires a matrix passed in */
97: KSPSetOperators(ksp[i],mat[i],mat[i]);
98: /*
99: We override the matrix passed in by forcing it to use Richardson with
100: a user provided application. This is non-standard and this practice
101: should be avoided.
102: */
103: PCShellSetApplyRichardson(pc,gauss_seidel);
104: if (use_jacobi) {
105: PCShellSetApplyRichardson(pc,jacobi_smoother);
106: }
107: KSPSetType(ksp[i],KSPRICHARDSON);
108: KSPSetInitialGuessNonzero(ksp[i],PETSC_TRUE);
109: KSPSetTolerances(ksp[i],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,smooths);
111: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
113: X[levels - 1 - i] = x;
114: if (i > 0) {
115: PCMGSetX(pcmg,levels - 1 - i,x);
116: }
117: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
119: B[levels -1 - i] = x;
120: if (i > 0) {
121: PCMGSetRhs(pcmg,levels - 1 - i,x);
122: }
123: VecCreateSeq(PETSC_COMM_SELF,N[i],&x);
125: R[levels - 1 - i] = x;
127: PCMGSetR(pcmg,levels - 1 - i,x);
128: }
129: /* create coarse level vectors */
130: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
131: PCMGSetX(pcmg,0,x); X[0] = x;
132: VecCreateSeq(PETSC_COMM_SELF,N[levels-1],&x);
133: PCMGSetRhs(pcmg,0,x); B[0] = x;
135: /* create matrix multiply for finest level */
136: MatCreateShell(PETSC_COMM_WORLD,N[0],N[0],N[0],N[0],(void*)0,&fmat);
137: MatShellSetOperation(fmat,MATOP_MULT,(void (*)(void))amult);
138: KSPSetOperators(kspmg,fmat,fmat);
140: CalculateSolution(N[0],&solution);
141: CalculateRhs(B[levels-1]);
142: VecSet(X[levels-1],0.0);
144: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
145: CalculateError(solution,X[levels-1],R[levels-1],e);
146: PetscPrintf(PETSC_COMM_SELF,"l_2 error %g max error %g resi %g\n",(double)e[0],(double)e[1],(double)e[2]);
148: KSPSolve(kspmg,B[levels-1],X[levels-1]);
149: KSPGetIterationNumber(kspmg,&its);
150: residual((Mat)0,B[levels-1],X[levels-1],R[levels-1]);
151: CalculateError(solution,X[levels-1],R[levels-1],e);
152: 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]);
154: PetscFree(N);
155: VecDestroy(&solution);
157: /* note we have to keep a list of all vectors allocated, this is
158: not ideal, but putting it in MGDestroy is not so good either*/
159: for (i=0; i<levels; i++) {
160: VecDestroy(&X[i]);
161: VecDestroy(&B[i]);
162: if (i) {VecDestroy(&R[i]);}
163: }
164: for (i=0; i<levels-1; i++) {
165: MatDestroy(&mat[i]);
166: }
167: MatDestroy(&cmat);
168: MatDestroy(&fmat);
169: KSPDestroy(&kspmg);
170: PetscFinalize();
171: return ierr;
172: }
174: /* --------------------------------------------------------------------- */
175: PetscErrorCode residual(Mat mat,Vec bb,Vec xx,Vec rr)
176: {
177: PetscInt i,n1;
178: PetscErrorCode ierr;
179: PetscScalar *x,*r;
180: const PetscScalar *b;
183: VecGetSize(bb,&n1);
184: VecGetArrayRead(bb,&b);
185: VecGetArray(xx,&x);
186: VecGetArray(rr,&r);
187: n1--;
188: r[0] = b[0] + x[1] - 2.0*x[0];
189: r[n1] = b[n1] + x[n1-1] - 2.0*x[n1];
190: for (i=1; i<n1; i++) r[i] = b[i] + x[i+1] + x[i-1] - 2.0*x[i];
191: VecRestoreArrayRead(bb,&b);
192: VecRestoreArray(xx,&x);
193: VecRestoreArray(rr,&r);
194: return(0);
195: }
196: PetscErrorCode amult(Mat mat,Vec xx,Vec yy)
197: {
198: PetscInt i,n1;
199: PetscErrorCode ierr;
200: PetscScalar *y;
201: const PetscScalar *x;
204: VecGetSize(xx,&n1);
205: VecGetArrayRead(xx,&x);
206: VecGetArray(yy,&y);
207: n1--;
208: y[0] = -x[1] + 2.0*x[0];
209: y[n1] = -x[n1-1] + 2.0*x[n1];
210: for (i=1; i<n1; i++) y[i] = -x[i+1] - x[i-1] + 2.0*x[i];
211: VecRestoreArrayRead(xx,&x);
212: VecRestoreArray(yy,&y);
213: return(0);
214: }
215: /* --------------------------------------------------------------------- */
216: 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)
217: {
218: PetscInt i,n1;
219: PetscErrorCode ierr;
220: PetscScalar *x;
221: const PetscScalar *b;
224: *its = m;
225: *reason = PCRICHARDSON_CONVERGED_ITS;
226: VecGetSize(bb,&n1); n1--;
227: VecGetArrayRead(bb,&b);
228: VecGetArray(xx,&x);
229: while (m--) {
230: x[0] = .5*(x[1] + b[0]);
231: for (i=1; i<n1; i++) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
232: x[n1] = .5*(x[n1-1] + b[n1]);
233: for (i=n1-1; i>0; i--) x[i] = .5*(x[i+1] + x[i-1] + b[i]);
234: x[0] = .5*(x[1] + b[0]);
235: }
236: VecRestoreArrayRead(bb,&b);
237: VecRestoreArray(xx,&x);
238: return(0);
239: }
240: /* --------------------------------------------------------------------- */
241: 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)
242: {
243: PetscInt i,n,n1;
244: PetscErrorCode ierr;
245: PetscScalar *r,*x;
246: const PetscScalar *b;
249: *its = m;
250: *reason = PCRICHARDSON_CONVERGED_ITS;
251: VecGetSize(bb,&n); n1 = n - 1;
252: VecGetArrayRead(bb,&b);
253: VecGetArray(xx,&x);
254: VecGetArray(w,&r);
256: while (m--) {
257: r[0] = .5*(x[1] + b[0]);
258: for (i=1; i<n1; i++) r[i] = .5*(x[i+1] + x[i-1] + b[i]);
259: r[n1] = .5*(x[n1-1] + b[n1]);
260: for (i=0; i<n; i++) x[i] = (2.0*r[i] + x[i])/3.0;
261: }
262: VecRestoreArrayRead(bb,&b);
263: VecRestoreArray(xx,&x);
264: VecRestoreArray(w,&r);
265: return(0);
266: }
267: /*
268: We know for this application that yy and zz are the same
269: */
270: /* --------------------------------------------------------------------- */
271: PetscErrorCode interpolate(Mat mat,Vec xx,Vec yy,Vec zz)
272: {
273: PetscInt i,n,N,i2;
274: PetscErrorCode ierr;
275: PetscScalar *y;
276: const PetscScalar *x;
279: VecGetSize(yy,&N);
280: VecGetArrayRead(xx,&x);
281: VecGetArray(yy,&y);
282: n = N/2;
283: for (i=0; i<n; i++) {
284: i2 = 2*i;
285: y[i2] += .5*x[i];
286: y[i2+1] += x[i];
287: y[i2+2] += .5*x[i];
288: }
289: VecRestoreArrayRead(xx,&x);
290: VecRestoreArray(yy,&y);
291: return(0);
292: }
293: /* --------------------------------------------------------------------- */
294: PetscErrorCode restrct(Mat mat,Vec rr,Vec bb)
295: {
296: PetscInt i,n,N,i2;
297: PetscErrorCode ierr;
298: PetscScalar *b;
299: const PetscScalar *r;
302: VecGetSize(rr,&N);
303: VecGetArrayRead(rr,&r);
304: VecGetArray(bb,&b);
305: n = N/2;
307: for (i=0; i<n; i++) {
308: i2 = 2*i;
309: b[i] = (r[i2] + 2.0*r[i2+1] + r[i2+2]);
310: }
311: VecRestoreArrayRead(rr,&r);
312: VecRestoreArray(bb,&b);
313: return(0);
314: }
315: /* --------------------------------------------------------------------- */
316: PetscErrorCode Create1dLaplacian(PetscInt n,Mat *mat)
317: {
318: PetscScalar mone = -1.0,two = 2.0;
319: PetscInt i,idx;
323: MatCreateSeqAIJ(PETSC_COMM_SELF,n,n,3,NULL,mat);
325: idx = n-1;
326: MatSetValues(*mat,1,&idx,1,&idx,&two,INSERT_VALUES);
327: for (i=0; i<n-1; i++) {
328: MatSetValues(*mat,1,&i,1,&i,&two,INSERT_VALUES);
329: idx = i+1;
330: MatSetValues(*mat,1,&idx,1,&i,&mone,INSERT_VALUES);
331: MatSetValues(*mat,1,&i,1,&idx,&mone,INSERT_VALUES);
332: }
333: MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);
334: MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);
335: return(0);
336: }
337: /* --------------------------------------------------------------------- */
338: PetscErrorCode CalculateRhs(Vec u)
339: {
341: PetscInt i,n;
342: PetscReal h,x = 0.0;
343: PetscScalar uu;
346: VecGetSize(u,&n);
347: h = 1.0/((PetscReal)(n+1));
348: for (i=0; i<n; i++) {
349: x += h; uu = 2.0*h*h;
350: VecSetValues(u,1,&i,&uu,INSERT_VALUES);
351: }
352: return(0);
353: }
354: /* --------------------------------------------------------------------- */
355: PetscErrorCode CalculateSolution(PetscInt n,Vec *solution)
356: {
358: PetscInt i;
359: PetscReal h,x = 0.0;
360: PetscScalar uu;
363: VecCreateSeq(PETSC_COMM_SELF,n,solution);
364: h = 1.0/((PetscReal)(n+1));
365: for (i=0; i<n; i++) {
366: x += h; uu = x*(1.-x);
367: VecSetValues(*solution,1,&i,&uu,INSERT_VALUES);
368: }
369: return(0);
370: }
371: /* --------------------------------------------------------------------- */
372: PetscErrorCode CalculateError(Vec solution,Vec u,Vec r,PetscReal *e)
373: {
377: VecNorm(r,NORM_2,e+2);
378: VecWAXPY(r,-1.0,u,solution);
379: VecNorm(r,NORM_2,e);
380: VecNorm(r,NORM_1,e+1);
381: return(0);
382: }
387: /*TEST
389: test:
391: TEST*/