Actual source code: ex51.c

petsc-3.7.2 2016-06-05
Report Typos and Errors
  2: static char help[] = "This example solves a linear system in parallel with KSP.  The matrix\n\
  3: uses arbitrary order polynomials for finite elements on the unit square.  To test the parallel\n\
  4: matrix assembly, the matrix is intentionally laid out across processors\n\
  5: differently from the way it is assembled.  Input arguments are:\n\
  6:   -m <size> -p <order> : mesh size and polynomial order\n\n";

  8: /* Contributed by Travis Austin <austin@txcorp.com>, 2010.
  9:    based on src/ksp/ksp/examples/tutorials/ex3.c
 10:  */

 12: #include <petscksp.h>

 14: /* Declare user-defined routines */
 15: static PetscReal      src(PetscReal,PetscReal);
 16: static PetscReal      ubdy(PetscReal,PetscReal);
 17: static PetscReal      polyBasisFunc(PetscInt,PetscInt,PetscReal*,PetscReal);
 18: static PetscReal      derivPolyBasisFunc(PetscInt,PetscInt,PetscReal*,PetscReal);
 19: static PetscErrorCode Form1DElementMass(PetscReal,PetscInt,PetscReal*,PetscReal*,PetscScalar*);
 20: static PetscErrorCode Form1DElementStiffness(PetscReal,PetscInt,PetscReal*,PetscReal*,PetscScalar*);
 21: static PetscErrorCode Form2DElementMass(PetscInt,PetscScalar*,PetscScalar*);
 22: static PetscErrorCode Form2DElementStiffness(PetscInt,PetscScalar*,PetscScalar*,PetscScalar*);
 23: static PetscErrorCode FormNodalRhs(PetscInt,PetscReal,PetscReal,PetscReal,PetscReal*,PetscScalar*);
 24: static PetscErrorCode FormNodalSoln(PetscInt,PetscReal,PetscReal,PetscReal,PetscReal*,PetscScalar*);
 25: static void leggaulob(PetscReal, PetscReal, PetscReal [], PetscReal [], int);
 26: static void qAndLEvaluation(int, PetscReal, PetscReal*, PetscReal*, PetscReal*);

 30: int main(int argc,char **args)
 31: {
 32:   PetscInt       p = 2, m = 5;
 33:   PetscInt       num1Dnodes, num2Dnodes;
 34:   PetscScalar    *Ke1D,*Ke2D,*Me1D,*Me2D;
 35:   PetscScalar    *r,*ue,val;
 36:   Vec            u,ustar,b,q;
 37:   Mat            A,Mass;
 38:   KSP            ksp;
 39:   PetscInt       M,N;
 40:   PetscMPIInt    rank,size;
 41:   PetscReal      x,y,h,norm;
 42:   PetscInt       *idx,indx,count,*rows,i,j,k,start,end,its;
 43:   PetscReal      *rowsx,*rowsy;
 44:   PetscReal      *gllNode, *gllWgts;

 47:   PetscInitialize(&argc,&args,(char*)0,help);
 48:   PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Options for p-FEM","");
 49:   PetscOptionsInt("-m","Number of elements in each direction","None",m,&m,NULL);
 50:   PetscOptionsInt("-p","Order of each element (tensor product basis)","None",p,&p,NULL);
 51:   PetscOptionsEnd();
 52:   N    = (p*m+1)*(p*m+1); /* dimension of matrix */
 53:   M    = m*m; /* number of elements */
 54:   h    = 1.0/m; /* mesh width */
 55:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 56:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 58:   /* Create stiffness matrix */
 59:   MatCreate(PETSC_COMM_WORLD,&A);
 60:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 61:   MatSetFromOptions(A);
 62:   MatSetUp(A);
 63:   start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
 64:   end   = start + M/size + ((M%size) > rank);

 66:   /* Create matrix  */
 67:   MatCreate(PETSC_COMM_WORLD,&Mass);
 68:   MatSetSizes(Mass,PETSC_DECIDE,PETSC_DECIDE,N,N);
 69:   MatSetFromOptions(Mass);
 70:   MatSetUp(Mass);
 71:   start = rank*(M/size) + ((M%size) < rank ? (M%size) : rank);
 72:   end   = start + M/size + ((M%size) > rank);

 74:   /* Allocate element stiffness matrices */
 75:   num1Dnodes = (p+1);
 76:   num2Dnodes = num1Dnodes*num1Dnodes;

 78:   PetscMalloc1(num1Dnodes*num1Dnodes,&Me1D);
 79:   PetscMalloc1(num1Dnodes*num1Dnodes,&Ke1D);
 80:   PetscMalloc1(num2Dnodes*num2Dnodes,&Me2D);
 81:   PetscMalloc1(num2Dnodes*num2Dnodes,&Ke2D);
 82:   PetscMalloc1(num2Dnodes,&idx);
 83:   PetscMalloc1(num2Dnodes,&r);
 84:   PetscMalloc1(num2Dnodes,&ue);

 86:   /* Allocate quadrature and create stiffness matrices */
 87:   PetscMalloc1(p+1,&gllNode);
 88:   PetscMalloc1(p+1,&gllWgts);
 89:   leggaulob(0.0,1.0,gllNode,gllWgts,p); /* Get GLL nodes and weights */
 90:   Form1DElementMass(h,p,gllNode,gllWgts,Me1D);
 91:   Form1DElementStiffness(h,p,gllNode,gllWgts,Ke1D);
 92:   Form2DElementMass(p,Me1D,Me2D);
 93:   Form2DElementStiffness(p,Ke1D,Me1D,Ke2D);

 95:   /* Assemble matrix */
 96:   for (i=start; i<end; i++) {
 97:      indx = 0;
 98:      for (k=0;k<(p+1);++k) {
 99:        for (j=0;j<(p+1);++j) {
100:          idx[indx++] = p*(p*m+1)*(i/m) + p*(i % m) + k*(p*m+1) + j;
101:        }
102:      }
103:      MatSetValues(A,num2Dnodes,idx,num2Dnodes,idx,Ke2D,ADD_VALUES);
104:      MatSetValues(Mass,num2Dnodes,idx,num2Dnodes,idx,Me2D,ADD_VALUES);
105:   }
106:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
107:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
108:   MatAssemblyBegin(Mass,MAT_FINAL_ASSEMBLY);
109:   MatAssemblyEnd(Mass,MAT_FINAL_ASSEMBLY);

111:   PetscFree(Me1D);
112:   PetscFree(Ke1D);
113:   PetscFree(Me2D);
114:   PetscFree(Ke2D);

116:   /* Create right-hand-side and solution vectors */
117:   VecCreate(PETSC_COMM_WORLD,&u);
118:   VecSetSizes(u,PETSC_DECIDE,N);
119:   VecSetFromOptions(u);
120:   PetscObjectSetName((PetscObject)u,"Approx. Solution");
121:   VecDuplicate(u,&b);
122:   PetscObjectSetName((PetscObject)b,"Right hand side");
123:   VecDuplicate(u,&q);
124:   PetscObjectSetName((PetscObject)q,"Right hand side 2");
125:   VecDuplicate(b,&ustar);
126:   VecSet(u,0.0);
127:   VecSet(b,0.0);
128:   VecSet(q,0.0);

130:   /* Assemble nodal right-hand-side and soln vector  */
131:   for (i=start; i<end; i++) {
132:     x    = h*(i % m);
133:     y    = h*(i/m);
134:     indx = 0;
135:     for (k=0;k<(p+1);++k) {
136:       for (j=0;j<(p+1);++j) {
137:         idx[indx++] = p*(p*m+1)*(i/m) + p*(i % m) + k*(p*m+1) + j;
138:       }
139:     }
140:     FormNodalRhs(p,x,y,h,gllNode,r);
141:     FormNodalSoln(p,x,y,h,gllNode,ue);
142:     VecSetValues(q,num2Dnodes,idx,r,INSERT_VALUES);
143:     VecSetValues(ustar,num2Dnodes,idx,ue,INSERT_VALUES);
144:   }
145:   VecAssemblyBegin(q);
146:   VecAssemblyEnd(q);
147:   VecAssemblyBegin(ustar);
148:   VecAssemblyEnd(ustar);

150:   PetscFree(idx);
151:   PetscFree(r);
152:   PetscFree(ue);

154:   /* Get FE right-hand side vector */
155:   MatMult(Mass,q,b);

157:   /* Modify matrix and right-hand-side for Dirichlet boundary conditions */
158:   PetscMalloc1(4*p*m,&rows);
159:   PetscMalloc1(4*p*m,&rowsx);
160:   PetscMalloc1(4*p*m,&rowsy);
161:   for (i=0; i<p*m+1; i++) {
162:     rows[i]          = i; /* bottom */
163:     rowsx[i]         = (i/p)*h+gllNode[i%p]*h;
164:     rowsy[i]         = 0.0;
165:     rows[3*p*m-1+i]  = (p*m)*(p*m+1) + i; /* top */
166:     rowsx[3*p*m-1+i] = (i/p)*h+gllNode[i%p]*h;
167:     rowsy[3*p*m-1+i] = 1.0;
168:   }
169:   count = p*m+1; /* left side */
170:   indx  = 1;
171:   for (i=p*m+1; i<(p*m)*(p*m+1); i+= (p*m+1)) {
172:     rows[count]    = i;
173:     rowsx[count]   = 0.0;
174:     rowsy[count++] = (indx/p)*h+gllNode[indx%p]*h;
175:     indx++;
176:   }
177:   count = 2*p*m; /* right side */
178:   indx  = 1;
179:   for (i=2*p*m+1; i<(p*m)*(p*m+1); i+= (p*m+1)) {
180:     rows[count]    = i;
181:     rowsx[count]   = 1.0;
182:     rowsy[count++] = (indx/p)*h+gllNode[indx%p]*h;
183:     indx++;
184:   }
185:   for (i=0; i<4*p*m; i++) {
186:     x    = rowsx[i];
187:     y    = rowsy[i];
188:     val  = ubdy(x,y);
189:     VecSetValues(b,1,&rows[i],&val,INSERT_VALUES);
190:     VecSetValues(u,1,&rows[i],&val,INSERT_VALUES);
191:   }
192:   MatZeroRows(A,4*p*m,rows,1.0,0,0);
193:   PetscFree(rows);
194:   PetscFree(rowsx);
195:   PetscFree(rowsy);

197:   VecAssemblyBegin(u);
198:   VecAssemblyEnd(u);
199:   VecAssemblyBegin(b);
200:   VecAssemblyEnd(b);

202:   /* Solve linear system */
203:   KSPCreate(PETSC_COMM_WORLD,&ksp);
204:   KSPSetOperators(ksp,A,A);
205:   KSPSetInitialGuessNonzero(ksp,PETSC_TRUE);
206:   KSPSetFromOptions(ksp);
207:   KSPSolve(ksp,b,u);

209:   /* Check error */
210:   VecAXPY(u,-1.0,ustar);
211:   VecNorm(u,NORM_2,&norm);
212:   KSPGetIterationNumber(ksp,&its);
213:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g Iterations %D\n",(double)(norm*h),its);

215:   PetscFree(gllNode);
216:   PetscFree(gllWgts);

218:   KSPDestroy(&ksp);
219:   VecDestroy(&u);
220:   VecDestroy(&b);
221:   VecDestroy(&q);
222:   VecDestroy(&ustar);
223:   MatDestroy(&A);
224:   MatDestroy(&Mass);

226:   PetscFinalize();
227:   return 0;
228: }

