Actual source code: petscgll.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

  2:  #include <petscgll.h>
  3:  #include <petscviewer.h>
  4:  #include <petscblaslapack.h>
  5:  #include <petsc/private/petscimpl.h>


  8: static void qAndLEvaluation(PetscInt n, PetscReal x, PetscReal *q, PetscReal *qp, PetscReal *Ln)
  9: /*
 10:   Compute the polynomial q(x) = L_{N+1}(x) - L_{n-1}(x) and its derivative in
 11:   addition to L_N(x) as these are needed for computing the GLL points via Newton's method.
 12:   Reference: "Implementing Spectral Methods for Partial Differential Equations: Algorithms
 13:   for Scientists and Engineers" by David A. Kopriva.
 14: */
 15: {
 16:   PetscInt k;

 18:   PetscReal Lnp;
 19:   PetscReal Lnp1, Lnp1p;
 20:   PetscReal Lnm1, Lnm1p;
 21:   PetscReal Lnm2, Lnm2p;

 23:   Lnm1  = 1.0;
 24:   *Ln   = x;
 25:   Lnm1p = 0.0;
 26:   Lnp   = 1.0;

 28:   for (k=2; k<=n; ++k) {
 29:     Lnm2  = Lnm1;
 30:     Lnm1  = *Ln;
 31:     Lnm2p = Lnm1p;
 32:     Lnm1p = Lnp;
 33:     *Ln   = (2.*((PetscReal)k)-1.)/(1.0*((PetscReal)k))*x*Lnm1 - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm2;
 34:     Lnp   = Lnm2p + (2.0*((PetscReal)k)-1.)*Lnm1;
 35:   }
 36:   k     = n+1;
 37:   Lnp1  = (2.*((PetscReal)k)-1.)/(((PetscReal)k))*x*(*Ln) - (((PetscReal)k)-1.)/((PetscReal)k)*Lnm1;
 38:   Lnp1p = Lnm1p + (2.0*((PetscReal)k)-1.)*(*Ln);
 39:   *q    = Lnp1 - Lnm1;
 40:   *qp   = Lnp1p - Lnm1p;
 41: }

 43: /*@C
 44:    PetscGLLCreate - creates a set of the locations and weights of the Gauss-Lobatto-Legendre (GLL) nodes of a given size
 45:                       on the domain [-1,1]

 47:    Not Collective

 49:    Input Parameter:
 50: +  n - number of grid nodes
 51: -  type - PETSCGLL_VIA_LINEARALGEBRA or PETSCGLL_VIA_NEWTON

 53:    Output Parameter:
 54: .  gll - the nodes

 56:    Notes:
 57:     For n > 30  the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
 58:           close enough to the desired solution

 60:    These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes

 62:    See  http://epubs.siam.org/doi/abs/10.1137/110855442  http://epubs.siam.org/doi/abs/10.1137/120889873 for better ways to compute GLL nodes

 64:    Level: beginner

 66: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLIntegrate(), PetscGLLElementLaplacianCreate(), PetscGLLElementLaplacianDestroy(), 
 67:           PetscGLLElementGradientCreate(), PetscGLLElementGradientDestroy(), PetscGLLElementAdvectionCreate(), PetscGLLElementAdvectionDestroy()

 69: @*/
 70: PetscErrorCode PetscGLLCreate(PetscInt n,PetscGLLCreateType type,PetscGLL *gll)
 71: {

 75:   if (n < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
 76:   PetscMalloc2(n,&gll->nodes,n,&gll->weights);

 78:   if (type == PETSCGLL_VIA_LINEARALGEBRA) {
 79:     PetscReal      *M,si;
 80:     PetscBLASInt   bn,lierr;
 81:     PetscReal      x,z0,z1,z2;
 82:     PetscInt       i,p = n - 1,nn;

 84:     gll->nodes[0]   =-1.0;
 85:     gll->nodes[n-1] = 1.0;
 86:     if (n-2 > 0){
 87:       PetscMalloc1(n-1,&M);
 88:       for (i=0; i<n-2; i++) {
 89:         si  = ((PetscReal)i)+1.0;
 90:         M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
 91:       }
 92:       PetscBLASIntCast(n-2,&bn);
 93:       PetscMemzero(&gll->nodes[1],bn*sizeof(gll->nodes[1]));
 94:       PetscFPTrapPush(PETSC_FP_TRAP_OFF);
 95:       x=0;
 96:       PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&gll->nodes[1],M,&x,&bn,M,&lierr));
 97:       if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
 98:       PetscFPTrapPop();
 99:       PetscFree(M);
100:     }
101:     if ((n-1)%2==0) {
102:       gll->nodes[(n-1)/2]   = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
103:     }

105:     gll->weights[0] = gll->weights[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
106:     z2 = -1.;                      /* Dummy value to avoid -Wmaybe-initialized */
107:     for (i=1; i<p; i++) {
108:       x  = gll->nodes[i];
109:       z0 = 1.0;
110:       z1 = x;
111:       for (nn=1; nn<p; nn++) {
112:         z2 = x*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
113:         z0 = z1;
114:         z1 = z2;
115:       }
116:       gll->weights[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
117:     }
118:   } else {
119:     PetscInt  j,m;
120:     PetscReal z1,z,q,qp,Ln;
121:     PetscReal *pt;
122:     PetscMalloc1(n,&pt);

124:     if (n > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGLL_VIA_NEWTON produces incorrect answers for n > 30");
125:     gll->nodes[0]     = -1.0;
126:     gll->nodes[n-1]   = 1.0;
127:     gll->weights[0]   = gll->weights[n-1] = 2./(((PetscReal)n)*(((PetscReal)n)-1.0));;
128:     m  = (n-1)/2; /* The roots are symmetric, so we only find half of them. */
129:     for (j=1; j<=m; j++) { /* Loop over the desired roots. */
130:       z = -1.0*PetscCosReal((PETSC_PI*((PetscReal)j)+0.25)/(((PetscReal)n)-1.0))-(3.0/(8.0*(((PetscReal)n)-1.0)*PETSC_PI))*(1.0/(((PetscReal)j)+0.25));
131:       /* Starting with the above approximation to the ith root, we enter */
132:       /* the main loop of refinement by Newton's method.                 */
133:       do {
134:         qAndLEvaluation(n-1,z,&q,&qp,&Ln);
135:         z1 = z;
136:         z  = z1-q/qp; /* Newton's method. */
137:       } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
138:       qAndLEvaluation(n-1,z,&q,&qp,&Ln);

140:       gll->nodes[j]       = z;
141:       gll->nodes[n-1-j]   = -z;      /* and put in its symmetric counterpart.   */
142:       gll->weights[j]     = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln);  /* Compute the weight */
143:       gll->weights[n-1-j] = gll->weights[j];                 /* and its symmetric counterpart. */
144:       pt[j]=qp;
145:     }

147:     if ((n-1)%2==0) {
148:       qAndLEvaluation(n-1,0.0,&q,&qp,&Ln);
149:       gll->nodes[(n-1)/2]   = 0.0;
150:       gll->weights[(n-1)/2] = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln);
151:     }
152:     PetscFree(pt);
153:   }
154:   gll->n = n;
155:   return(0);
156: }

