Actual source code: ex5.c

petsc-3.14.4 2021-02-03
Report Typos and Errors

  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*/