230: /* --------------------------------------------------------------------- */

234: /* 1d element stiffness mass matrix  */
235: static PetscErrorCode Form1DElementMass(PetscReal H,PetscInt P,PetscReal *gqn,PetscReal *gqw,PetscScalar *Me1D)
236: {
237:   PetscInt i,j,k;
238:   PetscInt indx;

241:   for (j=0; j<(P+1); ++j) {
242:     for (i=0; i<(P+1); ++i) {
243:       indx       = j*(P+1)+i;
244:       Me1D[indx] = 0.0;
245:       for (k=0; k<(P+1);++k) {
246:         Me1D[indx] += H*gqw[k]*polyBasisFunc(P,i,gqn,gqn[k])*polyBasisFunc(P,j,gqn,gqn[k]);
247:       }
248:     }
249:   }
250:   return(0);
251: }

253: /* --------------------------------------------------------------------- */

257: /* 1d element stiffness matrix for derivative */
258: static PetscErrorCode Form1DElementStiffness(PetscReal H,PetscInt P,PetscReal *gqn,PetscReal *gqw,PetscScalar *Ke1D)
259: {
260:   PetscInt i,j,k;
261:   PetscInt indx;

264:   for (j=0;j<(P+1);++j) {
265:     for (i=0;i<(P+1);++i) {
266:       indx = j*(P+1)+i;
267:       Ke1D[indx] = 0.0;
268:       for (k=0; k<(P+1);++k) {
269:         Ke1D[indx] += (1./H)*gqw[k]*derivPolyBasisFunc(P,i,gqn,gqn[k])*derivPolyBasisFunc(P,j,gqn,gqn[k]);
270:       }
271:     }
272:   }
273:   return(0);
274: }

