Actual source code: petscgll.c
petsc-3.10.5 2019-03-28
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: }