Actual source code: petscgll.c
petsc-3.8.4 2018-03-24
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: PetscMalloc2(n,&gll->nodes,n,&gll->weights);
76: if (type == PETSCGLL_VIA_LINEARALGEBRA) {
77: PetscReal *M,si;
78: PetscBLASInt bn,lierr;
79: PetscReal x,z0,z1,z2;
80: PetscInt i,p = n - 1,nn;
82: gll->nodes[0] =-1.0;
83: gll->nodes[n-1] = 1.0;
84: if (n-2 > 0){
85: PetscMalloc1(n-1,&M);
86: for (i=0; i<n-2; i++) {
87: si = ((PetscReal)i)+1.0;
88: M[i]=0.5*PetscSqrtReal(si*(si+2.0)/((si+0.5)*(si+1.5)));
89: }
90: PetscBLASIntCast(n-2,&bn);
91: PetscMemzero(&gll->nodes[1],bn*sizeof(gll->nodes[1]));
92: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
93: x=0;
94: PetscStackCallBLAS("LAPACKsteqr",LAPACKREALsteqr_("N",&bn,&gll->nodes[1],M,&x,&bn,M,&lierr));
95: if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in STERF Lapack routine %d",(int)lierr);
96: PetscFPTrapPop();
97: PetscFree(M);
98: }
99: if ((n-1)%2==0) {
100: gll->nodes[(n-1)/2] = 0.0; /* hard wire to exactly 0.0 since linear algebra produces nonzero */
101: }
103: gll->weights[0] = gll->weights[p] = 2.0/(((PetscReal)(p))*(((PetscReal)p)+1.0));
104: z2 = -1.; /* Dummy value to avoid -Wmaybe-initialized */
105: for (i=1; i<p; i++) {
106: x = gll->nodes[i];
107: z0 = 1.0;
108: z1 = x;
109: for (nn=1; nn<p; nn++) {
110: z2 = x*z1*(2.0*((PetscReal)nn)+1.0)/(((PetscReal)nn)+1.0)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.0));
111: z0 = z1;
112: z1 = z2;
113: }
114: gll->weights[i]=2.0/(((PetscReal)p)*(((PetscReal)p)+1.0)*z2*z2);
115: }
116: } else {
117: PetscInt j,m;
118: PetscReal z1,z,q,qp,Ln;
119: PetscReal *pt;
120: PetscMalloc1(n,&pt);
122: if (n > 30) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"PETSCGLL_VIA_NEWTON produces incorrect answers for n > 30");
123: gll->nodes[0] = -1.0;
124: gll->nodes[n-1] = 1.0;
125: gll->weights[0] = gll->weights[n-1] = 2./(((PetscReal)n)*(((PetscReal)n)-1.0));;
126: m = (n-1)/2; /* The roots are symmetric, so we only find half of them. */
127: for (j=1; j<=m; j++) { /* Loop over the desired roots. */
128: 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));
129: /* Starting with the above approximation to the ith root, we enter */
130: /* the main loop of refinement by Newton's method. */
131: do {
132: qAndLEvaluation(n-1,z,&q,&qp,&Ln);
133: z1 = z;
134: z = z1-q/qp; /* Newton's method. */
135: } while (PetscAbs(z-z1) > 10.*PETSC_MACHINE_EPSILON);
136: qAndLEvaluation(n-1,z,&q,&qp,&Ln);
138: gll->nodes[j] = z;
139: gll->nodes[n-1-j] = -z; /* and put in its symmetric counterpart. */
140: gll->weights[j] = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln); /* Compute the weight */
141: gll->weights[n-1-j] = gll->weights[j]; /* and its symmetric counterpart. */
142: pt[j]=qp;
143: }
145: if ((n-1)%2==0) {
146: qAndLEvaluation(n-1,0.0,&q,&qp,&Ln);
147: gll->nodes[(n-1)/2] = 0.0;
148: gll->weights[(n-1)/2] = 2.0/(((PetscReal)n)*(((PetscReal)n)-1.)*Ln*Ln);
149: }
150: PetscFree(pt);
151: }
152: gll->n = n;
153: return(0);
154: }
156: /*@C
157: PetscGLLDestroy - destroys a set of GLL nodes and weights
159: Not Collective
161: Input Parameter:
162: . gll - the nodes
164: Level: beginner
166: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLView()
168: @*/
169: PetscErrorCode PetscGLLDestroy(PetscGLL *gll)
170: {
174: PetscFree2(gll->nodes,gll->weights);
175: gll->n = 0;
176: return(0);
177: }
179: /*@C
180: PetscGLLView - views a set of GLL nodes
182: Not Collective
184: Input Parameter:
185: + gll - the nodes
186: . viewer - the viewer
188: Level: beginner
190: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()
192: @*/
193: PetscErrorCode PetscGLLView(PetscGLL *gll,PetscViewer viewer)
194: {
195: PetscErrorCode ierr;
196: PetscBool iascii;
198: PetscInt i;
201: if (!viewer) {
202: PetscViewerASCIIGetStdout(PETSC_COMM_SELF,&viewer);
203: }
205: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
206: if (iascii) {
207: PetscViewerASCIIPrintf(viewer,"%D Gauss-Lobatto-Legendre (GLL) nodes and weights\n",gll->n);
208: for (i=0; i<gll->n; i++) {
209: PetscViewerASCIIPrintf(viewer," %D %16.14e %16.14e\n",i,(double)gll->nodes[i],(double)gll->weights[i]);
210: }
211: }
212: return(0);
213: }
215: /*@C
216: PetscGLLIntegrate - Compute the L2 integral of a function on the GLL points
218: Not Collective
220: Input Parameter:
221: + gll - the nodes
222: . f - the function values at the nodes
224: Output Parameter:
225: . in - the value of the integral
227: Level: beginner
229: .seealso: PetscGLL, PetscGLLCreate(), PetscGLLDestroy()
231: @*/
232: PetscErrorCode PetscGLLIntegrate(PetscGLL *gll,const PetscReal *f,PetscReal *in)
233: {
234: PetscInt i;
237: *in = 0.;
238: for (i=0; i<gll->n; i++) {
239: *in += f[i]*f[i]*gll->weights[i];
240: }
241: return(0);
242: }
244: /*@C
245: PetscGLLElementLaplacianCreate - computes the Laplacian for a single 1d GLL element
247: Not Collective
249: Input Parameter:
250: . gll - the nodes
252: Output Parameter:
253: . A - the stiffness element
255: Level: beginner
257: Notes: Destroy this with PetscGLLElementLaplacianDestroy()
259: 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)
261: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
263: @*/
264: PetscErrorCode PetscGLLElementLaplacianCreate(PetscGLL *gll,PetscReal ***AA)
265: {
266: PetscReal **A;
267: PetscErrorCode ierr;
268: const PetscReal *nodes = gll->nodes;
269: const PetscInt n = gll->n, p = gll->n-1;
270: PetscReal z0,z1,z2 = 0,x,Lpj,Lpr;
271: PetscInt i,j,nn,r;
274: PetscMalloc1(n,&A);
275: PetscMalloc1(n*n,&A[0]);
276: for (i=1; i<n; i++) A[i] = A[i-1]+n;
278: for (j=1; j<p; j++) {
279: x = nodes[j];
280: z0 = 1.;
281: z1 = x;
282: for (nn=1; nn<p; nn++) {
283: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
284: z0 = z1;
285: z1 = z2;
286: }
287: Lpj=z2;
288: for (r=1; r<p; r++) {
289: if (r == j) {
290: A[j][j]=2./(3.*(1.-nodes[j]*nodes[j])*Lpj*Lpj);
291: } else {
292: x = nodes[r];
293: z0 = 1.;
294: z1 = x;
295: for (nn=1; nn<p; nn++) {
296: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
297: z0 = z1;
298: z1 = z2;
299: }
300: Lpr = z2;
301: A[r][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*Lpr*(nodes[j]-nodes[r])*(nodes[j]-nodes[r]));
302: }
303: }
304: }
305: for (j=1; j<p+1; j++) {
306: x = nodes[j];
307: z0 = 1.;
308: z1 = x;
309: for (nn=1; nn<p; nn++) {
310: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
311: z0 = z1;
312: z1 = z2;
313: }
314: Lpj = z2;
315: A[j][0] = 4.*PetscPowRealInt(-1.,p)/(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.+nodes[j])*(1.+nodes[j]));
316: A[0][j] = A[j][0];
317: }
318: for (j=0; j<p; j++) {
319: x = nodes[j];
320: z0 = 1.;
321: z1 = x;
322: for (nn=1; nn<p; nn++) {
323: z2 = x*z1*(2.*((PetscReal)nn)+1.)/(((PetscReal)nn)+1.)-z0*(((PetscReal)nn)/(((PetscReal)nn)+1.));
324: z0 = z1;
325: z1 = z2;
326: }
327: Lpj=z2;
329: A[p][j] = 4./(((PetscReal)p)*(((PetscReal)p)+1.)*Lpj*(1.-nodes[j])*(1.-nodes[j]));
330: A[j][p] = A[p][j];
331: }
332: A[0][0]=0.5+(((PetscReal)p)*(((PetscReal)p)+1.)-2.)/6.;
333: A[p][p]=A[0][0];
334: *AA = A;
335: return(0);
336: }
338: /*@C
339: PetscGLLElementLaplacianDestroy - frees the Laplacian for a single 1d GLL element
341: Not Collective
343: Input Parameter:
344: + gll - the nodes
345: - A - the stiffness element
347: Level: beginner
349: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate()
351: @*/
352: PetscErrorCode PetscGLLElementLaplacianDestroy(PetscGLL *gll,PetscReal ***AA)
353: {
357: PetscFree((*AA)[0]);
358: PetscFree(*AA);
359: *AA = NULL;
360: return(0);
361: }
363: /*@C
364: PetscGLLElementGradientCreate - computes the gradient for a single 1d GLL element
366: Not Collective
368: Input Parameter:
369: . gll - the nodes
371: Output Parameter:
372: . AA - the stiffness element
373: - AAT - the transpose of AA (pass in NULL if you do not need this array)
375: Level: beginner
377: Notes: Destroy this with PetscGLLElementGradientDestroy()
379: You can access entries in these arrays with AA[i][j] but in memory it is stored in contiguous memory, row oriented
381: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
383: @*/
384: PetscErrorCode PetscGLLElementGradientCreate(PetscGLL *gll,PetscReal ***AA, PetscReal ***AAT)
385: {
386: PetscReal **A, **AT = NULL;
387: PetscErrorCode ierr;
388: const PetscReal *nodes = gll->nodes;
389: const PetscInt n = gll->n, p = gll->n-1;
390: PetscReal q,qp,Li, Lj,d0;
391: PetscInt i,j;
394: PetscMalloc1(n,&A);
395: PetscMalloc1(n*n,&A[0]);
396: for (i=1; i<n; i++) A[i] = A[i-1]+n;
398: if (AAT) {
399: PetscMalloc1(n,&AT);
400: PetscMalloc1(n*n,&AT[0]);
401: for (i=1; i<n; i++) AT[i] = AT[i-1]+n;
402: }
404: if (n==1) {A[0][0] = 0.;}
405: d0 = (PetscReal)p*((PetscReal)p+1.)/4.;
406: for (i=0; i<n; i++) {
407: for (j=0; j<n; j++) {
408: A[i][j] = 0.;
409: qAndLEvaluation(p,nodes[i],&q,&qp,&Li);
410: qAndLEvaluation(p,nodes[j],&q,&qp,&Lj);
411: if (i!=j) A[i][j] = Li/(Lj*(nodes[i]-nodes[j]));
412: if ((j==i) && (i==0)) A[i][j] = -d0;
413: if (j==i && i==p) A[i][j] = d0;
414: if (AT) AT[j][i] = A[i][j];
415: }
416: }
417: if (AAT) *AAT = AT;
418: *AA = A;
419: return(0);
420: }
422: /*@C
423: PetscGLLElementGradientDestroy - frees the gradient for a single 1d GLL element obtained with PetscGLLElementGradientCreate()
425: Not Collective
427: Input Parameter:
428: + gll - the nodes
429: . AA - the stiffness element
430: - AAT - the transpose of the element
432: Level: beginner
434: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()
436: @*/
437: PetscErrorCode PetscGLLElementGradientDestroy(PetscGLL *gll,PetscReal ***AA,PetscReal ***AAT)
438: {
442: PetscFree((*AA)[0]);
443: PetscFree(*AA);
444: *AA = NULL;
445: if (*AAT) {
446: PetscFree((*AAT)[0]);
447: PetscFree(*AAT);
448: *AAT = NULL;
449: }
450: return(0);
451: }
453: /*@C
454: PetscGLLElementAdvectionCreate - computes the advection operator for a single 1d GLL element
456: Not Collective
458: Input Parameter:
459: . gll - the nodes
461: Output Parameter:
462: . AA - the stiffness element
464: Level: beginner
466: Notes: Destroy this with PetscGLLElementAdvectionDestroy()
468: This is the same as the Gradient operator multiplied by the diagonal mass matrix
470: You can access entries in this array with AA[i][j] but in memory it is stored in contiguous memory, row oriented
472: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianDestroy()
474: @*/
475: PetscErrorCode PetscGLLElementAdvectionCreate(PetscGLL *gll,PetscReal ***AA)
476: {
477: PetscReal **D;
478: PetscErrorCode ierr;
479: const PetscReal *weights = gll->weights;
480: const PetscInt n = gll->n;
481: PetscInt i,j;
484: PetscGLLElementGradientCreate(gll,&D,NULL);
485: for (i=0; i<n; i++){
486: for (j=0; j<n; j++) {
487: D[i][j] = weights[i]*D[i][j];
488: }
489: }
490: *AA = D;
491: return(0);
492: }
494: /*@C
495: PetscGLLElementAdvectionDestroy - frees the advection stiffness for a single 1d GLL element
497: Not Collective
499: Input Parameter:
500: + gll - the nodes
501: - A - advection
503: Level: beginner
505: .seealso: PetscGLL, PetscGLLDestroy(), PetscGLLView(), PetscGLLElementLaplacianCreate(), PetscGLLElementAdvectionCreate()
507: @*/
508: PetscErrorCode PetscGLLElementAdvectionDestroy(PetscGLL *gll,PetscReal ***AA)
509: {
513: PetscFree((*AA)[0]);
514: PetscFree(*AA);
515: *AA = NULL;
516: return(0);
517: }