276: /* --------------------------------------------------------------------- */

280:    /* element mass matrix */
281: static PetscErrorCode Form2DElementMass(PetscInt P,PetscScalar *Me1D,PetscScalar *Me2D)
282: {
283:   PetscInt i1,j1,i2,j2;
284:   PetscInt indx1,indx2,indx3;

287:   for (j2=0;j2<(P+1);++j2) {
288:     for (i2=0; i2<(P+1);++i2) {
289:       for (j1=0;j1<(P+1);++j1) {
290:         for (i1=0;i1<(P+1);++i1) {
291:           indx1 = j1*(P+1)+i1;
292:           indx2 = j2*(P+1)+i2;
293:           indx3 = (j2*(P+1)+j1)*(P+1)*(P+1)+(i2*(P+1)+i1);
294:           Me2D[indx3] = Me1D[indx1]*Me1D[indx2];
295:         }
296:       }
297:     }
298:   }
299:   return(0);
300: }

302: /* --------------------------------------------------------------------- */

306: /* element stiffness for Laplacian */
307: static PetscErrorCode Form2DElementStiffness(PetscInt P,PetscScalar *Ke1D,PetscScalar *Me1D,PetscScalar *Ke2D)
308: {
309:   PetscInt i1,j1,i2,j2;
310:   PetscInt indx1,indx2,indx3;

313:   for (j2=0;j2<(P+1);++j2) {
314:     for (i2=0; i2<(P+1);++i2) {
315:       for (j1=0;j1<(P+1);++j1) {
316:         for (i1=0;i1<(P+1);++i1) {
317:           indx1 = j1*(P+1)+i1;
318:           indx2 = j2*(P+1)+i2;
319:           indx3 = (j2*(P+1)+j1)*(P+1)*(P+1)+(i2*(P+1)+i1);
320:           Ke2D[indx3] = Ke1D[indx1]*Me1D[indx2] + Me1D[indx1]*Ke1D[indx2];
321:         }
322:       }
323:     }
324:   }
325:   return(0);
326: }