158: /*@C
159:    PetscGLLDestroy - destroys a set of GLL nodes and weights

161:    Not Collective

163:    Input Parameter:
164: .  gll - the nodes

166:    Level: beginner

168: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLView()

170: @*/
171: PetscErrorCode PetscGLLDestroy(PetscGLL *gll)
172: {

176:   PetscFree2(gll->nodes,gll->weights);
177:   gll->n = 0;
178:   return(0);
179: }

181: /*@C
182:    PetscGLLView - views a set of GLL nodes

184:    Not Collective

186:    Input Parameter:
187: +  gll - the nodes
188: .  viewer - the viewer

190:    Level: beginner

192: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()

194: @*/
195: PetscErrorCode PetscGLLView(PetscGLL *gll,PetscViewer viewer)
196: {
197:   PetscErrorCode    ierr;
198:   PetscBool         iascii;

200:   PetscInt          i;

203:   if (!viewer) {
204:     PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);
205:   }
207:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
208:   if (iascii) {
209:     PetscViewerASCIIPrintf(viewer,"%D Gauss-Lobatto-Legendre (GLL) nodes and weights\n",gll->n);
210:     for (i=0; i<gll->n; i++) {
211:       PetscViewerASCIIPrintf(viewer,"  %D %16.14e %16.14e\n",i,(double)gll->nodes[i],(double)gll->weights[i]);
212:     }
213:   }
214:   return(0);
215: }

217: /*@C
218:    PetscGLLIntegrate - Compute the L2 integral of a function on the GLL points

220:    Not Collective

222:    Input Parameter:
223: +  gll - the nodes
224: .  f - the function values at the nodes

226:    Output Parameter:
227: .  in - the value of the integral

229:    Level: beginner

231: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()

233: @*/
234: PetscErrorCode PetscGLLIntegrate(PetscGLL *gll,const PetscReal *f,PetscReal *in)
235: {
236:   PetscInt          i;

239:   *in = 0.;
240:   for (i=0; i<gll->n; i++) {
241:     *in += f[i]*f[i]*gll->weights[i];
242:   }
243:   return(0);
244: }

