Actual source code: petscgll.c
petsc-3.9.4 2018-09-11
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: For n > 30 the Newton approach computes duplicate (incorrect) values for some nodes because the initial guess is apparently not
57: close enough to the desired solution
59: These are useful for implementing spectral methods based on Gauss-Lobatto-Legendre (GLL) nodes
61: 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
63: Level: beginner
65: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLIntegrate(), PetscGLLElementLaplacianCreate(), PetscGLLElementLaplacianDestroy(),
66: PetscGLLElementGradientCreate(), PetscGLLElementGradientDestroy(), PetscGLLElementAdvectionCreate(), PetscGLLElementAdvectionDestroy()
68: @*/
69: PetscErrorCode PetscGLLCreate(PetscInt n,PetscGLLCreateType type,PetscGLL *gll)
70: {
74: if (n < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide at least 2 grid points per element");
75: PetscMalloc2(n,&gll->nodes,n,&gll->weights);
77: if (type == PETSCGLL_VIA_LINEARALGEBRA) {
78: PetscReal *M,si;
79: PetscBLASInt bn,lierr;
80: PetscReal x,z0,z1,z2;
81: PetscInt i,p = n - 1,nn;
83: gll->nodes[0] =-1.0;
84: gll->nodes[n-1] = 1.0;
85: if (n-2 > 0){
86: PetscMalloc1(n-1,&M);
87: for (i=0; i<n-2; i++) {
88: si = ((PetscReal)i)+1.0;
89: M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
90: }
91: PetscBLASIntCast(n-2,&bn);
92: PetscMemzero(&gll->nodes[1],bn*sizeof(gll->nodes[1]));
93: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
94: x=0;
95: PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&gll->nodes[1],M,&x,&bn,M,&lierr));
96: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
97: PetscFPTrapPop();
98: PetscFree(M);
99: }
100: if ((n-1)%2==0) {
101: gll->nodes[(n-1)/2] = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
102: }
104: gll->weights[0] = gll->weights[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
105: z2 = -1.; /* Dummy value to avoid -Wmaybe-initialized */
106: for (i=1; i<p; i++) {
107: x = gll->nodes[i];
108: z0 = 1.0;
109: z1 = x;
110: for (nn=1; nn<p; nn++) {
111: z2 = x*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
112: z0 = z1;
113: z1 = z2;
114: }
115: gll->weights[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
116: }
117: } else {
118: PetscInt j,m;
119: PetscReal z1,z,q,qp,Ln;
120: PetscReal *pt;
121: PetscMalloc1(n,&pt);
123: if (n > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGLL_VIA_NEWTON produces incorrect answers for n > 30");
124: gll->nodes[0] = -1.0;
125: gll->nodes[n-1] = 1.0;
126: gll->weights[0] = gll->weights[n-1] = 2./(((PetscReal)n)*(((PetscReal)n)-1.0));;
127: m = (n-1)/2; /* The roots are symmetric, so we only find half of them. */
128: for (j=1; j<=m; j++) { /* Loop over the desired roots. */
129: 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));
130: /* Starting with the above approximation to the ith root, we enter */
131: /* the main loop of refinement by Newton's method. */
132: do {
133: qAndLEvaluation(n-1,z,&q,&qp,&Ln);
134: z1 = z;
135: z = z1-q/qp; /* Newton's method. */
136: } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
137: qAndLEvaluation(n-1,z,&q,&qp,&Ln);
139: gll->nodes[j] = z;
140: gll->nodes[n-1-j] = -z; /* and put in its symmetric counterpart. */
141: gll->weights[j] = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln); /* Compute the weight */
142: gll->weights[n-1-j] = gll->weights[j]; /* and its symmetric counterpart. */
143: pt[j]=qp;
144: }
146: if ((n-1)%2==0) {
147: qAndLEvaluation(n-1,0.0,&q,&qp,&Ln);
148: gll->nodes[(n-1)/2] = 0.0;
149: gll->weights[(n-1)/2] = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln);
150: }
151: PetscFree(pt);
152: }
153: gll->n = n;
154: return(0);
155: }
157: /*@C
158: PetscGLLDestroy - destroys a set of GLL nodes and weights
160: Not Collective
162: Input Parameter:
163: . gll - the nodes
165: Level: beginner
167: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLView()
169: @*/
170: PetscErrorCode PetscGLLDestroy(PetscGLL *gll)
171: {
175: PetscFree2(gll->nodes,gll->weights);
176: gll->n = 0;
177: return(0);
178: }
180: /*@C
181: PetscGLLView - views a set of GLL nodes
183: Not Collective
185: Input Parameter:
186: + gll - the nodes
187: . viewer - the viewer
189: Level: beginner
191: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()
193: @*/
194: PetscErrorCode PetscGLLView(PetscGLL *gll,PetscViewer viewer)
195: {
196: PetscErrorCode ierr;
197: PetscBool iascii;
199: PetscInt i;
202: if (!viewer) {
203: PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);
204: }
206: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
207: if (iascii) {
208: PetscViewerASCIIPrintf(viewer,"%D Gauss-Lobatto-Legendre (GLL) nodes and weights\n",gll->n);
209: for (i=0; i<gll->n; i++) {
210: PetscViewerASCIIPrintf(viewer," %D %16.14e %16.14e\n",i,(double)gll->nodes[i],(double)gll->weights[i]);
211: }
212: }
213: return(0);
214: }
216: /*@C
217: PetscGLLIntegrate - Compute the L2 integral of a function on the GLL points
219: Not Collective
221: Input Parameter:
222: + gll - the nodes
223: . f - the function values at the nodes
225: Output Parameter:
226: . in - the value of the integral
228: Level: beginner
230: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()
232: @*/
233: PetscErrorCode PetscGLLIntegrate(PetscGLL *gll,const PetscReal *f,PetscReal *in)
234: {
235: PetscInt i;
238: *in = 0.;
239: for (i=0; i<gll->n; i++) {
240: *in += f[i]*f[i]*gll->weights[i];
241: }
242: return(0);
243: }
245: /*@C
246: PetscGLLElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
248: Not Collective
250: Input Parameter:
251: . gll - the nodes
253: Output Parameter:
254: . A - the stiffness element
256: Level: beginner
258: Notes: Destroy this with PetscGLLElementLaplacianDestroy()
260: 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)
262: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
264: @*/
265: PetscErrorCode PetscGLLElementLaplacianCreate(PetscGLL *gll,PetscReal ***AA)
266: {
267: PetscReal **A;
268: PetscErrorCode ierr;
269: const PetscReal *nodes = gll->nodes;
270: const PetscInt n = gll->n, p = gll->n-1;
271: PetscReal z0,z1,z2 = -1,x,Lpj,Lpr;
272: PetscInt i,j,nn,r;
275: PetscMalloc1(n,&A);
276: PetscMalloc1(n*n,&A[0]);
277: for (i=1; i<n; i++) A[i] = A[i-1]+n;
279: for (j=1; j<p; j++) {
280: x = nodes[j];
281: z0 = 1.;
282: z1 = x;
283: for (nn=1; nn<p; nn++) {
284: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
285: z0 = z1;
286: z1 = z2;
287: }
288: Lpj=z2;
289: for (r=1; r<p; r++) {
290: if (r == j) {
291: A[j][j]=2./(3.*(1.-nodes[j]*nodes[j])*Lpj*Lpj);
292: } else {
293: x = nodes[r];
294: z0 = 1.;
295: z1 = x;
296: for (nn=1; nn<p; nn++) {
297: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
298: z0 = z1;
299: z1 = z2;
300: }
301: Lpr = z2;
302: A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(nodes[j]-nodes[r])*(nodes[j]-nodes[r]));
303: }
304: }
305: }
306: for (j=1; j<p+1; j++) {
307: x = nodes[j];
308: z0 = 1.;
309: z1 = x;
310: for (nn=1; nn<p; nn++) {
311: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
312: z0 = z1;
313: z1 = z2;
314: }
315: Lpj = z2;
316: A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+nodes[j])*(1.+nodes[j]));
317: A[0][j] = A[j][0];
318: }
319: for (j=0; j<p; j++) {
320: x = nodes[j];
321: z0 = 1.;
322: z1 = x;
323: for (nn=1; nn<p; nn++) {
324: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
325: z0 = z1;
326: z1 = z2;
327: }
328: Lpj=z2;
330: A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-nodes[j])*(1.-nodes[j]));
331: A[j][p] = A[p][j];
332: }
333: A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
334: A[p][p]=A[0][0];
335: *AA = A;
336: return(0);
337: }
339: /*@C
340: PetscGLLElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
342: Not Collective
344: Input Parameter:
345: + gll - the nodes
346: - A - the stiffness element
348: Level: beginner
350: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate()
352: @*/
353: PetscErrorCode PetscGLLElementLaplacianDestroy(PetscGLL *gll,PetscReal ***AA)
354: {
358: PetscFree((*AA)[0]);
359: PetscFree(*AA);
360: *AA = NULL;
361: return(0);
362: }
364: /*@C
365: PetscGLLElementGradientCreate - computes the gradient for a single 1d GLL element
367: Not Collective
369: Input Parameter:
370: . gll - the nodes
372: Output Parameter:
373: . AA - the stiffness element
374: - AAT - the transpose of AA (pass in NULL if you do not need this array)
376: Level: beginner
378: Notes: Destroy this with PetscGLLElementGradientDestroy()
380: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
382: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
384: @*/
385: PetscErrorCode PetscGLLElementGradientCreate(PetscGLL *gll,PetscReal ***AA, PetscReal ***AAT)
386: {
387: PetscReal **A, **AT = NULL;
388: PetscErrorCode ierr;
389: const PetscReal *nodes = gll->nodes;
390: const PetscInt n = gll->n, p = gll->n-1;
391: PetscReal q,qp,Li, Lj,d0;
392: PetscInt i,j;
395: PetscMalloc1(n,&A);
396: PetscMalloc1(n*n,&A[0]);
397: for (i=1; i<n; i++) A[i] = A[i-1]+n;
399: if (AAT) {
400: PetscMalloc1(n,&AT);
401: PetscMalloc1(n*n,&AT[0]);
402: for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
403: }
405: if (n==1) {A[0][0] = 0.;}
406: d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
407: for (i=0; i<n; i++) {
408: for (j=0; j<n; j++) {
409: A[i][j] = 0.;
410: qAndLEvaluation(p,nodes[i],&q,&qp,&Li);
411: qAndLEvaluation(p,nodes[j],&q,&qp,&Lj);
412: if (i!=j) A[i][j] = Li/(Lj*(nodes[i]-nodes[j]));
413: if ((j==i) && (i==0)) A[i][j] = -d0;
414: if (j==i && i==p) A[i][j] = d0;
415: if (AT) AT[j][i] = A[i][j];
416: }
417: }
418: if (AAT) *AAT = AT;
419: *AA = A;
420: return(0);
421: }
423: /*@C
424: PetscGLLElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGLLElementGradientCreate()
426: Not Collective
428: Input Parameter:
429: + gll - the nodes
430: . AA - the stiffness element
431: - AAT - the transpose of the element
433: Level: beginner
435: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()
437: @*/
438: PetscErrorCode PetscGLLElementGradientDestroy(PetscGLL *gll,PetscReal ***AA,PetscReal ***AAT)
439: {
443: PetscFree((*AA)[0]);
444: PetscFree(*AA);
445: *AA = NULL;
446: if (*AAT) {
447: PetscFree((*AAT)[0]);
448: PetscFree(*AAT);
449: *AAT = NULL;
450: }
451: return(0);
452: }
454: /*@C
455: PetscGLLElementAdvectionCreate - computes the advection operator for a single 1d GLL element
457: Not Collective
459: Input Parameter:
460: . gll - the nodes
462: Output Parameter:
463: . AA - the stiffness element
465: Level: beginner
467: Notes: Destroy this with PetscGLLElementAdvectionDestroy()
469: This is the same as the Gradient operator multiplied by the diagonal mass matrix
471: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
473: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
475: @*/
476: PetscErrorCode PetscGLLElementAdvectionCreate(PetscGLL *gll,PetscReal ***AA)
477: {
478: PetscReal **D;
479: PetscErrorCode ierr;
480: const PetscReal *weights = gll->weights;
481: const PetscInt n = gll->n;
482: PetscInt i,j;
485: PetscGLLElementGradientCreate(gll,&D,NULL);
486: for (i=0; i<n; i++){
487: for (j=0; j<n; j++) {
488: D[i][j] = weights[i]*D[i][j];
489: }
490: }
491: *AA = D;
492: return(0);
493: }
495: /*@C
496: PetscGLLElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
498: Not Collective
500: Input Parameter:
501: + gll - the nodes
502: - A - advection
504: Level: beginner
506: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()
508: @*/
509: PetscErrorCode PetscGLLElementAdvectionDestroy(PetscGLL *gll,PetscReal ***AA)
510: {
514: PetscFree((*AA)[0]);
515: PetscFree(*AA);
516: *AA = NULL;
517: return(0);
518: }
520: PetscErrorCode PetscGLLElementMassCreate(PetscGLL *gll,PetscReal ***AA)
521: {
522: PetscReal **A;
523: PetscErrorCode ierr;
524: const PetscReal *weights = gll->weights;
525: const PetscInt n = gll->n;
526: PetscInt i,j;
529: PetscMalloc1(n,&A);
530: PetscMalloc1(n*n,&A[0]);
531: for (i=1; i<n; i++) A[i] = A[i-1]+n;
532: if (n==1) {A[0][0] = 0.;}
533: for (i=0; i<n; i++) {
534: for (j=0; j<n; j++) {
535: A[i][j] = 0.;
536: if (j==i) A[i][j] = weights[i];
537: }
538: }
539: *AA = A;
540: return(0);
541: }
543: PetscErrorCode PetscGLLElementMassDestroy(PetscGLL *gll,PetscReal ***AA)
544: {
548: PetscFree((*AA)[0]);
549: PetscFree(*AA);
550: *AA = NULL;
551: return(0);
552: }