Actual source code: dt.c
petsc-3.4.5 2014-06-29
1: /* Discretization tools */
3: #include <petscconf.h>
4: #if defined(PETSC_HAVE_MATHIMF_H)
5: #include <mathimf.h> /* this needs to be included before math.h */
6: #endif
8: #include <petscdt.h> /*I "petscdt.h" I*/
9: #include <petscblaslapack.h>
10: #include <petsc-private/petscimpl.h>
11: #include <petscviewer.h>
15: /*@
16: PetscDTLegendreEval - evaluate Legendre polynomial at points
18: Not Collective
20: Input Arguments:
21: + npoints - number of spatial points to evaluate at
22: . points - array of locations to evaluate at
23: . ndegree - number of basis degrees to evaluate
24: - degrees - sorted array of degrees to evaluate
26: Output Arguments:
27: + B - row-oriented basis evaluation matrix B[point*ndegree + degree] (dimension npoints*ndegrees, allocated by caller) (or NULL)
28: . D - row-oriented derivative evaluation matrix (or NULL)
29: - D2 - row-oriented second derivative evaluation matrix (or NULL)
31: Level: intermediate
33: .seealso: PetscDTGaussQuadrature()
34: @*/
35: PetscErrorCode PetscDTLegendreEval(PetscInt npoints,const PetscReal *points,PetscInt ndegree,const PetscInt *degrees,PetscReal *B,PetscReal *D,PetscReal *D2)
36: {
37: PetscInt i,maxdegree;
40: if (!npoints || !ndegree) return(0);
41: maxdegree = degrees[ndegree-1];
42: for (i=0; i<npoints; i++) {
43: PetscReal pm1,pm2,pd1,pd2,pdd1,pdd2,x;
44: PetscInt j,k;
45: x = points[i];
46: pm2 = 0;
47: pm1 = 1;
48: pd2 = 0;
49: pd1 = 0;
50: pdd2 = 0;
51: pdd1 = 0;
52: k = 0;
53: if (degrees[k] == 0) {
54: if (B) B[i*ndegree+k] = pm1;
55: if (D) D[i*ndegree+k] = pd1;
56: if (D2) D2[i*ndegree+k] = pdd1;
57: k++;
58: }
59: for (j=1; j<=maxdegree; j++,k++) {
60: PetscReal p,d,dd;
61: p = ((2*j-1)*x*pm1 - (j-1)*pm2)/j;
62: d = pd2 + (2*j-1)*pm1;
63: dd = pdd2 + (2*j-1)*pd1;
64: pm2 = pm1;
65: pm1 = p;
66: pd2 = pd1;
67: pd1 = d;
68: pdd2 = pdd1;
69: pdd1 = dd;
70: if (degrees[k] == j) {
71: if (B) B[i*ndegree+k] = p;
72: if (D) D[i*ndegree+k] = d;
73: if (D2) D2[i*ndegree+k] = dd;
74: }
75: }
76: }
77: return(0);
78: }
82: /*@
83: PetscDTGaussQuadrature - create Gauss quadrature
85: Not Collective
87: Input Arguments:
88: + npoints - number of points
89: . a - left end of interval (often-1)
90: - b - right end of interval (often +1)
92: Output Arguments:
93: + x - quadrature points
94: - w - quadrature weights
96: Level: intermediate
98: References:
99: Golub and Welsch, Calculation of Quadrature Rules, Math. Comp. 23(106), 221--230, 1969.
101: .seealso: PetscDTLegendreEval()
102: @*/
103: PetscErrorCode PetscDTGaussQuadrature(PetscInt npoints,PetscReal a,PetscReal b,PetscReal *x,PetscReal *w)
104: {
106: PetscInt i;
107: PetscReal *work;
108: PetscScalar *Z;
109: PetscBLASInt N,LDZ,info;
112: /* Set up the Golub-Welsch system */
113: for (i=0; i<npoints; i++) {
114: x[i] = 0; /* diagonal is 0 */
115: if (i) w[i-1] = 0.5 / PetscSqrtReal(1 - 1./PetscSqr(2*i));
116: }
117: PetscMalloc2(npoints*npoints,PetscScalar,&Z,PetscMax(1,2*npoints-2),PetscReal,&work);
118: PetscBLASIntCast(npoints,&N);
119: LDZ = N;
120: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
121: PetscStackCallBLAS("LAPACKsteqr",LAPACKsteqr_("I",&N,x,w,Z,&LDZ,work,&info));
122: PetscFPTrapPop();
123: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"xSTEQR error");
125: for (i=0; i<(npoints+1)/2; i++) {
126: PetscReal y = 0.5 * (-x[i] + x[npoints-i-1]); /* enforces symmetry */
127: x[i] = (a+b)/2 - y*(b-a)/2;
128: x[npoints-i-1] = (a+b)/2 + y*(b-a)/2;
130: w[i] = w[npoints-1-i] = 0.5*(b-a)*(PetscSqr(PetscAbsScalar(Z[i*npoints])) + PetscSqr(PetscAbsScalar(Z[(npoints-i-1)*npoints])));
131: }
132: PetscFree2(Z,work);
133: return(0);
134: }
138: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
139: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
140: PETSC_STATIC_INLINE PetscErrorCode PetscDTFactorial_Internal(PetscInt n, PetscReal *factorial)
141: {
142: PetscReal f = 1.0;
143: PetscInt i;
146: for (i = 1; i < n+1; ++i) f *= i;
147: *factorial = f;
148: return(0);
149: }
153: /* Evaluates the nth jacobi polynomial with weight parameters a,b at a point x.
154: Recurrence relations implemented from the pseudocode given in Karniadakis and Sherwin, Appendix B */
155: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobi(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
156: {
157: PetscReal apb, pn1, pn2;
158: PetscInt k;
161: if (!n) {*P = 1.0; return(0);}
162: if (n == 1) {*P = 0.5 * (a - b + (a + b + 2.0) * x); return(0);}
163: apb = a + b;
164: pn2 = 1.0;
165: pn1 = 0.5 * (a - b + (apb + 2.0) * x);
166: *P = 0.0;
167: for (k = 2; k < n+1; ++k) {
168: PetscReal a1 = 2.0 * k * (k + apb) * (2.0*k + apb - 2.0);
169: PetscReal a2 = (2.0 * k + apb - 1.0) * (a*a - b*b);
170: PetscReal a3 = (2.0 * k + apb - 2.0) * (2.0 * k + apb - 1.0) * (2.0 * k + apb);
171: PetscReal a4 = 2.0 * (k + a - 1.0) * (k + b - 1.0) * (2.0 * k + apb);
173: a2 = a2 / a1;
174: a3 = a3 / a1;
175: a4 = a4 / a1;
176: *P = (a2 + a3 * x) * pn1 - a4 * pn2;
177: pn2 = pn1;
178: pn1 = *P;
179: }
180: return(0);
181: }
185: /* Evaluates the first derivative of P_{n}^{a,b} at a point x. */
186: PETSC_STATIC_INLINE PetscErrorCode PetscDTComputeJacobiDerivative(PetscReal a, PetscReal b, PetscInt n, PetscReal x, PetscReal *P)
187: {
188: PetscReal nP;
192: if (!n) {*P = 0.0; return(0);}
193: PetscDTComputeJacobi(a+1, b+1, n-1, x, &nP);
194: *P = 0.5 * (a + b + n + 1) * nP;
195: return(0);
196: }
200: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
201: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapSquareToTriangle_Internal(PetscReal x, PetscReal y, PetscReal *xi, PetscReal *eta)
202: {
204: *xi = 0.5 * (1.0 + x) * (1.0 - y) - 1.0;
205: *eta = y;
206: return(0);
207: }
211: /* Maps from [-1,1]^2 to the (-1,1) reference triangle */
212: PETSC_STATIC_INLINE PetscErrorCode PetscDTMapCubeToTetrahedron_Internal(PetscReal x, PetscReal y, PetscReal z, PetscReal *xi, PetscReal *eta, PetscReal *zeta)
213: {
215: *xi = 0.25 * (1.0 + x) * (1.0 - y) * (1.0 - z) - 1.0;
216: *eta = 0.5 * (1.0 + y) * (1.0 - z) - 1.0;
217: *zeta = z;
218: return(0);
219: }
223: static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w)
224: {
225: PetscInt maxIter = 100;
226: PetscReal eps = 1.0e-8;
227: PetscReal a1, a2, a3, a4, a5, a6;
228: PetscInt k;
233: a1 = pow(2, a+b+1);
234: #if defined(PETSC_HAVE_TGAMMA)
235: a2 = tgamma(a + npoints + 1);
236: a3 = tgamma(b + npoints + 1);
237: a4 = tgamma(a + b + npoints + 1);
238: #else
239: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");
240: #endif
242: PetscDTFactorial_Internal(npoints, &a5);
243: a6 = a1 * a2 * a3 / a4 / a5;
244: /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.
245: Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */
246: for (k = 0; k < npoints; ++k) {
247: PetscReal r = -cos((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;
248: PetscInt j;
250: if (k > 0) r = 0.5 * (r + x[k-1]);
251: for (j = 0; j < maxIter; ++j) {
252: PetscReal s = 0.0, delta, f, fp;
253: PetscInt i;
255: for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);
256: PetscDTComputeJacobi(a, b, npoints, r, &f);
257: PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);
258: delta = f / (fp - f * s);
259: r = r - delta;
260: if (fabs(delta) < eps) break;
261: }
262: x[k] = r;
263: PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);
264: w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);
265: }
266: return(0);
267: }
271: /*@C
272: PetscDTGaussJacobiQuadrature - create Gauss-Jacobi quadrature for a simplex
274: Not Collective
276: Input Arguments:
277: + dim - The simplex dimension
278: . npoints - number of points
279: . a - left end of interval (often-1)
280: - b - right end of interval (often +1)
282: Output Arguments:
283: + points - quadrature points
284: - weights - quadrature weights
286: Level: intermediate
288: References:
289: Karniadakis and Sherwin.
290: FIAT
292: .seealso: PetscDTGaussQuadrature()
293: @*/
294: PetscErrorCode PetscDTGaussJacobiQuadrature(PetscInt dim, PetscInt npoints, PetscReal a, PetscReal b, PetscReal *points[], PetscReal *weights[])
295: {
296: PetscReal *px, *wx, *py, *wy, *pz, *wz, *x, *w;
297: PetscInt i, j, k;
301: if ((a != -1.0) || (b != 1.0)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must use default internal right now");
302: switch (dim) {
303: case 1:
304: PetscMalloc(npoints * sizeof(PetscReal), &x);
305: PetscMalloc(npoints * sizeof(PetscReal), &w);
306: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, x, w);
307: break;
308: case 2:
309: PetscMalloc(npoints*npoints*2 * sizeof(PetscReal), &x);
310: PetscMalloc(npoints*npoints * sizeof(PetscReal), &w);
311: PetscMalloc4(npoints,PetscReal,&px,npoints,PetscReal,&wx,npoints,PetscReal,&py,npoints,PetscReal,&wy);
312: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
313: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
314: for (i = 0; i < npoints; ++i) {
315: for (j = 0; j < npoints; ++j) {
316: PetscDTMapSquareToTriangle_Internal(px[i], py[j], &x[(i*npoints+j)*2+0], &x[(i*npoints+j)*2+1]);
317: w[i*npoints+j] = 0.5 * wx[i] * wy[j];
318: }
319: }
320: PetscFree4(px,wx,py,wy);
321: break;
322: case 3:
323: PetscMalloc(npoints*npoints*3 * sizeof(PetscReal), &x);
324: PetscMalloc(npoints*npoints * sizeof(PetscReal), &w);
325: PetscMalloc6(npoints,PetscReal,&px,npoints,PetscReal,&wx,npoints,PetscReal,&py,npoints,PetscReal,&wy,npoints,PetscReal,&pz,npoints,PetscReal,&wz);
326: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 0.0, 0.0, px, wx);
327: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 1.0, 0.0, py, wy);
328: PetscDTGaussJacobiQuadrature1D_Internal(npoints, 2.0, 0.0, pz, wz);
329: for (i = 0; i < npoints; ++i) {
330: for (j = 0; j < npoints; ++j) {
331: for (k = 0; k < npoints; ++k) {
332: PetscDTMapCubeToTetrahedron_Internal(px[i], py[j], pz[k], &x[((i*npoints+j)*npoints+k)*3+0], &x[((i*npoints+j)*npoints+k)*3+1], &x[((i*npoints+j)*npoints+k)*3+2]);
333: w[(i*npoints+j)*npoints+k] = 0.125 * wx[i] * wy[j] * wz[k];
334: }
335: }
336: }
337: PetscFree6(px,wx,py,wy,pz,wz);
338: break;
339: default:
340: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Cannot construct quadrature rule for dimension %d", dim);
341: }
342: if (points) *points = x;
343: if (weights) *weights = w;
344: return(0);
345: }
349: /* Overwrites A. Can only handle full-rank problems with m>=n
350: * A in column-major format
351: * Ainv in row-major format
352: * tau has length m
353: * worksize must be >= max(1,n)
354: */
355: static PetscErrorCode PetscDTPseudoInverseQR(PetscInt m,PetscInt mstride,PetscInt n,PetscReal *A_in,PetscReal *Ainv_out,PetscScalar *tau,PetscInt worksize,PetscScalar *work)
356: {
358: PetscBLASInt M,N,K,lda,ldb,ldwork,info;
359: PetscScalar *A,*Ainv,*R,*Q,Alpha;
362: #if defined(PETSC_USE_COMPLEX)
363: {
364: PetscInt i,j;
365: PetscMalloc2(m*n,PetscScalar,&A,m*n,PetscScalar,&Ainv);
366: for (j=0; j<n; j++) {
367: for (i=0; i<m; i++) A[i+m*j] = A_in[i+mstride*j];
368: }
369: mstride = m;
370: }
371: #else
372: A = A_in;
373: Ainv = Ainv_out;
374: #endif
376: PetscBLASIntCast(m,&M);
377: PetscBLASIntCast(n,&N);
378: PetscBLASIntCast(mstride,&lda);
379: PetscBLASIntCast(worksize,&ldwork);
380: PetscFPTrapPush(PETSC_FP_TRAP_OFF);
381: LAPACKgeqrf_(&M,&N,A,&lda,tau,work,&ldwork,&info);
382: PetscFPTrapPop();
383: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xGEQRF error");
384: R = A; /* Upper triangular part of A now contains R, the rest contains the elementary reflectors */
386: /* Extract an explicit representation of Q */
387: Q = Ainv;
388: PetscMemcpy(Q,A,mstride*n*sizeof(PetscScalar));
389: K = N; /* full rank */
390: LAPACKungqr_(&M,&N,&K,Q,&lda,tau,work,&ldwork,&info);
391: if (info) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"xORGQR/xUNGQR error");
393: /* Compute A^{-T} = (R^{-1} Q^T)^T = Q R^{-T} */
394: Alpha = 1.0;
395: ldb = lda;
396: BLAStrsm_("Right","Upper","ConjugateTranspose","NotUnitTriangular",&M,&N,&Alpha,R,&lda,Q,&ldb);
397: /* Ainv is Q, overwritten with inverse */
399: #if defined(PETSC_USE_COMPLEX)
400: {
401: PetscInt i;
402: for (i=0; i<m*n; i++) Ainv_out[i] = PetscRealPart(Ainv[i]);
403: PetscFree2(A,Ainv);
404: }
405: #endif
406: return(0);
407: }
411: /* Computes integral of L_p' over intervals {(x0,x1),(x1,x2),...} */
412: static PetscErrorCode PetscDTLegendreIntegrate(PetscInt ninterval,const PetscReal *x,PetscInt ndegree,const PetscInt *degrees,PetscBool Transpose,PetscReal *B)
413: {
415: PetscReal *Bv;
416: PetscInt i,j;
419: PetscMalloc((ninterval+1)*ndegree*sizeof(PetscReal),&Bv);
420: /* Point evaluation of L_p on all the source vertices */
421: PetscDTLegendreEval(ninterval+1,x,ndegree,degrees,Bv,NULL,NULL);
422: /* Integral over each interval: \int_a^b L_p' = L_p(b)-L_p(a) */
423: for (i=0; i<ninterval; i++) {
424: for (j=0; j<ndegree; j++) {
425: if (Transpose) B[i+ninterval*j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
426: else B[i*ndegree+j] = Bv[(i+1)*ndegree+j] - Bv[i*ndegree+j];
427: }
428: }
429: PetscFree(Bv);
430: return(0);
431: }
435: /*@
436: PetscDTReconstructPoly - create matrix representing polynomial reconstruction using cell intervals and evaluation at target intervals
438: Not Collective
440: Input Arguments:
441: + degree - degree of reconstruction polynomial
442: . nsource - number of source intervals
443: . sourcex - sorted coordinates of source cell boundaries (length nsource+1)
444: . ntarget - number of target intervals
445: - targetx - sorted coordinates of target cell boundaries (length ntarget+1)
447: Output Arguments:
448: . R - reconstruction matrix, utarget = sum_s R[t*nsource+s] * usource[s]
450: Level: advanced
452: .seealso: PetscDTLegendreEval()
453: @*/
454: PetscErrorCode PetscDTReconstructPoly(PetscInt degree,PetscInt nsource,const PetscReal *sourcex,PetscInt ntarget,const PetscReal *targetx,PetscReal *R)
455: {
457: PetscInt i,j,k,*bdegrees,worksize;
458: PetscReal xmin,xmax,center,hscale,*sourcey,*targety,*Bsource,*Bsinv,*Btarget;
459: PetscScalar *tau,*work;
465: if (degree >= nsource) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Reconstruction degree %D must be less than number of source intervals %D",degree,nsource);
466: #if defined(PETSC_USE_DEBUG)
467: for (i=0; i<nsource; i++) {
468: if (sourcex[i] >= sourcex[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Source interval %D has negative orientation (%G,%G)",i,sourcex[i],sourcex[i+1]);
469: }
470: for (i=0; i<ntarget; i++) {
471: if (targetx[i] >= targetx[i+1]) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"Target interval %D has negative orientation (%G,%G)",i,targetx[i],targetx[i+1]);
472: }
473: #endif
474: xmin = PetscMin(sourcex[0],targetx[0]);
475: xmax = PetscMax(sourcex[nsource],targetx[ntarget]);
476: center = (xmin + xmax)/2;
477: hscale = (xmax - xmin)/2;
478: worksize = nsource;
479: PetscMalloc4(degree+1,PetscInt,&bdegrees,nsource+1,PetscReal,&sourcey,nsource*(degree+1),PetscReal,&Bsource,worksize,PetscScalar,&work);
480: PetscMalloc4(nsource,PetscScalar,&tau,nsource*(degree+1),PetscReal,&Bsinv,ntarget+1,PetscReal,&targety,ntarget*(degree+1),PetscReal,&Btarget);
481: for (i=0; i<=nsource; i++) sourcey[i] = (sourcex[i]-center)/hscale;
482: for (i=0; i<=degree; i++) bdegrees[i] = i+1;
483: PetscDTLegendreIntegrate(nsource,sourcey,degree+1,bdegrees,PETSC_TRUE,Bsource);
484: PetscDTPseudoInverseQR(nsource,nsource,degree+1,Bsource,Bsinv,tau,nsource,work);
485: for (i=0; i<=ntarget; i++) targety[i] = (targetx[i]-center)/hscale;
486: PetscDTLegendreIntegrate(ntarget,targety,degree+1,bdegrees,PETSC_FALSE,Btarget);
487: for (i=0; i<ntarget; i++) {
488: PetscReal rowsum = 0;
489: for (j=0; j<nsource; j++) {
490: PetscReal sum = 0;
491: for (k=0; k<degree+1; k++) {
492: sum += Btarget[i*(degree+1)+k] * Bsinv[k*nsource+j];
493: }
494: R[i*nsource+j] = sum;
495: rowsum += sum;
496: }
497: for (j=0; j<nsource; j++) R[i*nsource+j] /= rowsum; /* normalize each row */
498: }
499: PetscFree4(bdegrees,sourcey,Bsource,work);
500: PetscFree4(tau,Bsinv,targety,Btarget);
501: return(0);
502: }