328: /* --------------------------------------------------------------------- */

332: static PetscErrorCode FormNodalRhs(PetscInt P,PetscReal x,PetscReal y,PetscReal H,PetscReal* nds,PetscScalar *r)
333: {
334:   PetscInt i,j,indx;

337:   indx=0;
338:   for (j=0;j<(P+1);++j) {
339:     for (i=0;i<(P+1);++i) {
340:       r[indx] = src(x+H*nds[i],y+H*nds[j]);
341:       indx++;
342:     }
343:   }
344:   return(0);
345: }

347: /* --------------------------------------------------------------------- */

351: static PetscErrorCode FormNodalSoln(PetscInt P,PetscReal x,PetscReal y,PetscReal H,PetscReal* nds,PetscScalar *u)
352: {
353:   PetscInt i,j,indx;

356:   indx=0;
357:   for (j=0;j<(P+1);++j) {
358:     for (i=0;i<(P+1);++i) {
359:       u[indx] = ubdy(x+H*nds[i],y+H*nds[j]);
360:       indx++;
361:     }
362:   }
363:   return(0);
364: }

366: /* --------------------------------------------------------------------- */

370: static PetscReal polyBasisFunc(PetscInt order, PetscInt basis, PetscReal *xLocVal, PetscReal xval)
371: {
372:   PetscReal denominator = 1.;
373:   PetscReal numerator   = 1.;
374:   PetscInt  i           =0;

376:   for (i=0; i<(order+1); i++) {
377:     if (i!=basis) {
378:       numerator   *= (xval-xLocVal[i]);
379:       denominator *= (xLocVal[basis]-xLocVal[i]);
380:     }
381:   }
382:   return (numerator/denominator);
383: }

385: /* --------------------------------------------------------------------- */

389: static PetscReal derivPolyBasisFunc(PetscInt order, PetscInt basis, PetscReal *xLocVal, PetscReal xval)
390: {
391:   PetscReal denominator;
392:   PetscReal numerator;
393:   PetscReal numtmp;
394:   PetscInt  i=0,j=0;

396:   denominator=1.;
397:   for (i=0; i<(order+1); i++) {
398:     if (i!=basis) denominator *= (xLocVal[basis]-xLocVal[i]);
399:   }
400:   numerator = 0.;
401:   for (j=0;j<(order+1);++j) {
402:     if (j != basis) {
403:       numtmp = 1.;
404:       for (i=0;i<(order+1);++i) {
405:         if (i!=basis && j!=i) numtmp *= (xval-xLocVal[i]);
406:       }
407:       numerator += numtmp;
408:     }
409:   }

411:   return (numerator/denominator);
412: }

414: /* --------------------------------------------------------------------- */

416: static PetscReal ubdy(PetscReal x,PetscReal y)
417: {
418:   return x*x*y*y;
419: }

421: static PetscReal src(PetscReal x,PetscReal y)
422: {
423:   return -2.*y*y - 2.*x*x;
424: }
425: /* --------------------------------------------------------------------- */

