Actual source code: ex1.c
petsc-3.13.6 2020-09-29
1: static char help[] = "Tests 1D discretization tools.\n\n";
3: #include <petscdt.h>
4: #include <petscviewer.h>
5: #include <petsc/private/petscimpl.h>
6: #include <petsc/private/dtimpl.h>
8: static PetscErrorCode CheckPoints(const char *name,PetscInt npoints,const PetscReal *points,PetscInt ndegrees,const PetscInt *degrees)
9: {
11: PetscReal *B,*D,*D2;
12: PetscInt i,j;
15: PetscMalloc3(npoints*ndegrees,&B,npoints*ndegrees,&D,npoints*ndegrees,&D2);
16: PetscDTLegendreEval(npoints,points,ndegrees,degrees,B,D,D2);
17: PetscPrintf(PETSC_COMM_WORLD,"%s\n",name);
18: for (i=0; i<npoints; i++) {
19: for (j=0; j<ndegrees; j++) {
20: PetscReal b,d,d2;
21: b = B[i*ndegrees+j];
22: d = D[i*ndegrees+j];
23: d2 = D2[i*ndegrees+j];
24: if (PetscAbsReal(b) < PETSC_SMALL) b = 0;
25: if (PetscAbsReal(d) < PETSC_SMALL) d = 0;
26: if (PetscAbsReal(d2) < PETSC_SMALL) d2 = 0;
27: PetscPrintf(PETSC_COMM_WORLD,"degree %D at %12.4g: B=%12.4g D=%12.4g D2=%12.4g\n",degrees[j],(double)points[i],(double)b,(double)d,(double)d2);
28: }
29: }
30: PetscFree3(B,D,D2);
31: return(0);
32: }
34: typedef PetscErrorCode(*quadratureFunc)(PetscInt,PetscReal,PetscReal,PetscReal,PetscReal,PetscReal[],PetscReal[]);
36: static PetscErrorCode CheckQuadrature_Basics(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[])
37: {
38: PetscInt i;
41: for (i = 1; i < npoints; i++) {
42: if (x[i] <= x[i-1]) SETERRQ6(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Quadrature points not monotonically increasing, %D points, alpha = %g, beta = %g, i = %D, x[i] = %g, x[i-1] = %g\n",npoints, (double) alpha, (double) beta, i, x[i], x[i-1]);
43: }
44: for (i = 0; i < npoints; i++) {
45: if (w[i] <= 0.) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Quadrature weight not positive, %D points, alpha = %g, beta = %g, i = %D, w[i] = %g\n",npoints, (double) alpha, (double) beta, i, w[i]);
46: }
47: return(0);
48: }
50: static PetscErrorCode CheckQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, const PetscReal x[], const PetscReal w[], PetscInt nexact)
51: {
52: PetscInt i, j, k;
53: PetscReal *Pi, *Pj;
54: PetscReal eps;
58: eps = PETSC_SMALL;
59: PetscMalloc2(npoints, &Pi, npoints, &Pj);
60: for (i = 0; i <= nexact; i++) {
61: PetscDTJacobiEval(npoints, alpha, beta, x, 1, &i, Pi, NULL, NULL);
62: for (j = i; j <= nexact - i; j++) {
63: PetscReal I_quad = 0.;
64: PetscReal I_exact = 0.;
65: PetscReal err, tol;
66: PetscDTJacobiEval(npoints, alpha, beta, x, 1, &j, Pj, NULL, NULL);
68: tol = eps;
69: if (i == j) {
70: I_exact = PetscPowReal(2.0, alpha + beta + 1.) / (2.*i + alpha + beta + 1.);
71: #if defined(PETSC_HAVE_LGAMMA)
72: I_exact *= PetscExpReal(PetscLGamma(i + alpha + 1.) + PetscLGamma(i + beta + 1.) - (PetscLGamma(i + alpha + beta + 1.) + PetscLGamma(i + 1.)));
73: #else
74: {
75: PetscInt ibeta = (PetscInt) beta;
77: if ((PetscReal) ibeta != beta) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"lgamma() - math routine is unavailable.");
78: for (k = 0; k < ibeta; k++) I_exact *= (i + 1. + k) / (i + alpha + 1. + k);
79: }
80: #endif
81: tol = eps * I_exact;
82: }
83: for (k = 0; k < npoints; k++) I_quad += w[k] * (Pi[k] * Pj[k]);
84: err = PetscAbsReal(I_exact - I_quad);
85: PetscInfo7(NULL,"npoints %D, alpha %g, beta %g, i %D, j %D, exact %g, err %g\n", npoints, (double) alpha, (double) beta, i, j, (double) I_exact, (double) err);
86: if (err > tol) SETERRQ7(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Incorrectly integrated P_%D * P_%D using %D point rule with alpha = %g, beta = %g: exact %g, err %g\n", i, j, npoints, (double) alpha, (double) beta, (double) I_exact, (double) err);
87: }
88: }
89: PetscFree2(Pi, Pj);
90: return(0);
91: }
93: static PetscErrorCode CheckJacobiQuadrature(PetscInt npoints, PetscReal alpha, PetscReal beta, quadratureFunc func, PetscInt nexact)
94: {
95: PetscReal *x, *w;
99: PetscMalloc2(npoints, &x, npoints, &w);
100: (*func)(npoints, -1., 1., alpha, beta, x, w);
101: CheckQuadrature_Basics(npoints, alpha, beta, x, w);
102: CheckQuadrature(npoints, alpha, beta, x, w, nexact);
103: #if defined(PETSCDTGAUSSIANQUADRATURE_EIG)
104: /* compare methods of computing quadrature */
105: PetscDTGaussQuadratureNewton_Internal = (PetscBool) !PetscDTGaussQuadratureNewton_Internal;
106: {
107: PetscReal *x2, *w2;
108: PetscReal eps;
109: PetscInt i;
111: eps = PETSC_SMALL;
112: PetscMalloc2(npoints, &x2, npoints, &w2);
113: (*func)(npoints, -1., 1., alpha, beta, x2, w2);
114: CheckQuadrature_Basics(npoints, alpha, beta, x2, w2);
115: CheckQuadrature(npoints, alpha, beta, x2, w2, nexact);
116: for (i = 0; i < npoints; i++) {
117: PetscReal xdiff, xtol, wdiff, wtol;
119: xdiff = PetscAbsReal(x[i] - x2[i]);
120: wdiff = PetscAbsReal(w[i] - w2[i]);
121: xtol = eps * (1. + PetscMin(PetscAbsReal(x[i]),1. - PetscAbsReal(x[i])));
122: wtol = eps * (1. + w[i]);
123: PetscInfo6(NULL,"npoints %D, alpha %g, beta %g, i %D, xdiff/xtol %g, wdiff/wtol %g\n", npoints, (double) alpha, (double) beta, i, (double) xdiff/xtol, (double) wdiff/wtol);
124: if (xdiff > xtol) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatch quadrature point: %D points, alpha = %g, beta = %g, i = %D, xdiff = %g\n", npoints, (double) alpha, (double) beta, i, (double) xdiff);
125: if (wdiff > wtol) SETERRQ5(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Mismatch quadrature weight: %D points, alpha = %g, beta = %g, i = %D, wdiff = %g\n", npoints, (double) alpha, (double) beta, i, (double) wdiff);
126: }
127: PetscFree2(x2, w2);
128: }
129: /* restore */
130: PetscDTGaussQuadratureNewton_Internal = (PetscBool) !PetscDTGaussQuadratureNewton_Internal;
131: #endif
132: PetscFree2(x, w);
133: return(0);
134: }
136: int main(int argc,char **argv)
137: {
139: PetscInt degrees[1000],ndegrees,npoints,two;
140: PetscReal points[1000],weights[1000],interval[2];
141: PetscInt minpoints, maxpoints;
142: PetscBool flg;
144: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
145: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Discretization tools test options",NULL);
146: {
147: ndegrees = 1000;
148: degrees[0] = 0;
149: degrees[1] = 1;
150: degrees[2] = 2;
151: PetscOptionsIntArray("-degrees","list of degrees to evaluate","",degrees,&ndegrees,&flg);
153: if (!flg) ndegrees = 3;
154: npoints = 1000;
155: points[0] = 0.0;
156: points[1] = -0.5;
157: points[2] = 1.0;
158: PetscOptionsRealArray("-points","list of points at which to evaluate","",points,&npoints,&flg);
160: if (!flg) npoints = 3;
161: two = 2;
162: interval[0] = -1.;
163: interval[1] = 1.;
164: PetscOptionsRealArray("-interval","interval on which to construct quadrature","",interval,&two,NULL);
166: minpoints = 1;
167: PetscOptionsInt("-minpoints","minimum points for thorough Gauss-Jacobi quadrature tests","",minpoints,&minpoints,NULL);
168: maxpoints = 30;
169: #if defined(PETSC_USE_REAL_SINGLE)
170: maxpoints = 5;
171: #elif defined(PETSC_USE_REAL___FLOAT128)
172: maxpoints = 20; /* just to make test faster */
173: #endif
174: PetscOptionsInt("-maxpoints","maximum points for thorough Gauss-Jacobi quadrature tests","",maxpoints,&maxpoints,NULL);
175: }
176: PetscOptionsEnd();
177: CheckPoints("User-provided points",npoints,points,ndegrees,degrees);
179: PetscDTGaussQuadrature(npoints,interval[0],interval[1],points,weights);
180: PetscPrintf(PETSC_COMM_WORLD,"Quadrature weights\n");
181: PetscRealView(npoints,weights,PETSC_VIEWER_STDOUT_WORLD);
182: {
183: PetscReal a = interval[0],b = interval[1],zeroth,first,second;
184: PetscInt i;
185: zeroth = b - a;
186: first = (b*b - a*a)/2;
187: second = (b*b*b - a*a*a)/3;
188: for (i=0; i<npoints; i++) {
189: zeroth -= weights[i];
190: first -= weights[i] * points[i];
191: second -= weights[i] * PetscSqr(points[i]);
192: }
193: if (PetscAbs(zeroth) < 1e-10) zeroth = 0.;
194: if (PetscAbs(first) < 1e-10) first = 0.;
195: if (PetscAbs(second) < 1e-10) second = 0.;
196: PetscPrintf(PETSC_COMM_WORLD,"Moment error: zeroth=%g, first=%g, second=%g\n",(double)(-zeroth),(double)(-first),(double)(-second));
197: }
198: CheckPoints("Gauss points",npoints,points,ndegrees,degrees);
199: {
200: PetscInt i;
202: for (i = minpoints; i <= maxpoints; i++) {
203: PetscReal a1, b1, a2, b2;
205: #if defined(PETSC_HAVE_LGAMMA)
206: a1 = -0.6;
207: b1 = 1.1;
208: a2 = 2.2;
209: b2 = -0.6;
210: #else
211: a1 = 0.;
212: b1 = 1.;
213: a2 = 2.;
214: b2 = 0.;
215: #endif
216: CheckJacobiQuadrature(i, 0., 0., PetscDTGaussJacobiQuadrature, 2*i-1);
217: CheckJacobiQuadrature(i, a1, b1, PetscDTGaussJacobiQuadrature, 2*i-1);
218: CheckJacobiQuadrature(i, a2, b2, PetscDTGaussJacobiQuadrature, 2*i-1);
219: if (i >= 2) {
220: CheckJacobiQuadrature(i, 0., 0., PetscDTGaussLobattoJacobiQuadrature, 2*i-3);
221: CheckJacobiQuadrature(i, a1, b1, PetscDTGaussLobattoJacobiQuadrature, 2*i-3);
222: CheckJacobiQuadrature(i, a2, b2, PetscDTGaussLobattoJacobiQuadrature, 2*i-3);
223: }
224: }
225: }
226: PetscFinalize();
227: return ierr;
228: }
230: /*TEST
231: test:
232: suffix: 1
233: args: -degrees 1,2,3,4,5 -points 0,.2,-.5,.8,.9,1 -interval -.5,1
234: TEST*/