246: /*@C
247:    PetscGLLElementLaplacianCreate - computes the Laplacian for a single 1d GLL element

249:    Not Collective

251:    Input Parameter:
252: .  gll - the nodes

254:    Output Parameter:
255: .  A - the stiffness element

257:    Level: beginner

259:    Notes:
260:     Destroy this with PetscGLLElementLaplacianDestroy()

262:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented (the array is symmetric)

264: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()

266: @*/
267: PetscErrorCode PetscGLLElementLaplacianCreate(PetscGLL *gll,PetscReal ***AA)
268: {
269:   PetscReal        **A;
270:   PetscErrorCode  ierr;
271:   const PetscReal  *nodes = gll->nodes;
272:   const PetscInt   n = gll->n, p = gll->n-1;
273:   PetscReal        z0,z1,z2 = -1,x,Lpj,Lpr;
274:   PetscInt         i,j,nn,r;

277:   PetscMalloc1(n,&A);
278:   PetscMalloc1(n*n,&A[0]);
279:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

281:   for (j=1; j<p; j++) {
282:     x  = nodes[j];
283:     z0 = 1.;
284:     z1 = x;
285:     for (nn=1; nn<p; nn++) {
286:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
287:       z0 = z1;
288:       z1 = z2;
289:     }
290:     Lpj=z2;
291:     for (r=1; r<p; r++) {
292:       if (r == j) {
293:         A[j][j]=2./(3.*(1.-nodes[j]*nodes[j])*Lpj*Lpj);
294:       } else {
295:         x  = nodes[r];
296:         z0 = 1.;
297:         z1 = x;
298:         for (nn=1; nn<p; nn++) {
299:           z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
300:           z0 = z1;
301:           z1 = z2;
302:         }
303:         Lpr     = z2;
304:         A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(nodes[j]-nodes[r])*(nodes[j]-nodes[r]));
305:       }
306:     }
307:   }
308:   for (j=1; j<p+1; j++) {
309:     x  = nodes[j];
310:     z0 = 1.;
311:     z1 = x;
312:     for (nn=1; nn<p; nn++) {
313:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
314:       z0 = z1;
315:       z1 = z2;
316:     }
317:     Lpj     = z2;
318:     A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+nodes[j])*(1.+nodes[j]));
319:     A[0][j] = A[j][0];
320:   }
321:   for (j=0; j<p; j++) {
322:     x  = nodes[j];
323:     z0 = 1.;
324:     z1 = x;
325:     for (nn=1; nn<p; nn++) {
326:       z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
327:       z0 = z1;
328:       z1 = z2;
329:     }
330:     Lpj=z2;

332:     A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-nodes[j])*(1.-nodes[j]));
333:     A[j][p] = A[p][j];
334:   }
335:   A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
336:   A[p][p]=A[0][0];
337:   *AA = A;
338:   return(0);
339: }

341: /*@C
342:    PetscGLLElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element

344:    Not Collective

346:    Input Parameter:
347: +  gll - the nodes
348: -  A - the stiffness element

350:    Level: beginner

352: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate()

354: @*/
355: PetscErrorCode PetscGLLElementLaplacianDestroy(PetscGLL *gll,PetscReal ***AA)
356: {

360:   PetscFree((*AA)[0]);
361:   PetscFree(*AA);
362:   *AA  = NULL;
363:   return(0);
364: }