427: static void leggaulob(PetscReal x1, PetscReal x2, PetscReal x[], PetscReal w[], int n)
428: /*******************************************************************************
429: Given the lower and upper limits of integration x1 and x2, and given n, this
430: routine returns arrays x[0..n-1] and w[0..n-1] of length n, containing the abscissas
431: and weights of the Gauss-Lobatto-Legendre n-point quadrature formula.
432: *******************************************************************************/
433: {
434:   PetscInt    j,m;
435:   PetscReal z1,z,xm,xl,q,qp,Ln,scale;
436:   if (n==1) {
437:     x[0] = x1;   /* Scale the root to the desired interval, */
438:     x[1] = x2;   /* and put in its symmetric counterpart.   */
439:     w[0] = 1.;   /* Compute the weight */
440:     w[1] = 1.;   /* and its symmetric counterpart. */
441:   } else {
442:     x[0] = x1;   /* Scale the root to the desired interval, */
443:     x[n] = x2;   /* and put in its symmetric counterpart.   */
444:     w[0] = 2./(n*(n+1));;   /* Compute the weight */
445:     w[n] = 2./(n*(n+1));   /* and its symmetric counterpart. */
446:     m    = (n+1)/2; /* The roots are symmetric, so we only find half of them. */
447:     xm   = 0.5*(x2+x1);
448:     xl   = 0.5*(x2-x1);
449:     for (j=1; j<=(m-1); j++) { /* Loop over the desired roots. */
450:       z=-1.0*PetscCosReal((PETSC_PI*(j+0.25)/(n))-(3.0/(8.0*n*PETSC_PI))*(1.0/(j+0.25)));
451:       /* Starting with the above approximation to the ith root, we enter */
452:       /* the main loop of refinement by Newton's method.                 */
453:       do {
454:         qAndLEvaluation(n,z,&q,&qp,&Ln);
455:         z1 = z;
456:         z  = z1-q/qp; /* Newton's method. */
457:       } while (fabs(z-z1) > 3.0e-11);
458:       qAndLEvaluation(n,z,&q,&qp,&Ln);
459:       x[j]   = xm+xl*z;      /* Scale the root to the desired interval, */
460:       x[n-j] = xm-xl*z;      /* and put in its symmetric counterpart.   */
461:       w[j]   = 2.0/(n*(n+1)*Ln*Ln);  /* Compute the weight */
462:       w[n-j] = w[j];                 /* and its symmetric counterpart. */
463:     }
464:   }
465:   if (n%2==0) {
466:     qAndLEvaluation(n,0.0,&q,&qp,&Ln);
467:     x[n/2]=(x2-x1)/2.0;
468:     w[n/2]=2.0/(n*(n+1)*Ln*Ln);
469:   }
470:   /* scale the weights according to mapping from [-1,1] to [0,1] */
471:   scale = (x2-x1)/2.0;
472:   for (j=0; j<=n; ++j) w[j] = w[j]*scale;
473: }


476: /******************************************************************************/
477: static void qAndLEvaluation(int n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln)
478: /*******************************************************************************
479: Compute the polynomial qn(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in
480: addition to L_N(x) as these are needed for the GLL points.  See the book titled
481: "Implementing Spectral Methods for Partial Differential Equations: Algorithms
482: for Scientists and Engineers" by David A. Kopriva.
483: *******************************************************************************/
484: {
485:   PetscInt k;

487:   PetscReal Lnp;
488:   PetscReal Lnp1, Lnp1p;
489:   PetscReal Lnm1, Lnm1p;
490:   PetscReal Lnm2, Lnm2p;

492:   Lnm1  = 1.0;
493:   *Ln   = x;
494:   Lnm1p = 0.0;
495:   Lnp   = 1.0;

497:   for (k=2; k<=n; ++k) {
498:     Lnm2  = Lnm1;
499:     Lnm1  = *Ln;
500:     Lnm2p = Lnm1p;
501:     Lnm1p = Lnp;
502:     *Ln   = (2.*k-1.)/(1.0*k)*x*Lnm1 - (k-1.)/(1.0*k)*Lnm2;
503:     Lnp   = Lnm2p + (2.0*k-1)*Lnm1;
504:   }
505:   k     = n+1;
506:   Lnp1  = (2.*k-1.)/(1.0*k)*x*(*Ln) - (k-1.)/(1.0*k)*Lnm1;
507:   Lnp1p = Lnm1p + (2.0*k-1)*(*Ln);
508:   *q    = Lnp1 - Lnm1;
509:   *qp   = Lnp1p - Lnm1p;
510: }