366: /*@C
367:    PetscGLLElementGradientCreate - computes the gradient for a single 1d GLL element

369:    Not Collective

371:    Input Parameter:
372: .  gll - the nodes

374:    Output Parameter:
375: .  AA - the stiffness element
376: -  AAT - the transpose of AA (pass in NULL if you do not need this array)

378:    Level: beginner

380:    Notes:
381:     Destroy this with PetscGLLElementGradientDestroy()

383:    You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented

385: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()

387: @*/
388: PetscErrorCode PetscGLLElementGradientCreate(PetscGLL *gll,PetscReal ***AA, PetscReal ***AAT)
389: {
390:   PetscReal        **A, **AT = NULL;
391:   PetscErrorCode  ierr;
392:   const PetscReal  *nodes = gll->nodes;
393:   const PetscInt   n = gll->n, p = gll->n-1;
394:   PetscReal        q,qp,Li, Lj,d0;
395:   PetscInt         i,j;

398:   PetscMalloc1(n,&A);
399:   PetscMalloc1(n*n,&A[0]);
400:   for (i=1; i<n; i++) A[i] = A[i-1]+n;

402:   if (AAT) {
403:     PetscMalloc1(n,&AT);
404:     PetscMalloc1(n*n,&AT[0]);
405:     for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
406:   }

408:   if (n==1) {A[0][0] = 0.;}
409:   d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
410:   for  (i=0; i<n; i++) {
411:     for  (j=0; j<n; j++) {
412:       A[i][j] = 0.;
413:       qAndLEvaluation(p,nodes[i],&q,&qp,&Li);
414:       qAndLEvaluation(p,nodes[j],&q,&qp,&Lj);
415:       if (i!=j)             A[i][j] = Li/(Lj*(nodes[i]-nodes[j]));
416:       if ((j==i) && (i==0)) A[i][j] = -d0;
417:       if (j==i && i==p)     A[i][j] = d0;
418:       if (AT) AT[j][i] = A[i][j];
419:     }
420:   }
421:   if (AAT) *AAT = AT;
422:   *AA  = A;
423:   return(0);
424: }

426: /*@C
427:    PetscGLLElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGLLElementGradientCreate()

429:    Not Collective

431:    Input Parameter:
432: +  gll - the nodes
433: .  AA - the stiffness element
434: -  AAT - the transpose of the element

436:    Level: beginner

438: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()

440: @*/
441: PetscErrorCode PetscGLLElementGradientDestroy(PetscGLL *gll,PetscReal ***AA,PetscReal ***AAT)
442: {

446:   PetscFree((*AA)[0]);
447:   PetscFree(*AA);
448:   *AA  = NULL;
449:   if (*AAT) {
450:     PetscFree((*AAT)[0]);
451:     PetscFree(*AAT);
452:     *AAT  = NULL;
453:   }
454:   return(0);
455: }

457: /*@C
458:    PetscGLLElementAdvectionCreate - computes the advection operator for a single 1d GLL element

460:    Not Collective

462:    Input Parameter:
463: .  gll - the nodes

465:    Output Parameter:
466: .  AA - the stiffness element

468:    Level: beginner

470:    Notes:
471:     Destroy this with PetscGLLElementAdvectionDestroy()

473:    This is the same as the Gradient operator multiplied by the diagonal mass matrix

475:    You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented

477: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()

479: @*/
480: PetscErrorCode PetscGLLElementAdvectionCreate(PetscGLL *gll,PetscReal ***AA)
481: {
482:   PetscReal       **D;
483:   PetscErrorCode  ierr;
484:   const PetscReal  *weights = gll->weights;
485:   const PetscInt   n = gll->n;
486:   PetscInt         i,j;

489:   PetscGLLElementGradientCreate(gll,&D,NULL);
490:   for (i=0; i<n; i++){
491:     for (j=0; j<n; j++) {
492:       D[i][j] = weights[i]*D[i][j];
493:     }
494:   }
495:   *AA = D;
496:   return(0);
497: }

499: /*@C
500:    PetscGLLElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element

502:    Not Collective

504:    Input Parameter:
505: +  gll - the nodes
506: -  A - advection

508:    Level: beginner

510: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()

512: @*/
513: PetscErrorCode PetscGLLElementAdvectionDestroy(PetscGLL *gll,PetscReal ***AA)
514: {

518:   PetscFree((*AA)[0]);
519:   PetscFree(*AA);
520:   *AA  = NULL;
521:   return(0);
522: }

524: PetscErrorCode PetscGLLElementMassCreate(PetscGLL *gll,PetscReal ***AA)
525: {
526:   PetscReal        **A;
527:   PetscErrorCode  ierr;
528:   const PetscReal  *weights = gll->weights;
529:   const PetscInt   n = gll->n;
530:   PetscInt         i,j;

533:   PetscMalloc1(n,&A);
534:   PetscMalloc1(n*n,&A[0]);
535:   for (i=1; i<n; i++) A[i] = A[i-1]+n;
536:   if (n==1) {A[0][0] = 0.;}
537:   for  (i=0; i<n; i++) {
538:     for  (j=0; j<n; j++) {
539:       A[i][j] = 0.;
540:       if (j==i)     A[i][j] = weights[i];
541:     }
542:   }
543:   *AA  = A;
544:   return(0);
545: }

547:   PetscErrorCode PetscGLLElementMassDestroy(PetscGLL *gll,PetscReal ***AA)
548: {

552:   PetscFree((*AA)[0]);
553:   PetscFree(*AA);
554:   *AA  = NULL;
555:   return(0);
556: }